2024q1 Homework4 (quiz3+4)

contributed by <ShawnXuanc>

第三周測驗

測驗 1

版本一

使用 log 的方式來取得當前數 2 的次方數，即可知道 msb 為第幾個位元

版本二

以右移的方式來減少 log 的運算量

int i_sqrt(int N)
{
-   int msb = (int) log2(N);
+    int msb = 0;
+    int n = N;
+    while (n > 1) {
+        n >>= 1;
+        msb++;
+    }
+    int a = 1 << msb;
    int result = 0;
    while (a != 0) {
        if ((result + a) * (result + a) <= N)
            result += a;
        a >>= 1;
    }
    return result;
}

版本三

將

$N^2$ 以 Digit-by-digit calculation 所提及的方式進行拆解，

$N^2=(a_n+a_{n-1}+\cdots +a_0{)}^2,a_m=2^{m}or0$

乘開後的結果為

$N^2=(a_n+a_{n-1}+\cdots +a_0)(a_n+a_{n-1}+\cdots +a_0)$
用下方矩陣表示相乘之後所得結果

$\left[\begin{array}{cccc}{a}_0{a}_0& a_1{a}_0& \cdots & a_n{a}_0\\ a_0{a}_1& a_1{a}_1& \cdots & a_n{a}_1\\ a_0{a}_2& a_1{a}_2& \cdots & a_n{a}_2\\ ⋮& ⋮& \ddots & ⋮\\ a_0{a}_n& a_1{a}_n& ...& a_n{a}_n\end{array}\right]$

將其整理過後得到

$N^2=a_n^2+[2(a_n)+a_{n-1}]a_{n-1}+[2(a_n+a_{n-1})+a_{n-2}]a_{n-2}+...+[2(\sum _{i=1}^n{a}_i)+a_0]a_0$
並使用

$P_m=a_n+a_{n-1}+\cdots +a_m$ 來取代

$\sum _{i=m}^n{a}_i$的內容，得到下面式子

$N^2=a_n^2+[2P_n+a_{n-1}]a_{n-1}+[2P_{n-1}+a_{n-2}]a_{n-2}+\cdots +[2P_1+a_0]a_0$

也因為

$N=a_n+a_{n-1}+\cdots +a_0$ 因此

$P_0=a_n+a_{n-1}+\cdots +a_0$ 即為最後我們要求的結果

觀察遞迴式

$P_m=P_{m+1}+a_m$ 最後藉由遞迴式運算到

$m=0$ 得到的即我們要的答案，但是遞迴式中需要藉由

$P_m^2\le N^2$ 判斷

$P_m$ 是否要加上

$a_m$ 也就是一個

$2^m$ 的數

然而

$P_m^2\le N^2$ 的運算成本過高所以這邊藉由上一輪的結果

使用

$X_m=N^2-P_m^2$ 以及

$X_{m+1}=N^2-P_{m+1}^2$ 進行整理

1. $X_m=X_{m+1}-Y_m$
2. $Y_m=X_{m+1}-X_m=P_m^2-P_{m+1}^2=(P_{m+1}+2^m{)}^2-P_{m+1}^2$
3. $Y_m=2P_{m+1}{a}_m+a_m^2=2P_{m+1}{2}^m+(2^m{)}^2$

並將

$Y_m$ 拆解用

$c_m$ 以及

$d_m$ 表示 ,

$c_m=P_{m+1}{2}^{m+1}$,

$d_m=a_m^2$

$Y_m=$

$\{\begin{array}{l}{c}_m+d_{m}if{a}_m=2^m\\ 0if{a}_m=0\end{array}$

最後使用位元運算的方式推出下一輪的公式

$c_{m-1}=P_m{2}^m=(P_{m_{+}1}+a_m)2^m=P_{m+1}{2}^m+a_m{2}^m=\{\begin{array}{l}{c}_m/2+d_{m}if{a}_m=2^m\\ c_m/2if{a}_m=0\end{array}$

$d_{m-1}=\frac{d_m}{4}$

藉由拆解過後的公式可對應到 m 即為

$d_m$，z 為

$c_m$

int i_sqrt(int x)
{
    if (x <= 1) /* Assume x is always positive */
        return x;

    int z = 0;
    for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
        int b = z + m;
        z >>= 1;
        if (x >= b)
            x -= b, z += m;               
    }
    return z;
}

測驗 2

測驗二一開始在做的為對字串進行相加，而對其相加時需要使用到 % 以及 / 的運算，用來計算數字相加後的所得的個位數，以及使用 carry 來計算 10 是否需要進位

在過程中對於使用到的 % 以及 / 的操作會有較高的運算成本，因此藉由 bitwise 的操作來減少運算成本

因為需要使用 bitwise 進行運算但是 10 這個數沒辦法使用 2 的冪來做表示，而在進行 tmp 的計算時 tmp 的值最多只會到 19 也就是 9 + 9 + 1，經過式子的推導過後，再去以 19 尋找一個可以使用的近似值發現會介於 9.55 以及 10 之間

再以 bitwise 的式子來找到這個近似值最後得到

$\frac{128}{13}$ 為 9.84

d0 = q & 0b1;
d1 = q & 0b11;
d2 = q & 0b111;
q = ((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2) >> 7;

接著將我們進行運算的數 tmp 乘上

$\frac{13}{128}$ 來作為除以 10 的運算，也因為 13 為

$2^3+2^2+2^0$ ，所以讓

$(tmp\ast (2^3+2^2+2^0))/2^3$ 使其變成

$tmp+\frac{tmp}{2}+\frac{tmp}{2^3}$ 轉成 bitwise 的操作即為上面的(tmp >> 3) + (tmp >> 1) + tmp) 再將結果乘回 8 即為 * 13 的結果

而在過程中因為 bitwise 的操作會使得低位元的數值缺失，因此事先將其儲存起來再計算完後再加回去，並用右移 7 代表除以 128

r = tmp - (((q << 2) + q) << 1);

最後再以得到的除以 10 結果將其乘回 10，在用 tmp 減掉 10 位數的值剩下個位數來作為 mod 10 的結果

void divmod_10(uint32_t in, uint32_t *div, uint32_t *mod)
{
    uint32_t x = (in | 1) - (in >> 2); /* div = in/10 ==> div = 0.75*in/8 */
    uint32_t q = (x >> 4) + x;
    x = q;
    q = (q >> 8) + x;
    q = (q >> 8) + x;
    q = (q >> 8) + x;
    q = (q >> 8) + x;

    *div = (q >> 3);
    *mod = in - ((q & ~0x7) + (*div << 1));   
}

包裝過後的程式碼一開始將 x 為設定為 in 的值先加上 1 是為了避免除以 0 的情況發生並減掉右移 2 的結果也就是減掉自己的 1/4 所以最後的結果會只剩下 3/4 也就是乘上 0.75

而在下面多個對於 q 的計算是因為在前面先計算接近於除以 10 的部分也就是 x ，但這樣會導致精確度的缺失問題，所以利用不斷將值進行右移的方式來補回喪失掉的精確度

*div = (q >> 3);

在前面可以看到以 div = 0.75 * in / 8 的方式來得到除以 10 的結果而在前面計算了 0.75 * in 的部分在這邊補上前面提到的除以 8 讓最後得到的是除以 10 的結果

*mod = in - ((q & ~0x7) + (*div << 1));

同樣的跟一開始的方式一樣用減掉十位數的方式來計算出

測驗 3

先判斷是否有大於等於特定的值，當判斷完之後就直接右移藉由這樣的方式來減少右移的次數

static size_t ilog2(size_t i)
{
    size_t result = 0;
    while (i >= 65536) {
        result += 16;
        i >>= 16;
    }
    while (i >= 256) {
        result += 8;
        i >>= 8;
    }
    while (i >= 16) {
        result += 4;
        i >>= 4;
    }
    while (i >= 2) {
        result += 1;
        i >>= 1;
    }
    return result;
}

以 0000 0010 10000 0000 1111 0000 0000 0000 為例一開始判斷為大於 65536 的數直接右移 16 個 bit 會得到 0000 0000 0000 0000 0000 0010 1000 0000 之後會判斷為大於 256 的數右移 8 格 .... 0000 0000 0000 0010 在判斷大於等於 2 右移一格

以位移量遞減的方式大幅減少每次只右移一格所需花費的時間

__builtin_clz 其功能為會計算 most significant bit 後面有多少個 0 ，以 0001 0010 為例可以看到 msb 為由右往左數第 5 個 bit 而左邊有三個 0 即回傳 3

Returns the number of leading 0-bits in x, starting at the most significant bit position. If x is 0, the result is undefined.

int ilog32(uint32_t v)
{
    return (31 - __builtin_clz(v | 1));
}

在這邊 使用 v | 1 的原因是因為當 __builtin_clz 傳入 0 時會出現錯誤，所以當為 0 時將其變成 1 而藉由 __builtin_clz(1) 得到 0000 … 0001 左邊有 31 個 0 所以 __builtin_clz(v | 1) 所獲得的結果為 31 ，再用 31 相減就可以得到 0

而在非 0 的情況下藉由 __builtin_clz 來獲得 msb 之後的 0 藉此用 31 做相減就可得到目前最高位元 1 的所在位置

測驗 4

題目為計算指數加權移動平均 EWMA ， 意義在於當資料的使用時間越久其權重也會隨著時間降低

$S_t=\{\begin{array}{c}{Y}_0\\ \alpha Y_t+(1-\alpha )\cdot S_{t-1}\end{array}$

struct ewma {
    unsigned long internal;
    unsigned long factor;
    unsigned long weight;
};

在一開始時使用結構的方式儲存 ewma 運算時所會使用到的內容

void ewma_init(struct ewma *avg, unsigned long factor, unsigned long weight)
{
    if (!is_power_of_2(weight) || !is_power_of_2(factor))
        assert(0 && "weight and factor have to be a power of two!");

    avg->weight = ilog2(weight);
    avg->factor = ilog2(factor);
    avg->internal = 0;
}

並對結構中的內容進行初始化，為了要讓運算的效能可以最佳化所以在一開始時判斷 weight 以及 factor 是否為 2 的冪，並對 weight 以及 factor 取對數以便之後進行 bitwise 的使用

struct ewma *ewma_add(struct ewma *avg, unsigned long val)
{
    avg->internal = avg->internal
                        ? (((avg->internal << avg->weight) - avg->internal) +
                           (val << avg->factor)) >> avg->weight
                        : (val << avg->factor);
    return avg;
}

在計算這邊先以 avg->internal 的值做判斷，如果不為 0 就使用公式更新，若為 0 就進行初始值的賦予，但這邊會發現題目所使用的公式跟上面給予的參考公式有所出入

所以試著將公式進行整理，步驟如下

1. $internal\ast 2^{weight}-internal+val\ast factor/2^{weight}$
2. $internal(2^{weight}-1)+val\ast factor/2^{weight}$
3. $internal(1-\frac{1}{2^{weight}})2^{weight}+val\ast factor\ast \frac{1}{2^{weight}}$
4. $\frac{1}{2^{weight}}\cdot (val\ast factor)+(1-\frac{1}{2^{weight}})\cdot internal\ast 2^{weight}$

到了這邊可以發現將式子整理之後會變得跟前面是一樣的，可以得到

$\alpha$ 為

$\frac{1}{2^{weight}}$

測驗 5

int ceil_ilog2(uint32_t x)
{
    uint32_t r, shift;

    x--;
    r = (x > 0xFFFF) << 4;              // 1 << 4 or 0 => r += 16 or 0  
    x >>= r;                            // x = x >> 16 or 0
    shift = (x > 0xFF) << 3;            // 1 << 8 or 0
    x >>= shift;                        // x = x >> 8 or 0
    r |= shift;                         // r += 8 or 0
    shift = (x > 0xF) << 2;             // 1 << 4 or 0
    x >>= shift;                        // x >> 4 or 0
    r |= shift;                         // r += 4 or 0
    shift = (x > 0x3) << 1;             // 1 << 1 or 0
    x >>= shift;                        // x >> 1 or 0
    return (r | shift | x > GGG) + 1;   // r += 1 and rounding up
}

將 ceil 以及 log2 的功能進行結合

對於 r 一開始的計算藉由是否大於 65535 來獲得 1 並藉由將 1 像左移 4 來獲得 16 以這樣的方式就如同測驗 3 的版本來進行數值的判斷，進行較大幅度的位移並將 x 向右 16 個 bit

而下方的 shift 也是用同樣的方式進行位移的判斷，並且在過程中都使用 | 的方式來作為結果的計算

以最初的 r 為例，若一開始 x > 0xFFFF 成立則將 r + 16 ， x > 0xFF 則將 r + 8 以此類推，以這樣的方式計算最後要回傳的結果並切在最後進行無條件進位

第四週測驗

測驗1

unsigned popcount_naive(unsigned v)
{
    unsigned n = 0;
    while (v)
        v &= (v - 1), n = -(~n);
    return n;
}

在一開始 popcount_naive 的實作上使用的是將 v 以及 v - 1 進行 & 操作皆由這樣的方式來將最右邊的 bit 也就 lsb 的部分消除，並使用 n = -(~n) 的方式來進行 n++ ，為什麼這樣可以實現 n++ 即藉由

$-n=\sim n+1$ 進行移位成為

$-(\sim n)=n+1$ 來實現，但這樣的缺點是會需要線性時間的複雜度來計算 lsb 的數量

unsigned popcount_branchless(unsigned v)
{
    unsigned n;
    n = (v >> 1) & 0x77777777;
    v -= n;
    n = (n >> 1) & 0x77777777;
    v -= n;
    n = (n >> 1) & 0x77777777;
    v -= n;

    v = (v + (v >> 4)) & 0x0F0F0F0F;
    v *= 0x01010101;                                     

    return v >> 24;
}

對於 branchless 的方式則是使用數學的方式改寫，以每 4 個位元為一組進行計算

$popcount(x)=x-\lfloor \frac{x}{2}\rfloor -\lfloor \frac{x}{4}\rfloor -\lfloor \frac{x}{8}\rfloor$

n = (v >> 1) & 0x77777777;
v -= n;
n = (n >> 1) & 0x77777777;
v -= n;
n = (n >> 1) & 0x77777777;
v -= n;

以 x >> 1 來表示

$\lfloor \frac{x}{2}\rfloor$ 的計算以此類推，並搭配 0x77777777 來將原先右移之前的最低位元去除，以這樣的方式才不會造成在計算每 4 個位元為一組中的 1 的個數時因為右移導致左邊集合的 1 跑了過來，以正確的取得 4 個位元為一組的集合中所要的位元個數

最後會得到

$Bn$ 如下，

$Bn$ 代表原先 4 個 bit 中 1 的個數

B7 B6 B5 B4 B3 B2 B1 B0

再利用 v + x >> 4 & 0x0F0F0F0F 的方式來漂亮的得到集合由右往左 0, 2, 4, 6 為下方的結果

0 (B7+B6) 0 (B5+B4) 0 (B3+B2) 0 (B1+B0)

並運用 v *= 0x01010101 的技巧來將所有結果相加得到 B0 + B1 + B2 + B3 + B4 + B5 + B6 + B7 也就是一開始 v 中有多少個位元為 1

__builtin_popcount(x) 即為回傳 x 中位元 1 的個數

Returns the number of 1-bits in x.

Hamming distance 為兩個數字中 2 進位不同位元的個數

即可以用互斥或的方式來取的兩個數不同的位元數有幾個並用__builtin_popcount(x)來計算以達到 branchless 的辦法否則會需要一個位元一個位元的查看需要線性的時間

int totalHammingDistance(int* nums, int numsSize)
{
    int total = 0;;
    for (int i = 0;i < numsSize;i++)
        for (int j = 0; j < numsSize;j++)
            total += __builtin_popcount(nums[i] ^ nums[j]); 
    return total >> 1;
}

在這題中可以看到題目是在計算 Hamming distance ，也就是給定一個陣列，將陣列中的數每兩個數之間的 Hamming distance 計算出來並相加，在計算過程中使用互斥或的方式來得到兩個數的位元差也就是有多少個位元是不同的，在用前面提到的 __builtin_popcount 來計算 1 的個數

但這邊用的是兩個 for 迴圈的方式，以這樣來去走訪每個數字會導致有重複計算的問題，所以說最後 total 所計算出來的結果會是正確答案的兩倍，因此在最後補上 >> 1 來將結果除以 2 ，以去除掉重複計算的部分

測驗 2

根據 Hacker's Delight 中 10-18 的內容，以定理

$a+c\equiv b+d(modm)$ 以及

$ac\equiv bd(modm)$ 以 mod 3 為例子，當以

$2^k$ 的數進行 mod 時 k 為偶數時會得到

$2^k\equiv 1$ ， 而 k 為奇數時

$2^k\equiv -1$

因此在

$mod3$ 之下

$n=..{.2}^3\cdot b_3+2^2\cdot b_2+2^1\cdot b_1+2^0\cdot b_0\equiv b_k...-b_3+b_2-b_1+b_0$ 當 k 奇數時即為負的

所以以 population count 的方式來寫成 n = pop(n & 0x55555555) - pop(n &0xAAAAAAAA) 再根據文中公式 (34)

$pop(x\mathrm{\&}\bar{m})–pop(x\mathrm{\&}m)=pop(x\oplus m)-pop(m)$

將上述 n = pop(n & 0x55555555) - pop(n &0xAAAAAAAA) 轉換成 pop(n ^ 0xAAAAAAAA) - 16 藉由這樣的方式可以讓 n = 0~2 之間但是當 n 為負數時需要將 n 加上一個足夠大的 3 的倍數做運算文中使用39，而使用這樣的方式需要花費 11 個指令來作運算

int remu3(unsigned n) {
 n = pop(n ^ 0xAAAAAAAA) + 23; // Now 23 <= n <= 55.
 n = pop(n ^ 0x2A) - 3; // Now -3 <= n <= 2.
 return n + (((int)n >> 31) & 3);
}

經過加上 39 在轉換會得到 n 介於 -3 以及 2 之間，所以用 (n >> 31) & 3 即 n 如果是負的就加上 3 將其轉換成正數

最後使用查表的方式來節省指令的花費

int mod3(unsigned n) {
    static char table[33] = {2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1 };
    n = popcount(n ^ 0xAAAAAAAA);
    return table[n];
}

測驗 3

xtree 為一棵可以自平衡的二元搜尋樹，平衡方式為使用類似 AVL tree 的高度差的方式，以 hint 來作為平衡的標準，藉由 hint 判斷 xt_update 更新的方式

xt_balance 計算左右節點的 hint 的差，而 hint 代表著從該節點的子節點開始計算的最長的樹高，所以在計算時需補上自身的高度

在 xt_update 時藉由 b 來判斷當前節點的左子右子 hint 的差，將當前節點的 hint 值儲存下來，並藉由 b 的值來判斷是否要進行左旋，或是右旋，在最後判斷當前節點 hint 值是否與原先相同或是否為 0 來決定要不要對親代節點進行更新

xt_max_hint 判斷左右子 hint 較大者

xt_insert 進行節點插入，一開始會先進行節點的搜尋，藉由 _xt_find 判斷節點是否存在，若不存在則尋找等等需要插入的正確位置 xt_dir 以及其親代節點 p ，若存在則退出，在使用 __xt_insert 進行節點插入，並使用 xt_update 更新狀態

xt_rotate_left, xt_rotate_left 用來進行左右旋轉的操作

__xt_remove 進行節點的刪除的操作，將刪除的節點以右子樹進行不斷向左查找的直到左子為空的節點取代原先刪除節點，並對取代節點的右子節點進行更新，若缺少右子節點，則對刪除節點的左子樹進行的向右查找並用其取代刪除節點，最後更新取代節點的左子，若在左右子皆缺少的情況下，則直接將其移除，並更新刪除節點的親代節點