# 2024q1 Homework4 (quiz3+4) contributed by < [`Wufangni`]("https://github.com/Wufangni") > ## 第三週測驗 ### 測驗 1 `i_sqrt` 為找出平方根的函式,下面為版本一的作法: ``` c #include <math.h> int i_sqrt(int N) { int msb = (int) log2(N); int a = 1 << msb; int result = 0; while (a != 0) { if ((result + a) * (result + a) <= N) result += a; a >>= 1; } return result; } ``` 使用 `log2` 找出最高位元數 `msb` ,接著利用 `if ((result + a) * (result + a) <= N)` 判別本次迴圈的結果取平方是否小於原數,若不是則持續往下找,直到找出平方可小於原數的結果即為他的平方根。 版本二使用迴圈方式替代 `log2` 函式找出最高位元 `msb`: ``` c int i_sqrt(int 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; } ``` ``` c 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 ``` c carry = tmp / 10; tmp = tmp % 10; ``` 因 `mod` 的計算量較大,根據除法原理利用已計算出來的 `tmp / 10` 可將 `%` 改寫成下面形式: ``` c carry = tmp / 10; tmp = tmp - carry * 10; ``` 接著著手改寫 `divmod_10` 方法: ``` c #include <stdint.h> 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)); } ``` `(in | 1) - (in >> 2)` 可看成 $in-({in \over 4})$ ,化簡開可得 $0.75in$,因此 `q` 的計算結果為 ${3 \over 4}in \cdot {1 \over 16} + {3 \over 4}in$, 再整理後可得 ${102 \over 128}in$,做了四次 `q = (q >> 8) + x;` 讓 `q` 的結果更趨近於 $0.8in$ , 最後可計算 `div` 結果,${8 \over 10}in \cdot {1 \over 8}={1 \over 10}in$。 `mod` 計算則利用 `(q & ~0x7)` 方式,將 `~0x7` 看待成 `0x11...1000` 與 `q` 做 `&` 計算,相等於 `div` 乘上 $8$ ,再與 `(*div << 1)` 相加可得結果為 $div \cdot 8 + div \cdot 2 = div \cdot 10$,根據上述提到的 `mod` 簡化算法 `tmp = tmp - carry * 10` 代入數值,`mod` 即為 $in - div \cdot 10$ 結果。 ### 測驗 3 定義一個 `ilog2` 用來找出以 2 為底的對數,每次迴圈中 i 的二進位值都會右移一格,用 `log` 作為計數器計算有幾個 $2^n, n\ge0$ 。 ``` c int ilog2(int i) { int log = -1; while (i) { i >>= 1; log++; } return log; } ``` 上述方法因每個 `bit` 都需掃過一次做計算,因此改良成下面版本,在迴圈中新增判別式判斷 i 是否大於某個 $2^n, n\ge0$ ,若達成條件則一次移動 n 格 `bit` ,以減少總迴圈次數。 ``` c 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; } ``` 使用 GNU 中提供的 `__builtin_clz` ,因輸入不能為零所以設置 `v | 1` 作為輸入。 ``` c int ilog32(uint32_t v) { return (31 - __builtin_clz(v | 1)); } ``` ### 測驗 4 預先定義 `ewma` 結構, `internal` 為當前 `ewma` 的值, `factor` 為縮放因子,`weight` 用來決定衰減的 $\alpha$ 大小。 ``` c struct ewma { unsigned long internal; unsigned long factor; unsigned long weight; }; ``` 初始化 `ewma` 結構體,因 `factor` 及 `weight` 必須為 2 的冪,使用 `is_power_of_2` 排除可能報錯的值,接著使用測驗 3 提及的 `ilog2` 取出 2 的次方。 ``` c 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; } ``` 創建 `*ewma_add` 更新 `ewma` ,添加新資料 `val` ,若 `avg->internal` 為空值代表 `val` 為 $Y_0$,直接填入 `val << avg->factor` ,其餘情況則按照 $\alpha Y_t \ + \ (1-\alpha)\cdot S_{t-1}$ 更新 $S_{t}$。 $\alpha$ 以$\cfrac{1}{2^{avg->weight}}$ 表示,$\alpha Y_t$ 則可以以 `(val << avg->factor)) >> avg->weight` 實現。 $(1 -\alpha)$ 可寫成 $\cfrac{2^{avg->weight}-1} {2^{avg->weight}}$,乘回 `avg->internal` 展開可得 $[(2^{avg->weight})(avg\to internal)-(avg\to internal)]\ /\ 2^{avg->weight}$,前半部分可以 `((avg->internal << avg->weight) - avg->internal)` 實現,最後因 `(val << avg->factor)) >> avg->weight` 也需 `>> avg->weight` 因此將此移至算式尾端當分母。 ``` c 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; } ``` ## 第四周測驗題 ### 測驗 1 計算數值的二進位表示法中有被設置的 `bit` 數量, `popcount_naive` 為其中一種方式: ``` c unsigned popcount_naive(unsigned v) { unsigned n = 0; while (v) v &= (v - 1), n = -(~n); return n; } ``` 利用 `v &= (v - 1)` 重置數值中 `LSB` 位元,並使用 `n = -(~n)` 方式對 `n` 加一,直到所有 `LSB` 位元重置為止。 以另一種方式做改寫: ``` c 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; } ``` 以 32bit 舉例,population count 可以下列數學式表示: $popcount(x)=x-\lfloor \cfrac{x}{2} \rfloor - \lfloor \cfrac{x}{4}\rfloor-...-\lfloor\cfrac{x}{32} \rfloor$ 為了避免每個 bit 都做一次計算,可以以每 4 個 bit 同時做計算,即:$popcount(x)=x-\lfloor \cfrac{x}{2} \rfloor - \lfloor \cfrac{x}{4}\rfloor-\lfloor\cfrac{x}{8} \rfloor$ `n = (v >> 1) & 0x77777777` 得出$\lfloor \cfrac{x}{2} \rfloor$,持續累加三次得到以 4 bits 為單位的 `set bit` 數量。 ``` c v = (v + (v >> 4)) & 0x0F0F0F0F; v *= 0x01010101; return v >> 24; ``` 利用一次位移一個 4 bits 單位做相加乘上 `0x0F0F0F0F` 過濾多餘的內容,可得以下數列, $Bn$ 為第 n 個單位 `set bit` 數量。 ```c 0 (B7+B6) 0 (B5+B4) 0 (B3+B2) 0 (B1+B0) ``` 最後乘上 `0x01010101` 並取出在原 A6 位置的結果 (B7+B6+B5+B4+B3+B2+B1+B0),即為總 `set bit` 數目。 ``` c int hammingDistance(int x, int y) { return __builtin_popcount(x ^ y); } ``` Hamming distance 為兩數列不相同位元之數目,透過 GCC 內建函式庫中提供 `__builtin_popcount(x)` 可計算 `set bit` 數量,可利用 $XOR$ 特性(相同為0,不同為1)計算 $A\oplus B$ ,得出數列 A 與數列 B 相異位元的數列,即為 `A ^ B`。 ``` c 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; } ``` ### 測驗 3 從 `xtree.h` 中先觀察 `xt_node` 結構: ``` c struct xt_node { short hint; struct xt_node *parent; struct xt_node *left, *right; }; ``` `*parent` 作為指向父節點的指標,`*left` 及 `*right` 分別指向左右子樹的根節點,`hint` 則作為判斷是否處於平衡的變數。 利用 `treeint_xt.c` 中提供對 `xt_tree` 的基本操作,完成 `__xt_remove` 及 `xt_insert` 等節點操作。 ``` c if (xt_right(del)) { struct xt_node *least = xt_first(xt_right(del)); if (del == *root) *root = least; xt_replace_right(del, least); xt_update(root, xt_right(least)); return; } ```