# 2024q1 Homework 4 (quiz3+4) **contributed by <`MiohitoKiri5474`>** ## 第三週 ### 測驗一 根據 [digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) 原理,假設欲求平方根為 $N$,那我們可以把 $N^2$ 化為 2 的冪總和: $N^2 = (a_0+a_1+ \dots + a_n)^2, a_m = 2^m(or \ 2^m = 0 ), m \in \{0, 1, \dots, n \}$ 已知 $( a + b ) ^ 2 = a^2 + 2ab + b^2$,則有以下遞迴關係: $$N ^ 2 = ((a_0 + a_1 + \dots + a_{n - 1}) + a_n)^2 = {a_n}^2 +2a_n \cdot \sum^{n - 1}_{k = 0}a_k+(a_0 + a_1 + \dots + a_{n - 1})^2$$ $$N ^ 2= {a_n}^2 +(2P_n + a_{n - 1} + \dots + (2P_1 + a_0)a_0$$ 同時令平方根逼近值 $P_m = P_m + a_{m - 1}$,因此對於每個 $m \in \{0, 1, \dots, n \}$,於 $P_{m - 1}$ 加上 $a_m$,由於是 $2$ 的冪,因此只會有以下兩種可能: $$P_{m - 1} = \begin{cases} P_m + a_{m - 1}, & \text{if} \ {P_m}^2 < N^2 \\ P_m , & \text{otherwise} \end{cases}$$ 但我們的程式碼不會用這種比較直白的方式判斷是否加上 $a_{m - 1}$,而是利用差值的方式: 令 $X_{n + 1} = N^2$(考慮從 $n$ 開始),$Y_m = (2P_m + a_{m - 1})a_{m - 1}$ 而對於每個 $m$:$X_m = X_{m + 1} - Y_m = N^2 - {P_m}^2$ 但到了這邊,程式碼到這邊又進一步將 $Y_m$ 拆解成 $c_m, d_m$,令人有些費解,再經過一陣子思考後,發現只是單純將加法的兩項拆開來看,也就是說: $$\begin{cases} Y_m = c_m + d_m \\ c_m = P_{m+1} \cdot 2^{m + 1} \\ d_m = (a_m)^2 \end{cases}$$ 而每次迴圈更新則是: $c_{m - 1} = P_m \cdot 2^m = (P_{m + 1} + a_m) \cdot 2^m = \begin{cases} \frac{c_m}{2} + d_m, if \ a_m = 2^m \\ \frac{c_m}{2}, if \ a_m = 0 \end{cases}$ $d_{m - 1} = \frac{d_m}{4}$ 這個計算會不斷進行下去,直到 $c_{-1} = P_0 = N, d_0 = 0$ 的時候停止。 ### 測驗二 為了避免使用 / 和 % 來計算除以 $10$ 的商和餘,因此我們要找到 $\frac{1}{10}$ 的二進位表示法(除以 $10$ 可以看作乘以 $\frac{1}{10}$、同時得出除以 $10$ 的結果之後就能得到餘數),但 $10$ 不符合 $2^k \pm 1$ 的格式,因此無法以二進位或 digit recurssion 來表示,只能以一個非常接近 $\frac{1}{10}$ 的數字 $\frac{a}{2^n}$ 來表示 $10$ 的倒數。 題目這邊挑選 $2^N = 128, a = 13, \frac{128}{13} \approx 9.84,因此我們可以把 `x / 10` 的計算轉換成 `x * 13 / 128`,再進一步拆解成 `x * 13 / 8 / 16`。 其中 `x * 13 / 8` 可以透過 `(x >> 3) + (x >> 1) + x` 得到,為了將其還原成 `x * 13` 會再進行一次左移,不過這樣會遺失最低三位元的值,因此在操作前先將最低的三為元的值先存起來、操作完成後再加回去,得到 `x * 13` 之後再左移七位(除以 $128$)。 > 這邊其實也可以直接再 `(x >> 3) + (x >> 1) + x` 後直接右移四,計算結果是一樣的。 > 反正丟失的三位最後也是會被右移七吃掉,在這個情況下可以不用考慮儲存起來。 ```c unsigned d0, d1, d2, q; d0 = x & 0b1; d1 = x & 0b11; d3 = x & 0b111; q = ((((x >> 3 ) + (x >> 1) + x) << 3) + d0 + d1 + d2) >> 7; q = ((x >> 3) + (x >> 1) + x) >> 4; ``` 而餘數再使用餘數定理而得: ```c unsigned r = x - (((q << 2) + q) << 1); ``` 完整的函數可以寫成這樣: ```c void calc_div_mod_10 ( uint32_t in, uint32_t *div, uint32_t *mod ) { unsigned d0, d1, d2; d0 = in & 0b1; d1 = in & 0b11; d3 = in & 0b111; *div = ((((in >> 3 ) + (in >> 1) + in) << 3) + d0 + d1 + d2) >> 7; *mod = in - (((*div << 2) + *div) << 1); } ``` 最後再來看一下包裝的結果: ```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)); } ``` 看似跟原先我們計算出來的結果很不相同,好像沒辦法很順利的轉換;那麼我們先來分析一下。 我們檢視前面的計算,發現在函數第一行將 $\frac{3}{4} \cdot in$ 計算出來了,後面再加上一個 $\frac{3}{64} \cdot in$ 也就是 $\frac{51}{64} \cdot in$,而 $\frac{51}{64} \approx 0.797$。 後面再加上 `q >> 8` 也就是 $(\frac{51}{64} \cdot in) \cdot \frac{257}{256} = \frac{13107}{16384} \cdot in \approx 0.7999$,經過額外三次的操作後 $q = \frac{219902325555}{274877906944} \approx 0.8 \cdot in$。 也就是說,到最後要真正計算出 `div` 之前,前面已經將 $in \cdot \frac{13}{16}$ 的部分計算完了,所以最後只要再將最後的 $8$ 除掉就好了。 而計算 `mod` 時,使用 `q & ~0x7` 用意是將最後三個位元的值清除,作用和 `*div << 3` (`*div * 8`)相同,從這邊可以推敲出 `*div << DDDD` 就是要將剩下的值補上。 ### 測驗三 原先的作法有點類似二分搜,從高位開始檢查,如果輸入大於 $2^16$ 會將輸入的數字左移 16 位、並將回傳值加上 16;再來檢查是否大於 $2^8$ 依此類推。 另外我們可以發現 `ilog2(n)` 就是在計算 n 從 msb 數來第一個 set bit 的位置,因此就能使用 `__builtin_clz(n)` 先計算 `n` 的 leading zeros,接著用 31 減掉它也就是結果。 ### 測驗四 EMWA 公式為: $$S_t = \begin{cases} \alpha \cdot Y_t + ( 1 - \alpha ) \cdot S_{t - 1}, & t > 0 \\ Y_0, & t = 0 \end{cases}$$ 從 `ewma_init` 函數中可以發現,一開始初始化的時候會將 `weight`, `factor` 先對數化,所以平移等量的次數可以視為乘除。 而 `ewma_add` 則是對應到 EWMA 的公式,其中比較特別的是每次資料點會先乘上(左移) `factor` 倍才帶入公式中。 `EEEE` 為 `avg->weight`、`FFFF` 為 `avg->factor`。 ### 測驗五 經過觀察發現,測驗三的每個 while 迴圈其實只會執行一次,也就是說我們可以改成 if 的形式實作。 > 舉例來說,如果 $x \ge 2^8$ 這一項會執行兩次,代表 $x \ge 2^8\cdot 2^8 = 2^{16}$,這種情況應該會在前一項 $x \ge 2^{16}$ 時就被篩掉了。 而改為 if 的形式,則代表說我們可以改寫成不會造成分支的寫法,也就是原先題目中給的樣式,並將位移量保留在變數 `r` 中。 ## 第四週 ### 測驗一 總之,要計算出所有組合的 Hamming Distance 和。 ```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; } ``` 觀察一下會發現,所有組合都會被重複算到(舉例來說,`nums[1], nums[2]` 和 `nums[2], nums[1]` 都會被分別計算一次),因此最後要將答案除以二,也就是右移一。 另外我們可以把 `j` 的起始值設成 `i+1`,這樣就能夠避免重複計算的問題了。 #### 避免使用 popcount 的改進 即便我們從根本上減少一半的計算數量,但還是沒有將複雜度下降(一樣為 $O(N^2)$,因此需要尋找其他作法。 我們換個角度思考,如果我們改成計算出 **某位數為一的數字**,同時也能算出 **某位數不為一的數字**(用全部減掉該位數為一的數量),將這兩個數字相乘,就會是該位數的 Hamming Distance,並將所有位數的 Hamming Distance 加總在一起即可。 所以程式碼可以改寫成這樣: ```c int HammingDistantLinear(int* nums, int numsSize) { int res = 0; for (int i = 0, cnt = 0, it = 1; i < 32; i++, cnt = 0, it <<= 1) { for (int j = 0; j < numsSize; j++) if (nums[j] & it) cnt++; res += (numsSize - cnt) * cnt; } return res; } ``` ### 測驗二 tictactoe 的棋局一共會有九種可能的選擇,因此我們可以將棋盤以二進位的形式紀錄(對於該玩家而言,哪一位置有他的記號),並且將整個棋盤一共八種路線(三條橫、三條竪、兩條對角線)分別紀錄下來,這樣只需要 24 位元就能夠很詳盡的紀錄完所有資訊,不需要額外做判斷(是否能連成線等等)。 但實際上實作過後才會發現,這樣還是不便於紀錄,那麼我們將其中一條連線的狀態取出觀察,會發現三個位元會是以 $(111)_2$ 的方式呈現,那麼我們便可以為每一條路線新增一位,這樣將這條路線的值 +1 的時候就會成為 $(1000)_2$ 的狀態,可以更方便的做運算;且即便新增一位,總位元數也還在 `uint32_t` 的範圍內,同時每四個位元為一組可以更輕易的轉換成 hex 編碼。 接著往下看,在 `mod7` 函數中,其中 $(24924925)_{16} \approx \frac{2^{32}}{7}$,從 [Hacker's Delight](https://doc.lagout.org/security/Hackers%20Delight.pdf) 中知道可以利用 $8 \equiv 1 \pmod{7}$ 此數學關係是,可以推導出 $n \pmod{7} = popcount(n)$。 例如書中的其中一段程式碼: ```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); n = (n >> 9) + (n & 0x001FF); n = (n >> 6) + (n & 0x0003F); return table[n]; } ``` 可以發現是利用 $n \pmod{7} = popcount(n)$ 不斷的進行 digit sum,每一步都在重複求 modulo 7。 回來觀察題目中的程式碼,可以發現該程式碼中另外運用了 $n \pmod{7} \equiv \lfloor \frac{8}{7} n \rfloor \pmod{8}$,而我們可以利用 $\lfloor \frac{2^{32}}{7} \rfloor$ 再右移 29 位已取得 $\lfloor \frac{8}{7} n \rfloor \pmod{8}$ 的值,雖然取 floor 會有誤差(例如題目中的程式碼以 `0x24924925` 作為 $\lceil \frac{2^{32}}{7} \rceil$ 的值),不過乘上這個數字後許多高位位元會 overflow 被遺棄掉;同時按照操作將最後的數字右移 29 位元後即是想要的答案。 ### 測驗三 按照題目要求補完程式碼,推測是要將目前的節點刪除,並且如果右子樹並非為空則從右子樹中找出適當的節點;如果右子樹為空則從左子樹中尋找;如果兩者皆為空則視此子樹為空。 而按照 AVL 的性質,該節點的右子樹中的所有節點的值、皆要比該節點的值還大,因此從右子樹中需要找出最小的值。 同樣按照 AVL 的性質,該節點中的左子樹中的所有節點的值、皆要比該節點的值還要小,因此從左子樹中需要找出最大的值。 並且將當前節點以上述節點替代時,需要重新平衡該子樹。 ```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; } if (xt_left(del)) { struct xt_node *most = xt_last(xt_left(del)); if (del == *root) *root = most; xt_replace_left(del, most); xt_update(root, xt_left(most)); return; } ``` ```c /* empty node */ struct xt_node *parent = xt_parent(del); if (xt_left(parent) == del) xt_left(parent) = 0; else xt_right(parent) = 0; xt_update(root, parent); ```