# 2024q1 Homework4 (quiz3+4) contributed by < [ChenFuhuangKye](https://github.com/kyehuang) > ## 2024q1 第 3 週測驗題 ### 測驗 `1` #### 根據 [Digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) 提出直式開方法的變形 假設欲求的平方根的 $N^2$ 拆成 \begin{aligned} N^2 = (a_n + a_{n-1} + a_{n-2} + ... + a_0)^2 \\ a_m = 2^m\quad or \quad a_m = 0 \end{aligned} 將其乘開 \begin{aligned} N^2 = a_n^2 + [2P_n+a_{n-1}]a_{n-1} + [2P_{n-1}+a_{n-2}]a_{n-2} + ... + [2P_1+a_0]a_0 \\ P_m = a_n+a_{n-1}+a_{n-2}+...+a_m \end{aligned} 而 $P_0 = a_n+a_{n-1}+a_{n-2}+...+a_0$ 即為所求的平方根 $N^2$ 可以知道 \begin{aligned} P_m = P_{m+1}+a_m \end{aligned} 我們需要判斷 $a_m$ 是 $2^m$ 或是 $0$。 $a_m$ 可以透過試驗 $P_m^2 <= N^2$ 是否成立,然而這個做法成本太高。 另一種想法就是從上一輪的差值 $X_{m+1}$ 減去 $Y_m$ 得到 $X_m$ * $X_m = N^2 - P_m^2 = X_{m+1} - Y_m$ * $Y_m = P_m^2 -P_{m+1}^2=(P_{m+1}+a_m)^2-P_{m+1}^2=2P_{m+1}a_m+a_m^2$ 接著將 $Y_m$ 拆分成 $c_m$ 及 $d_m$ * $c_m=2P_{m+1}a_m = P_{m+1}2^{m+1}$ * $d_m=a_m^2=(2^m)^2$ * $Y_m = \begin{cases} c_m+d_m & if \ a_m = 2^m \\ 0 & if \ a_m = 0 \end{cases}$ 可以發現 $c_0=P_02^1 = 2P_0$ , 因此求出 $c_0/2$ 即為所求的平方根 $N^2$。 而 * $c_{m-1}=P_m2^m=(P_{m+1}+a_m)2^m=P_{m+1}2^m+a_m2^m = \begin{cases} c_m/2+dm & if \ a_m =2^m\\ c_m/2 & if \ a_m=0 \end{cases}$ * $d_{m-1} = d_m/4$ 可以知道每一輪 $c$ 需要除以2 且 $d$ 需要除以4。 因此程式碼可以改寫成 ```c int i_sqrt(int x) { if (x <= 1) /* Assume x is always positive */ return x; int c = 0; for (int d = 1UL << ((31 - __builtin_clz(x)) & ~1UL); d; d >>= 2) { int y = c + d; c >>= 1; if (x >= y) x -= y, c += d; } return c; } ``` __builtin_clz(x) 是求出 x 的最高位之前的 0 的個數,跟 ffs (find first bit in word) 功能一樣,因此可以將程式碼改成。 ```diff - for (int d = 1UL << ((31 - __builtin_clz(x)) & ~1UL); d; d >>= 2) { + for (int d = 1UL << ((31 - __ffs(x)) & ~1UL); d; d >>= 2) { ``` ### 測驗 `2` 參照 [Hacker's Delight](https://web.archive.org/web/20180517023231/http://www.hackersdelight.org/divcMore.pdf) 以及 [公式推導](https://hackmd.io/@sysprog/linux2024-quiz3#%E6%B8%AC%E9%A9%97-2) 可以知道只要找到一個除數 $x$ 且 $9.55<=x<=10$ 。假設 $n=ab$ ( a 是十位數、 b 是個位數), $n$ 除以 $x$ 的商即為 $n$ 除以 10 的商,如此一來可以透過 bitwise operation 取得商跟餘數。 此外也可以透過 $0.8n$ 除以 $8$ 得到 $0.1n$ , 而 $0.1n$ 既為 div 10 的結果,因此找到一個 $x$ , 使得 $x*n\cong0.8n$。 ```c uint32_t q = (x >> 4) + x; ``` $q=\frac{3}{4*16}n + \frac{3}{4}n = \frac{51}{64}n =0.796875n$ , $x=0.796875$ 透過 4 輪的逼近。 ```c q = (q >> 8) + x; ``` | 第幾輪 | x | |:------:|:------------------ | | 1 | 0.79998779296875 | | 2 | 0.7999999523162842 | | 3 | 0.7999999998137355 | | 4 | 0.7999999999992724 | 我們可以得到 $q\cong0.8n$ ,將 $q$ 除以 $8$ 可以得到 $div$。 $div = \frac{q}{8} \cong \frac{0.8n}{8} \cong 0.1n$ ```c (q & ~0x7) ``` 清除 q 的後三位元,相當於 div 乘 8,再加上 div 乘 2 , 可以算出 $mod=n-(8*div+2*div)=n-10*div$。 因此程式可以改寫成 ```diff 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 >> CCCC); - *mod = in - ((q & ~0x7) + (*div << DDDD)); + *div = (q >> 3); + *mod = in - ((q & ~0x7) + (*div << 1)); } ``` ```c uint32_t x = (in | 1) - (in >> 2); ``` 要避免出現 20 的倍數,因為 20 乘上 x 會出現接近 16 的值 15。 當 $q=15$ , 20 的 div 為 1 。 所以將每一個數字轉換成奇數,避免出現 20 的倍數。 ### 測驗 `3` 版本一的作法透過不斷把目標整數 i 除以 2 直到其變為 0 , 同時使用一個計數器 log 來記錄右移的次數, 這個過程是計算 $log_2(i)$ 的整數部分。因此目標整數出現第一個 1 的位數,對應 $log_2(i)$ 的整數部分。 版本二的做法是確認目標整數 i 第一個 1 的位數。 result 是一個計數器計算目前右移的次數,假設有一個整數 x 且 x 是 2 的冪,如果 i 大於 x 時,將 result 加上 $log2(x)$ 並把 i 右移 $log2(x)$ 次。 ```diff static size_t ilog2(size_t i) { size_t result = 0; - while (i >= AAAA) { + while (i >= 65536) { result += 16; i >>= 16; } - while (i >= BBBB) { + while (i >= 256) { result += 8; i >>= 8; } - while (i >= CCCC) { + while (i >= 16) { result += 4; i >>= 4; } while (i >= 2) { result += 1; i >>= 1; } return result; } ``` 版本三直接利用 `__builtin_clz` 直接取得目標整數 i 第一個 1 的位數,為了處理 i 為 0 的情況, 將 i 與 1 進行 `i | 1` ,確保 i 不為 0 。 ### 測驗 `4` 程式碼的資料型態是 `unsigned long` ,然而平均值會是小數,這邊是利用定點數,自行定義小數點的位置,如此一來就可以利用 `bitwise` 來計算小數。 程式碼中透過 `ewma_read` 函數讀取 `EWMA` , 發現讀取時會先把 S 右移 `factor` 次 , 因此在 `ewma_add` 所有新加入的值都會先左移 `factor` 次 。 因此在 `ewma_add` 中 `EWMA` 的公式可以看成 $f*S_t =\begin{cases} f*Y_0 & t=0\\ f*\alpha Y_t+(1-\alpha)*S_{t-1} & t>0\end{cases}$ >$f$ 為 $2^{factor}$ 、 $\alpha$ 為 $2^{-weight}$。 當 $t=0$ 時 , 將 `val` 左移 `factor` 次即可。 當 $t>0$ 時 , 這邊先將 $\alpha$ 提出去 $S_t=f*\alpha Y_t+(1-\alpha)*S_{t-1} = \alpha [f*Y_t + (\frac{1}{\alpha} -1)S_{t-1}] = \alpha [f*Y_t + (\frac{1}{\alpha}*S_{t-1} -S_{t-1})]$ * $\alpha=2^{-weight}$ , 代表右移`weight` 次 * $\frac1\alpha={2^{weight}}$ , 代表左移 `weight` 次 ```diff struct ewma *ewma_add(struct ewma *avg, unsigned long val) { avg->internal = avg->internal - ? (((avg->internal << EEEE) - avg->internal) + - (val << FFFF)) >> avg->weight + ? (((avg->internal << avg->weight) - avg->internal) + + (val << avg->factor)) >> avg->weight : (val << avg->factor); return avg; } ``` ### 測驗 `5` 程式碼中 `r` 、 `shift` 和 `x>1` 可以視為 `ilog2` 計算 $log_2(i)$ 的整數部分的結果 , [ceil](https://man7.org/linux/man-pages/man3/ceil.3.html) 為找到不小於 x 的整數值,如何利用 `ilog2` 的結果,計算出 $[log_2(i)]$ 需要考慮以下兩點。 第一點,如果 $log_2(i)$ 有小數 , 將 `ilog2(i)` 的結果加 1 即為答案。 第二點,i 為 2 的冪也就是 $log_2(i)$ 只有整數,不能直接將 `ilog2` 的結果加 1 ,因此透過改變 i 的值,先將 i 減 1 ,再將`ilog2(i)` 的結果加 1 即為答案。 ## 2024q1 第 4 週測驗題 ### 測驗 `1` 假設有三個數字 $A$ $B$ $C$ 。 $total = (A\bigoplus A+A\bigoplus B+A\bigoplus C)+(B\bigoplus A+B\bigoplus B+B\bigoplus C)+(C\bigoplus A+C\bigoplus B+C\bigoplus C)$ 因為 $A\bigoplus A =0$ $total = (A\bigoplus B+A\bigoplus C)+(B\bigoplus A+B\bigoplus C)+(C\bigoplus A+C\bigoplus B)$ $total = 2*(A\bigoplus B+A\bigoplus C+B\bigoplus C)$ 所以須將 total 除以2。 ```diff 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 >> AAAA; + return total >> 1; } ``` ### 測驗 `2` `popcount(n ^ 0xAAAAAAAA) - 16` 的結果介於 -16~16 ,為了避免負數,選擇加上足夠大的 3 的倍數 39 , 改成 `n = popcount(n ^ 0xAAAAAAAA) + 23` ,此時的 n 介於 23~55 ,需要再進行一次化簡。 55 出現第一個 1 的位數為 5 。 因此本次的 m 要改成 `0x2AA` , `n = popcount(n ^ 0x2A) - 3` 的結果介於 -2~2 , 當結果是負數需要加 3 。因此 `return n + ((n >> 31) & 3)` 有誤,應改成 `return n + (((int)n >> 31) & 3)`。 #### 棋盤解析 在棋盤中,有 8 種不同的連線方式可以構成勝利條件,這包括 3 行、 3 列以及 2 個對角線的連線,因此每 4 個位元 (nibble) 可以作為一的單位來判斷棋盤上是否成功形成連線。 ```c static const uint32_t move_masks[9] = { 0x40040040, 0x20004000, 0x10000404, 0x04020000, 0x02002022, 0x01000200, 0x00410001, 0x00201000, 0x00100110, }; ``` 假設第一行成功形成連線,也就是第 0 個、第 1 個以及第 2 個元素做 or 運算,可以得到 `0x70044444` ,假設第三行成功形成連線,也就是第 6 個、第 7 個以及第 8 個元素做 or 運算,可以得到 `0x00711111`。 因此可以看出如果有一組 nibble 的和為 7 代表連線成功,所以將每一組 nibble 加 1 ,跟 8 作 and 運算時就不會是 0。 ```diff - return (player_board + BBBB) & 0x88888888; + return (player_board + 0x11111111) & 0x88888888; ``` [Hacker's Delight](https://web.archive.org/web/20180517023231/http://www.hackersdelight.org/divcMore.pdf) 中有提及 mod 7 的作法。 ```c int remu7(unsigned n) { static char table[75] = {0,1,2,3,4,5,6, 0,1,2,3,4,5,6, 0,1,2,3,4,5,6, 0,1,2,3,4,5,6, 0,1,2,3,4,5,6, 0,1,2,3,4,5,6, 0,1,2,3,4,5,6, 0,1,2,3,4,5,6, 0,1,2,3,4,5,6, 0,1,2,3,4,5,6, 0,1,2,3,4}; n = (n >> 15) + (n & 0x7FFF); // Max 0x27FFE. n = (n >> 9) + (n & 0x001FF); // Max 0x33D. n = (n >> 6) + (n & 0x0003F); // Max 0x4A. return table[n]; } ``` 可以發現 `n = (n >> 15) + (n & 0x7FFF)` 為作一次化簡的步驟,但詳細的推導過程還在釐清。 ```diff - x = (x >> CCCC) + (x & UINT32_C(0x7FFF)); + x = (x >> 15) + (x & UINT32_C(0x7FFF)); ``` ### 測驗 `3`