# 2024q1 Homework4 (quiz3+4) contributed by < [`padaray`](https://github.com/padaray) > ## quiz3 ### 測驗 `1` #### 計算平方根 **版本一：** 對輸入的參數 `N` 以 2 為底數取 log，意義是當我們將 `N` 取 log2 時，可以得到 most significant bit (msb)。 對 1 左移 msb 個位元後存入變數 `a` ，意義是找出一個必定大於 `N` 的平方根的數，利用此數來逼近找到 `N` 的平方根。 最後使用 while 迴圈，計算 `a + result` 的平方是否大於 `N`，若小於 `N` 則將此數存入 `result` 變數，對 `a` 右移 1 bit，繼續計算，直到 `a` = 0。 </br> **版本二：** 版本二的差別是不使用 `math.h` 的 `log2` 函式，直接對參數 `N` 進行右移，右移至參數為 0 時，則代表找到 most significant bit。此方法將程式的可移植性增強，並且不影響原本的函式原型 (prototype) 和精度 (precision) ```diff - int msb = (int) log2(N); + int msb = 0; + int n = N; + while (n > 1) { + n >>= 1; + msb++; + } ``` </br> **版本三：** 完整程式碼(無遮蔽處) ```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; } ``` `__builtin_clz(x)` ：是 gcc 內建處理二進位的其中一個函式巨集，此函式會回傳 counting leading zero，幫助我們找到 most significant bit。 在 for 迴圈中左移 `31 - __builtin_clz(x)` bit 數得到 msb，`& ~1UL` 可以確保最後一個 bit 為 0，理由是變數 `m` 必須是偶數對應到 $d_n = (2^n)^2$。而 m 右移兩個 bit 是對應到 $d_{m-1} = d_m/4$ </br> ### 測驗 `2` #### 針對正整數在相鄰敘述進行 mod 10 和 div 10 操作，如何減少運算成本？ 除法原理：$f=g \times Q + r$ * $f$: 被除數 * $g$: 除數 * $Q$: 商 * $r$: 餘數 程式碼如下： ```c carry = tmp / 10; // carry 代表商(Q) tmp = tmp - carry * 10; // tmp 代表餘數(r) ``` 如果我們想要使用 bitwise operation 來實作除法，如果被除數包含了除了 2 以外的因數，就會造成結果不精確，所以使用以下方式。 假設除數為 10，我們必須找到 $\dfrac{1}{10} \approx \dfrac{a}{2^N}$，最後發現是 $\approx \dfrac{13}{128}$ 程式碼如下： ```c unsigned d2, d1, d0, q, r; d0 = q & 0b1; d1 = q & 0b11; d2 = q & 0b111; q = ((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2) >> 7; r = tmp - (((q << 2) + q) << 1); ``` `(((tmp >> 3) + (tmp >> 1) + tmp) << 3)` 使用 bitwise operation 取得一個分子有 13 的分數，因此用 $\dfrac{tmp}{8} + \dfrac{tmp}{2} + tmp = \dfrac{13tmp}{8}$，之後再左移三位消掉分母 8。 `d0`、`d1`、`d2` 是為了將 `tmp` 右移時損失的 bit 數加回，加回後再右移 7 位也就是除以 128 得到商數。 `(((q << 2) + q) << 1)` 相當於 $(q \times 4 + q) \times 2$，也就是最一開始程式碼的 `carry * 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)); } ``` </br> ### 測驗 `3` #### 考慮 ilog2 可計算以 2 為底的對數，且其輸入和輸出皆為整數。 **版本一：** 初始化變數 `log` 為負一，直接對參數 `i` 進行右移，每次右移將 `log++`，右移至參數 `i` 為 0 時 </br> **版本二：** 完整程式碼(無遮蔽處) ```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; } ``` 為了避免數字過大每次位移一個 bit 太慢，例如：65536 需要位移 16 個 bit，因此加入條件判斷，若大於特定 2 的冪次方，直接移動該次方的位元數。 </br> **版本三：** ```c int ilog32(uint32_t v) { return (31 - __builtin_clz(v)); } ``` </br> ### 測驗 `4` #### 解釋 EWMA 的實作 Exponentially Weighted Moving Average (EWMA) 是一種統計方法，透過對每筆資料給予不同的權重計算平均，通常越新的資料權重就越高，越舊的資料權重越小，好處是可以適應新的資料變化，又保留一定的歷史資料影響。 數學定義： $S_t = \begin{cases} Y_0, & t=0 \\ \alpha Y_t + (1-\alpha)\cdot S_{t-1}, & t > 0 \end{cases}$ $\alpha$ : 舊資料加權降低的程度，越高代表降低越快。介於 0 ~ 1 之間 $Y_t$ : 時間 t 時的資料 $S_t$ : 時間 t 時算出的 EWMA 以下為程式碼實作： ```c * ewma_add() - Exponentially weighted moving average (EWMA) * @avg: Average structure * @val: Current value * * Add a sample to the average. */ 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; } ``` `internal` 代表目前計算的 ewma (數學式中的 $S_t$) ， `(avg->internal << avg->weight)` 代表數學式中的 $\alpha = \dfrac{1}{2^{weight}}$ (**但不知道原因為何**)，因此 $1- \alpha = \dfrac{2^{weight} - 1}{2^{weight}}$。 當 `internal` 為 0 時，`internal` 代入 `(val << avg->factor)` 寫成公式為 $2^{factor} \cdot val$。 當 `internal` 不為 0 時，`internal` 代入 `(((avg->internal << avg->weight) - avg->internal) + (val << avg->factor)) >> avg->weight` 寫成公式為 $\dfrac{2^{weight}-1}{2^{weight}} \cdot internal + \dfrac {2^{factor}}{2^{weight}} \cdot val$，轉換一下可以得到 $(1-\alpha) \cdot internal + {2^{factor}} \cdot \alpha \cdot val$ </br> ### 測驗 `5` #### 延伸測驗 `3` 使用下方程式碼計算 $\lceil \log_2 x\rceil$ ```c= int ceil_ilog2(uint32_t x) { uint32_t r, shift; x--; r = (x > 0xFFFF) << 4; x >>= r; shift = (x > 0xFF) << 3; x >>= shift; r |= shift; shift = (x > 0xF) << 2; x >>= shift; r |= shift; shift = (x > 0x3) << 1; x >>= shift; return (r | shift | x > 1) + 1; } ``` 第五、六行判斷參數 `x` 是否大於 0xFFFF，意義是確定 `x` 的 msb 是否位於高位的 16 bit，若為真，將 `x` 右移 16 bit。 第七、八、九行判斷參數 `x` 在這剩下的 16 bit 中，msb 是否位於高位的 8 bit，若為真，將 `x` 右移 8 bit，將右移的位數紀錄在 `r` 中。 第十行到第十四行依此類推。 第四、十五行將 `x` 先 -1 後 + 1 回傳是為了處理 2 的冪次方，若不處理會讓 2 的冪次方多計算一次，例：$2^7$ 經過無條件進位法，會變成 8，因此要先減 1 </br></br> ## quiz4 ### 測驗 `1` #### [popcount(population count)](https://en.wikichip.org/wiki/population_count)