# 2024q1 Homework4 (quiz3+4) contributed by < [gawei1206](https://github.com/gawei1206) > ## 第三周測驗題 ### 測驗 1 以 Digit-by-digit calculation 方式求得 $x$ 的平方根 　　　　$x$ = $N^2$ = $(a_n + a_{n-1} + a_{n-2} + ... + a_0)^2$, $a_m$ = $2^n$ $or$ $0$ $P_m$ = $a_n + a_{n-1} + a_{n-2} + ... + a_{m+1} + a_m$，而 $P_0$ 即為 $x$ 的平方根 $N$ 再來可以觀察到： 　　　　$P_m$ = $P_{m+1}$ + $a_m$ 因為現在的目的是要求出 $P_0$，也就是要從 $a_n$ 累加到 $a_0$，在每回合判斷是否加入 $a_m$ \begin{cases} P_m = P_{m+1} + 2^m, & \text{if $P_m^2 \leq N^2$}\\ P_m = P_{m+1}, & \text{otherwise} \end{cases} 同先前版本實作的效果 ```c if ((result + a) * (result + a) <= N) ``` 由於每輪計算 $N^2 - P_m^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} = 2P_{m+1}a_m+a_m^2$ 也就是紀錄上一輪的 $P_{m+1}$ 來計算 $Y_m$。 再進一步將 $Y_m$ 拆成兩個變數相加，$2P_{m+1}a_m+a_m^2 = P_{m+1}2^{m+1} + (2^m)^2$ $$ c_m = P_{m+1}2^{m+1}, \quad d_m = (2^m)^2, \quad Y_m = \begin{cases} c_m+d_m \quad &if\quad a_m = 2^m \\ 0 \quad &if\quad a_m = 0 \end{cases} $$ 藉由 $c_m, d_m$ 算出 $c_{m-1}, d_{m-1}$ $c_{m-1} = P_{m}2^m = (P_{m+1} + a_m)2^m = P_{m+1}2^m+a_m2^m = \begin{cases} c_m/2+d_m & \text{if }a_m=2^m \\ c_m/2 & \text{if }a_m=0 \end{cases}, \quad$ $d_{m-1}=2^{2m-2}=\dfrac{d_m}{4}$ 最後以上述的推導求出 $P_0$ 即從 $a_n$ 累加到 $a_0$ 的和，以下為參數的對應： z : $c_m$ m : $d_m$ b : $Y_m$ x : $X_{m+1}$ ```c 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; } ``` 其中 `if (x >= b)` 就是以 $X_{m+1}$ 和 $Y_m$ 判斷 $X_{m}$，即 $N^2 - P_m^2$ #### 延伸問題 以 `ffs` 取代 `__builtin_clz`: 先前的方法是以 `(31 - __builtin_clz(x))` 來計算 MSB，`ffs()`的回傳是目前 LSB 的位置，以每次求得的 LSB 位置累加並右移，也可以計算出 MSB。 :::info 但不知道不依賴 GNU extension 會有什麼優點 ::: ```c int tmp = x, msb = 0; while (tmp) { int s = ffs(tmp); msb += s; tmp >>= s; } msb--; for (int m = 1UL << (msb & ~1UL); m; m >>= 2) ``` 提供分支和無分支 (branchless) 的實作: 參考 `LindaTing0106` 同學的作法，將迴圈中的判別式改寫，達到無分支的實作。 ```diff= -if (x >= b) - x -= b, z += m; +int mask = (x >= b); +x -= b * mask; +z += m * mask; ``` ### 測驗 2 採用 bitwise operation 來實作除法，但是除數無法直接用 2 的冪來表示的話結果會是不精確的，因此我們的目標是計算新的除數，使 `tmp / 10` 直到小數點後第一位都是精確的。 題目的除數為 `10`，再來因為`a[i]`, `b[i]` 的範圍是 `0 ~ 9`，`carry` 則是 `0 ~ 1`，，所以 `tmp` 的範圍是 `0 ~ 19`。 ```c int tmp = (b[i] - '0') + (a[i] - '0') + carry; carry = tmp / 10; tmp = tmp - carry * 10; ``` $x = \frac{2^N}{a}$ 為這邊想要求得的新除數，可以用 2 的冪來表示且 $\frac{tmp}{x}$ 仍然會在精確度內，其中 $N$ 為非負整數 $a$ 為正整數。 證明猜想成立後，接下來只需要針對 `19` 來做討論即可。 $a.b \le \frac{n}{x} \le a.b9 \quad\Rightarrow\quad 1.9 \le \frac{19}{x} \le 1.99 \quad\Rightarrow\quad 9.55 \le x \le 10$ 接下來需要找出可以符合條件的 $N$ 和 $a$，題目以 $2^N = 128，a = 13$ 為其中一解，使得 $\frac{128}{13}\approx9.84$ 符合條件。 嘗試透過 bitwise operation 來配對數值，首先要得到 $\frac{13tmp}{8} = \frac{tmp}{8}+\frac{tmp}{2}+tmp$ ```c (tmp >> 3) + (tmp >> 1) + tmp ``` 接下來再左移 3 bit 後就可以得到 $13tmp$，但在右移時會有部分位元被捨棄，因此要另行處理。 ```c ((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2) ``` :::info 但在這部分我有些疑惑，因為上述程式碼所計算出來的結果有時可能不是 $13tmp$，而只是接近 $13tmp$ ::: 雖然最後要再右移 7 bit 結果應該是不影響，但如果真的要將被捨棄的位元加回去，應該如下較為正確： ```c q = ((((x >> 3) + (x >> 1) + x) << 3) + (d0 << 2) + d2) >> 7; r = tmp - (((q << 2) + q) << 1); ``` 再透過 $r = f - g\cdot Q$ 算出餘數。 最後包裝成題目的程式碼如下 ```c 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) - (in >> 2)` 算出 $x = \frac{3}{4} \cdot in$，其中 `(in | 1)` 表示當 `in` 為偶數時會將它加一，再來 `q = (x >> 4) + x` 的意思是 $\frac{1}{16} \cdot \frac{3}{4} \cdot in + \frac{3}{4} \cdot in = \frac{51}{64} \cdot in \approx 0.79 \cdot in$，後續四行重複的操作是為了將結果逼近 $0.8 \cdot in$，最後右移 3 bit 得到 $\frac{1}{10} \cdot in$。 一樣透過 $r = f - g\cdot Q$ 算出餘數，`q & ~0x7` 與 `*div << 3` 是一樣意思的，再加上 `*div << 1` 求得 $10 \cdot div$，即可計算出餘數。 ### 測驗 3 實作 `ilog2` 的目標是找出二進位表示的 msb，在第二種實作方法中的想法是以二元搜尋快速找到 msb，相較於方法一逐位元的判斷來說更有效率。 運用 GNU extension `__builtin_clz` 來改寫，為了避免 `0` 的例外狀況需將輸入與 `1` 做 or。 ```c int ilog32(uint32_t v) { return (31 - __builtin_clz(v | 1)); } ``` ### 測驗 4 在對 `ewma` 結構體初始化時，會確保 `weight` 及 `factor` 的值是 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; } ``` 在後續的程式碼中需要更新 $S_t$ 的值，在等號右邊的 `internal` 代表的是 $S_{t-1}$，而 `val` 代表的是 $Y_t$，由下方的程式碼可以發現，在 `internal` 是 0 的時候程式會執行 `(val << avg->factor)`，所以這邊推測上方的 FFFF 也是做相同的操作，再來觀察 `(((avg->internal << EEEE) - avg->internal) + (val << avg->factor)) >> avg->weight` 這段程式碼與 $\alpha Y_t + (1 - \alpha)\cdot S_{t-1}$ 的關係，可以轉換成 $((S_{t-1} \cdot2^E - S_{t-1})+Y_t) \times \ 2^{-weight}$，可以發現到 $2^{-weight}$ 就是 $\alpha$ ，再將 $2^{-weight}$ 乘進去想要算出 $(1 - \alpha)\cdot S_{t-1}$，所以 EEEE 就是 `avg->weight`。 ```c struct ewma *ewma_add(struct ewma *avg, unsigned long val) { avg->internal = avg->internal ? (((avg->internal << EEEE) - avg->internal) + (val << FFFF)) >> avg->weight : (val << avg->factor); return avg; } ``` ### 測驗 5 與第三題類似，目的是要算出 MSB，但這邊要求出 $\lceil log_2(x) \rceil$，所以當 `x` 為二的冪時需要先將 `x--` ，使得最後結果正確。 ```c ... r = (x > 0xFFFF) << 4; x >>= r; shift = (x > 0xFF) << 3; x >>= shift; r |= shift; ... ``` 後續步驟與第三題找出 MSB 的部分相同，不過是以 branchless 的方式實作，以變數 `r` 紀錄位移量，因為位移量都是二的冪，所以 `r |= shift` 可以看做 `r += shift`， 最後 `return (r | shift | x > 1) + 1;`可以看做兩部分，`r | shift` 是將前一次的 `shift` 加上，而在最後一次位移後 `x` 的值可能會是 0x3 或是 0x2，所以後半部會判斷 `x` 是否大於 1 並將它加上，最後再加上 1 達到取 ceil 的效果。 #### 延伸問題 改進程式碼，使其得以處理 x = 0 的狀況，並仍是 branchless： 我們需要判斷 `x` 是否大於 0 再將其減一，在之前的作業中有看過 likely 還有 unlikely，其中這兩個巨集的實作中有用到 `!!(x)` 的用法，目的是透過兩次 NOT 來確保值一定是 0 或 1，所以可將程式中 `x--` 的部分修改如下： ```diff + x -= !!(x); - x--; ```