# 2024q1 Homework4 (quiz3+4)
contributed by < [Petakuo](https://github.com/Petakuo) >
## 第三週測驗題
## [測驗1](https://hackmd.io/@sysprog/linux2024-quiz3#%E6%B8%AC%E9%A9%97-1)
### 版本一
首先，利用 `log2()` 函式找到最高位的位元，接著由該位元開始測試，逐漸逼近答案，而逼近的方法為利用 `if` 函示判斷 `(result + a)` 的平方是否大於被要求開平方的值 `N` ，如果有，則繼續往低位進行逼近，如果沒有，則將原本的答案加上 `a` 進行更新，並同樣往低位進行逼近，如此做到 `a` 為 0 時即可得到 `N` 的平方根。
### 版本二
```
int msb = 0;
int n = N;
while (n > 1) {
n >>= 1;
msb++;
}
```
與版本一大致相同，唯改變了取最高位元的做法，將 `log2()` 用以上程式碼呈現，首先判斷 input 的 `n` 是否大於 1 ，有則將其右移一個位元，並將所求 +1 ，如此進行直到 `n` 只剩下 1 即代表該位元為最高位元，而所求也隨著每一輪迴圈的 +1 而得出。
### 版本三
利用 [Digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) 的方法計算，先將要求平方根的 $N^2$ 拆解成 ：
$N^2 = (a_n+a_{n-1}+a_{n-2}+...+a_1+a_0)^2$
並且先觀察較簡單的例子 ：
$(a_n+a_{n-1})^2 = a_n^2+2a_na_{n-1}+a_{n-1}^2 = a_n^2+[2a_n+a_{n-1}]a_{n-1}$
$(a_n+a_{n-1}+a_{n-2})^2 = a_n^2+[2a_n+a_{n-1}]a_{n-1}+[2(a_n+a_{n-1})+a_{n-2}]a_{n-2}$
依此模式進行，我們最終就可以得到 ：
$N^2 = a_n^2+[2a_n+a_{n-1}]a_{n-1}+[2(a_n+a_{n-1})+a_{n-2}]a_{n-2}+...+[2(\sum_{i=1}^na_i)+a_0]a_0$
此時再令 $P_m = a_n+a_{n-1}+...+a_m$ ，則上述式子可以進一步被改寫為 ：
$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_0 = a_n+a_{n-1}+...+a_0$ 即為所求 $N$
為求 $P_0$ ，我們可以由 $P_m = P_{m-1}+a_m$ 進行迭代，每一輪皆去判斷 $P_m^2 \leq N^2$ 是否成立，若成立則將該位元設為 1 ，反之則設為 0 。
但此時由於去計算每一輪 $P_m^2$ 的成本太高，因此我們透過利用上一輪的資訊來做計算，首先定義
$X_m = N^2-P_m^2 = X_{m+1}-Y_m$
並且
$Y_m = P_m^2-P_{m+1}^2 = 2P_{m+1}a_m+a_m^2$
此時再將 $Y_m$ 改寫為 $c_m+d_m$ ，其中 $c_m = P_{m+1}a_{m+1} , d_m = a_m^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+d_m &\text{if }a_m = 2^m\\
c_m/2 &\text{if }a_m = 0\\
\end{cases}$
$d_{m-1} = d_m/4$
到這裡可以發現，將 m = 0 代入，可得 $c_{-1} = P_02^0 = P_0 = a_n+a_{n-1}+...+a_0 = N$ ，即為所求。
```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 >>= AAAA) {
int b = z + m;
z >>= BBBB;
if (x >= b)
x -= b, z += m;
}
return z;
}
```
在此程式碼中，利用 `__builtin_clz(x)` 找出 msb ，且 b 為上述推導過程中的 $Y_m$ ，因此 z 等同於 $c_m$ 而 m 等同於 $d_m$ 。由於每一輪 $d_m$ 都會除以 4 ，因此 `AAAA` 為 2 ，同樣道理， $c_m$ 每一輪會除以 2 ，所以 `BBBB` 為 1 。
延伸問題 ： 利用 `ffs` 取代 `__builtin_clz` ，使程式不依賴 GNU extension
```c
static inline unsigned long __ffs(unsigned long word)
{
int num = 0;
if ((word & 0xffff) == 0) {
num += 16;
word >>= 16;
}
if ((word & 0xff) == 0) {
num += 8;
word >>= 8;
}
if ((word & 0xf) == 0) {
num += 4;
word >>= 4;
}
if ((word & 0x3) == 0) {
num += 2;
word >>= 2;
}
if ((word & 0x1) == 0)
num += 1;
return num;
}
int i_sqrt(int x)
{
if (x <= 1) /* Assume x is always positive */
return x;
int z = 0;
for (int m = 1UL << ((31 - __ffs(x)) & ~1UL); m; m >>= 2) {
int b = z + m;
z >>= 1;
if (x >= b)
x -= b, z += m;
}
return z;
}
```
## [測驗2](https://hackmd.io/@sysprog/linux2024-quiz3#%E6%B8%AC%E9%A9%97-2)
```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);
```
此段程式碼嘗試利用 bitwise 操作來取代正整數的除法和 mod 運算，這裡以 10 為例進行解說。
從數學的角度切入，除以 10 等同於乘以 $\cfrac{1}{10}$ ，因此我們必須找到在二進制中 $\cfrac{1}{10}$ 的表示法，而因為不是所有數都能夠利用有限個 bits 來表達，所以只能取近似值來做運算，在此例中使用了 $\cfrac{128}{13}\approx 9.84$ 來計算，接著，就是要想辦法透過 bit operation 來得出原數乘上 $\cfrac{13}{128}$ 的結果。
首先， $\cfrac{13}{128}$ 可以被拆解成 $\cfrac{13}{8}$ 再除以 16 ，而 $\cfrac{13}{8}$ 又可以被拆解成 $1+\cfrac{1}{2}+\cfrac{1}{8}$ ，如此一來，原數乘上 $\cfrac{13}{8}$ 就可以用 `q = (tmp >> 3) + (tmp >> 1) + tmp` 計算出來，最後再透過 `q = q << 4` 來除以 16 即可。但若考慮到 mod 的運算，則在利用位元右移計算乘以 2 的冪時所被捨棄的位元就會造成計算上的錯誤，因此該程式碼利用了以下三行將原數最右邊的三個位元先保留下來。
```c
d0 = q & 0b1;
d1 = q & 0b11;
d2 = q & 0b111;
```
在計算完 $\cfrac{13tmp}{8}$ 後，再將其乘上 8 (`q = q << 3`)，目的是在於要加回被捨棄掉的 bits ，而在加完之後，再直接除以 128 (`q = q >> 7`) 即可。
至於 mod 的部分，可以直接利用除法原理得到，因為除法原理告訴我們 $tmp = q\times 10 + r$ ，所以餘數可由 $r = tmp - q\times 10$ 得出，因此在程式碼的最後寫了 `((q << 2) + q) << 1` ，即 $(q\times 4 + q)\times 2 = q\times 10$。
將此概念包裝後，可以下方程式碼重新呈現 ：
```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 >> CCCC);
*mod = in - ((q & ~0x7) + (*div << DDDD));
}
```
這段程式碼做的事情為將 `in` 這個數以 bitwise 操作同時除以和 mod 10 ，具體的想法為除以 10 等於乘以 0.1 ，因此可先乘上 0.8 再除以 8 。
第一行 `(in | 1) - (in >> 2)` 強制將 `in` 轉為奇數（不知道原因，且不管有沒有那個 `+1` 經過測試後的皆是正確的）並減去 $\cfrac{in}{4}$ 得到 $\cfrac{3in}{4}$ ，而第二行 `(x >> 4) + x` 則是在計算 $\cfrac{3in}{64} + \cfrac{3in}{4}$ ，此計算結果為 $\cfrac{51in}{64} \approx 0.79in$ 即代表前兩行做的事情為令一變數 `q` 等於 `0.8x` 。
而後續的 `q = (q >> 8) + x` 則是在對 0.8 逼近，做越多次就會越接近，最後需再將結果除以 8 ，因此 `CCCC` 為 3 。
餘數的部份，靠除法原理可知 `*mod = in - 10*div` ，因此 `(q & ~0x7) + (*div << DDDD)` 就會等於 `10*div` ，而 `(q & ~0x7)` 做的事情為將 `q` 最低位的 3 個位元清零，此做法等同於 `(q >> 3) << 3` ，即 `8*div` ，因此 `*div << DDDD` 就必須等於 `2*div` ， 最終可知 `DDDD` 為 1 。
延伸問題 ： 練習撰寫不依賴任何除法指令的 % 9 (modulo 9) 和 % 5 (modulo 5) 程式碼，並證明在 uint32_t 的有效範圍可達到等效。
```c
void divmod_9(uint32_t in, uint32_t *div, uint32_t *mod)
{
uint32_t x = in - (in >> 3); /*in = (in*8/9)/8 8/9 = 7/8 + 7/512 */
uint32_t q = (x >> 6) + x;
x = q;
q = (q >> 12) + x;
q = (q >> 12) + x;
q = (q >> 12) + x;
q = (q >> 12) + x;
*div = (q >> 3);
*mod = in - (q & ~0x7) - *div;
if(*mod == 10 || *mod == 9)
{
*div++;
*mod -= 9;
}
}
```
測試程式碼如下，首先撰寫出 `randnum` 來造一個隨機的 32 位元無號數，接下來利用實際的除法與 mod 運算來確認 `divmod_9` 是否正確 ：
```c
uint32_t randnum()
{
uint32_t ran = 0;
ran |= (uint32_t)(rand() & 0xFFFF) << 16;
ran |= (uint32_t)rand() & 0xFFFF;
return ran;
}
int main()
{
uint32_t div, mod;
for(int i = 0; i < 10; i++){
uint32_t ran = randnum();
uint32_t a = ran/9;
uint32_t b = ran%9;
divmod_9(ran, &div, &mod);
printf("a = %d\ndiv = %u\nb = %d\nmod = %u\n", a, div, b, mod);
}
return 0;
}
```
列印出的 10 組數據如下 ：
```
a = 129376363
div = 129376363
b = 3
mod = 3
a = 284115184
div = 284115184
b = 3
mod = 3
a = 410702193
div = 410702193
b = 6
mod = 6
a = 276433377
div = 276433377
b = 3
mod = 3
a = 58090291
div = 58090291
b = 2
mod = 2
a = 165404435
div = 165404435
b = 0
mod = 0
a = 122931853
div = 122931853
b = 6
mod = 6
a = 316698205
div = 316698205
b = 1
mod = 1
a = 905749
div = 905749
b = 5
mod = 5
a = 15525887
div = 15525887
b = 1
mod = 1
```
看起來是正確的，接著再改寫 main 的程式碼，把測試資料放大到 1000000 組 ：
```c
int main()
{
uint32_t div, mod;
int correct = 0;
for(int i = 0; i < 1000000; i++){
uint32_t ran = randnum();
uint32_t a = ran/9;
uint32_t b = ran%9;
divmod_9(ran, &div, &mod);
if(a == div && b == mod)
{
correct++;
}
}
printf("%d/1000000\n", correct);
return 0;
}
```
輸出的結果為 ：
```
1000000/1000000
```
因此可驗證 `divmod_9` 為可行的程式碼。
而對於 mod5 也可以用同樣的方式驗證，程式碼如下 ：
```c
uint32_t randnum()
{
uint32_t ran = 0;
ran |= (uint32_t)(rand() & 0xFFFF) << 16;
ran |= (uint32_t)rand() & 0xFFFF;
return ran;
}
void divmod_5(uint32_t in, uint32_t *div, uint32_t *mod)
{
uint32_t x = in - (in >> 2); /*in = (in*4/5)/4 4/5 = 3/4 + 3/64 */
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 >> 2);
*mod = in - (q & ~0x3) - *div;
if(*mod == 6 || *mod == 5)
{
(*div)++;
*mod -= 5;
}
}
int main()
{
uint32_t div, mod;
int correct = 0;
for(int i = 0; i < 1000000; i++){
uint32_t ran = randnum();
uint32_t a = ran/5;
uint32_t b = ran%5;
divmod_5(ran, &div, &mod);
if(a == div && b == mod)
{
correct++;
}
}
printf("%d/1000000\n", correct);
return 0;
}
```
輸出的結果為 ：
```
1000000/1000000
```
## [測驗3](https://hackmd.io/@sysprog/linux2024-quiz3#%E6%B8%AC%E9%A9%97-3)
### 版本一
首先將 `log` 設置為 -1 ，這是考慮到 input 為 0 的情況，因為 $2^0$ 為 1 ，接著進入迴圈，若 input 不為 0 ，則將 `log` + 1 ，並且同時將 input 右移一個位元，也就是除以 2 ，在離開迴圈時，答案就會是最高位元的 1 位置，也等於對 input 取 log 的整數值。
### 版本二
```c
static size_t ilog2(size_t i)
{
size_t result = 0;
while (i >= AAAA) {
result += 16;
i >>= 16;
}
while (i >= BBBB) {
result += 8;
i >>= 8;
}
while (i >= CCCC) {
result += 4;
i >>= 4;
}
while (i >= 2) {
result += 1;
i >>= 1;
}
return result;
}
```
與版本一的概念相同，都是利用迴圈判斷是否將 input 右移，只是在這個版本利用多個判斷式子來加速進行，而判斷的依據為是否大於 $2^{16}$ $2^8$ $2^4$ 以及 $2^2$ ，因此 `AAAA` , `BBBB` , `CCCC` 分別為 65536 , 256 以及 16 。
#### 版本三
```c
int ilog32(uint32_t v)
{
return (31 - __builtin_clz(DDDD));
}
```
版本三利用了 GNU extension 中的函式 `__builtin_clz` 來改寫，該函式的功能為返回從最高位數來為 0 的個數，因此若利用總位元數減去該函數之結果就可以得到最高位為 1 的位置，即為所求，所以直覺上可以直接把 `v` 放入該函式中，但還必需注意到該函式在 input 為 0 時未定義，因此需將 `v` 改為 `v|1` ，已確保 `v` 不為 0 ，藉此避免未定義的情況發生，因此 `DDDD` 為 `v|1` 。
## [測驗4](https://hackmd.io/@sysprog/linux2024-quiz3#%E6%B8%AC%E9%A9%97-4)
```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;
}
```
此程式碼用來計算 [EWMA](https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average) ，具體的數學定義為 ：
$S_t = \begin{cases}
Y_0 & t = 0\\
\alpha Y_t + (1 - \alpha)S_{t-1} & t>0
\end{cases}$
$S_t$ 為程式碼中的 `internal` ， $Y_t$ 為 `val << avg->factor` ，會隨著 `factor` 值而改變，而 $\alpha$ 為 $\frac{1}{2^{avg->weight}}$ 。
首先，該程式碼應用了三元運算來對 `t` 是否為 0 的情況分類， 並且先將整個式子除以 $\alpha$ 做運算，因此在程式碼中若 `t` 不為 0 ，則需做 `((avg->internal << EEEE) - avg->internal) + (val << FFFF)` ，即表示 $\frac{1}{\alpha}S_{t-1} - S_{t-1} + Y_t$ ，因此 `EEEE` 為 `avg->weight` ， `FFFF` 為 `avg->factor` ，算完後還需乘回原本的 $\alpha$ ，因此最後會有 `>> avg->weight` 。
## [測驗5](https://hackmd.io/@sysprog/linux2024-quiz3#%E6%B8%AC%E9%A9%97-5)
```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 > GGGG) + 1;
}
```
根據[測驗3](https://hackmd.io/e9bG8Gx9Qia7Un611VwFaw?view#%E6%B8%AC%E9%A9%973)我們可以知道，若要求某數取 log2 的值只要找到該數最高位元的 1 即可，此程式碼也在做相同的事情，利用 `r` 和 `shift` 記錄一共位移了幾個 bits ，最後再將他們相加得出結果。
首先，`(x > 0xFFFF) << 4` 是在判斷 `x` 的左半邊 16 個 bits 中是否有 1 存在，若有則將 `r` 設為 16 並且將 `x` 右半邊 16 個 bits 捨棄，若沒有則將 `r` 設為 0 ， `x` 保持不變。後續就是不斷的對剩餘的位元數做相同的事情，也就是將該數切半並判斷，最後利用 or 運算將 `r` 和 `shift` 以及 `(x > 0x1)` 相加，因此 `GGGG` 為 1 。
在程式碼的最後寫了 `+1` 是對應到了程式碼一開始的 `x--` ，這裡考量了 input 本身為 2 的羃的情況，要注意到當沒有 `x--` 且 input 為 2 的羃時，最後的結果會多 1 ，因此最後利用 `+1` 來修正 ： 若 input 為 2 的羃，則加減 1 後答案不變，若不為 2 的羃，則加 1 可達到取 ceiling 的效果。
延伸問題 ： 改進程式碼，使其得以處理 x = 0 的狀況，並仍是 branchless 。
```c
int ceil_ilog2(uint32_t x)
{
uint32_t r, shift;
x -= !!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 > GGGG) + 1;
}
```
將 `x--` 改寫為 `x -= !!x` ，利用兩次的邏輯運算 `!` 來判斷 `x` 的值是否為 0 ，如此一來若 `x` 為 0 則 - 0 ，若 `x` 為 1 則 - 1 。
## 第四週測驗題
## [測驗1](https://hackmd.io/@sysprog/linux2024-quiz4#%E6%B8%AC%E9%A9%97-1)
```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;
}
```
這段程式碼用來計算一個 32 bits 無號整數中 1 的個數，也就是在實作 popcount() 函式。
首先，假設一 32 bits 無號整數 $x = b_{31}b_{30}\cdots b_{2}b_{1}b_{0}$ ，則該數中 1 的個數可以將每個位元相加而得，也就是 $\sum_{n=0}^{31} b_n = \sum_{n=0}^{31} (2^n - \sum_{i=0}^{n-1} 2^i) b_n$ 。
而另一種做法為利用數學式去推導，可以得知 $x - \lfloor \cfrac{x}{2} \rfloor - \lfloor \cfrac{x}{4} \rfloor - \cdots - \lfloor \cfrac{x}{2^{31}} \rfloor$ 也是一種算法，該程式碼就是利用這個觀念求得。
直觀的做法為不斷的將 `v` 除以二，並且用原本的 `v` 去扣除它，就如上方程式碼中的 `n = (v >> 1)` 和 `v -= n` ，不過在 32 bits 的數值上運作就必須重複做 32 次循環，因此該程式碼將其分為 8 組，也就是將每 4 個位元 (nibble) 視為一個單位去做計算，如此本該重複 32 次迭代的程式碼就能簡寫為上方的形式。
首先從公式的部份討論每組 nibble 該如何計算，因為只有 4 個 bits ，所以可以將公式改寫為$\sum_{n=0}^{3}b_n = \sum_{n=0}^{3}(2^n-\sum_{i=0}^{n-1}2^i)b_n$
展開後可得
$(2^3-2^2-2^1-2^0)b_3+(2^2-2^1-2^0)b_2+(2^1-2^0)b_1+2^0b_0$
而若是以數學是推導，可以得知在 4 個 bits 的情況下可由 $x - \lfloor \cfrac{x}{2} \rfloor - \lfloor \cfrac{x}{4} \rfloor- \lfloor \cfrac{x}{8} \rfloor$ 得出答案，且該式等同於 $(2^3b_3+2^2b_2+2^1b_1+2^0b_0)-(2^2b_3+2^1b_2+2^0b_1)-(2^1b_3+2^0b_2)-2^0b_3$ 並經過改寫後可得 $(2^3-2^2-2^1-2^0)b_3+(2^2-2^1-2^0)b_2+(2^1-2^0)b_1+2^0b_0$ ，與上述推導出的公式相同，驗證了其正確性。
此程式碼便是以這個概念進行運算，首先，`v >> 1` 是在進行 $\lfloor \cfrac{x}{2} \rfloor$ 的運算，而每多做一次就會再除以 2 並取其底，由於我們讓每 4 個 bits 為一組，因此需要做 3 次至 $\lfloor \cfrac{x}{2^3} \rfloor$ ，而因為每次在除以 2 時皆會將 `v` 右移一位，如此做法會造成前三組 nibble 的最高位元在移動後都會等於下一組的最低位元，計算出的結果就會錯誤。
直接列出來解釋如下 :
```
b31 b30 b29 b28 ... b3 b2 b1 b0 // v
0 b31 b30 b29 ... b4 b3 b2 b1 // v >> 1
```
右移一位後， b4 會成為第一組 nibble 的最高位元，但我們在運算時不需要 b4 ，因此需將其設定為 0 ，這裡的做法是和 `0x77777777` 做 `&` 運算，如此一來便可以同時將八組 nibble 的最高位歸零，不過這裡使用 `0xF7777777` 也可以達到相同效果，因為右移後最高位本身就會被 0 補上。
```
b31 b30 b29 b28 ... b3 b2 b1 b0 // v
0 b31 b30 b29 ... 0 b3 b2 b1 // (v >> 1) & 0x77777777
```
重複做四次後，可以得到一組新的 32 bits 表示，令為 $B_7B_6\cdots B_1B_0$ 其中每一 $B_n$ 皆表示一組 nibble ，接著就是要將八組 nibble 所求出來 1 的個數相加，因此我們想要得到的數值為 $B_7+B_6+B_5+B_4+B_3+B_2+B_1+B_0$ ，實際做法如下 :
```
B7 B6 B5 B4 B3 B2 B1 B0 // v
0 B7 B6 B5 B4 B3 B2 B1 // v >> 4
(B7+0) (B6+B7) (B5+B6) (B4+B5) (B3+B4) (B2+B3) (B1+B2) (B0+B1) // v + (v >> 4)
```
先將該數右移一個 nibble 的長度，接著與原數相加，從得到的數觀察可以發現我們想要的數值為偶數組的 nibble 相加，因此再與 `0x0F0F0F0F` 做 `&` 運算，可得 :
```
// & 0x0F0F0F0F
0 (B7+B6) 0 (B5+B4) 0 (B3+B2) 0 (B1+B0)
```
最後，為了把這些數值相加，我們讓 `v` 去和 `0x01010101` 相乘，以直式展開結果如下 :
```
0 (B7+B6) 0 (B5+B4) 0 (B3+B2) 0 (B1+B0)
x 0 1 0 1 0 1 0 1
------------------------------------------------------------------------------
0 (B7+B6) 0 (B5+B4) 0 (B3+B2) 0 (B1+B0)
0 (B7+B6) 0 (B5+B4) 0 (B3+B2) 0 (B1+B0)
0 (B7+B6) 0 (B5+B4) 0 (B3+B2) 0 (B1+B0)
0 (B7+B6) 0 (B5+B4) 0 (B3+B2) 0 (B1+B0)
------------------------------------------------------------------------------
↑________(B7+B6)+(B5+B4)+(B3+B2)+(B1+B0)
```
到此可以看到我們需要的結果在相乘後數值的第七個 nibble 的位置，也就是我們可以將前面 24 bits 捨棄掉，因此在程式碼的最後用了 `v >> 24` 來進行此一操作，最終就能夠得到 $(2^3-2^2-2^1-2^0)b_{31}+(2^2-2^1-2^0)b_{30}+(2^1-2^0)b_{29}+2^0b_{28}+\cdots+(2^3-2^2-2^1-2^0)b_3+(2^2-2^1-2^0)b_2+(2^1-2^0)b_1+2^0b_0$ ，即為所求。
## [測驗2](https://hackmd.io/@sysprog/linux2024-quiz4#%E6%B8%AC%E9%A9%97-2)
本測驗告訴我們如何不用除法就能得出某兩數相除之後的餘數，以下用 mod 3 當作例子進行操作和解釋 :
由數學上同餘 (modulo) 的特性可以知道，當 $a\equiv b(mod\ \ m)$ 且 $c\equiv d(mod\ \ m)$ ，則 $\ \ a+c\equiv b+d(mod\ \ m)$ 且 $ac\equiv bd(mod\ \ m)$ ，因此若我們有 $1\equiv 1(mod\ \ 3)$ 且 $2\equiv -1(mod\ \ 3)$ ，則可以推出 $2^k\equiv\begin{cases}
1(mod\ \ 3) &if\ \ k\ \ is\ \ even\\
-1(mod\ \ 3) &if\ \ k\ \ is\ \ odd\\
\end{cases}$ 。
假設一 32 bits 無號整數 $n = b_{31}b_{30}\cdots b_{2}b_{1}b_{0}$ ，則 n 的值就會等於 $\sum_{i=0}^{31}b_i2^i$ ，由上述推論可得 $n = \sum_{i=0}^{31}b_i(-1)^i\ \ (mod\ \ 3)$ 。到這裡我們可以發現對於以二進位表示的數來講， mod 3 其實就是將偶數位元相加再減去奇數位元，因此可以利用 population count 類型的函式來進行操作，其中一種做法為 `n = popcount(n & 0x55555555) - popcount(n & 0xAAAAAAAA)` ，因為 5 為 0101 ， A 為 1010 ，所以如此運算便能將奇偶穿插著相減，達到預期的效果。
而上述做法其實還有別的寫法，那就是利用以下定理將其簡化 :
$popcount(x\& \overline m)-popcount(x\&m) = popcount(x \oplus m)-popcount(m)$
因此我們就能將其改寫為 `n = popcount(n ^ 0xAAAAAAAA) - 16` ，但這裡必須注意到 `popcount(n ^ 0xAAAAAAAA)` 這個值的範圍為 0 < n < 32 ， -16 這個動作可能會導致算出來的值為負數，因此需透過加上 3 的倍數避免這樣的情況發生，而該數只要大於 16 即可，以下利用 39 做說明。
```c
int mod3(unsigned n)
{
n = popcount(n ^ 0xAAAAAAAA) + 23;
n = popcount(n ^ 0x2A) - 3;
return n + ((n >> 31) & 3);
}
```
程式碼第一行是由上述推倒而來的，也就是 `n = popcount(n ^ 0xAAAAAAAA) - 16 + 39` ，但此時 n 的範圍就會變成 23 < n < 55 ，而我們希望最終的結果會落在 0 < n < 2 (因為 mod 3 )，因此可以再次利用上述定理縮小其範圍。
在經過第一行的運算後，其範圍已經變為 23 < n < 55 ，只有最低位的 6 個 bits 可能會有 1 存在，因此根據定理，第二行即為 `popcount(n ^ 0x2A) - 3` (因為 0x2A 為 101010 )，此時 n 的範圍為 -3 < n < 2 ，我們只需再處理負數的部分即可。
處理的方法為得到負數就將其 + 3 ，首先，透過右移 31 位來判斷正負，若為 1 則代表是負數，需要 + 3 ，於是和 3 做 `&` 運算，此時會發現要用全 1 去做 `&` 才會是 + 3 的結果，也就是右移需要是算術位移而非邏輯位移，但上方的 input 卻是 `unsigned n` ，因此該程式碼有個小 bug 存在，可以透過將 input 宣告成 `int` 來解決，或者是修改 return 的部分，修改過後的程式碼如下 :
```c
int mod3(unsigned n)
{
n = __builtin_popcount(n ^ 0xAAAAAAAA) + 23;
n = __builtin_popcount(n ^ 0x2A) - 3;
return n + (((n >> 31) | (n >> 30)) & 3);
}
```
在最後一行的地方，我將 `(n >> 31)` 與 `(n >> 30)` 做 or 運算，由於 3 的二進位表示法為 011 ，因此只須考慮最後 2 個 bits 即可，並且利用 or 來達到算術位移的效果，也就是將負數的高位元填滿 1 。
而事實上若改寫為這樣也不需要再 `&3` 了，因為若是負數則位移完會得到 `011` ，若是正數則為 `0` ，因此最終的程式碼為 ：
```c
int mod3(unsigned n)
{
n = __builtin_popcount(n ^ 0xAAAAAAAA) + 23;
n = __builtin_popcount(n ^ 0x2A) - 3;
return n + ((n >> 31) | (n >> 30));
}
```
還有另一種做法為利用 lookup table :
```c
int mod3(unsigned n)
{
static char table[33] = {2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1 };
n = popcount(n ^ 0xAAAAAAAA);
return table[n];
}
```
首先要建立一個表格，該表格的內容為 mod 3 後可能的值，也就是 {012} ，而大小則取決於 `n` 的範圍，這裡是 0 ~ 32 ，因此表格大小為 33 。至於表格內值的排序則是由 `n` 所決定的，在這裡因為 `n = popcount(n ^ 0xAAAAAAAA)` ，可以想成 `n = popcount(n & 0x55555555) - popcount(n & 0xAAAAAAAA) + 16` ，因此對於 table[0] 來講，要找到當程式碼內的 `n = 0` 時，實際上 mod 3 為多少，而這裡只要找出前 3 項，然後依序填完整個 table 即可 :
`n = 0` : -16 mod 3 為 2 ，所以 table[0] = 2 ；
`n = 1` : -15 mod 3 為 0 ，所以 table[1] = 0 ；
`n = 2` : -14 mod 3 為 1 ，所以 table[2] = 1 ；
如此填完便可以得到正確的值。
將上述的概念應用於 [tictactoe.c](https://gist.github.com/jserv/9088f6f72a35a1fcaaf1c932399a5624) 程式碼，該程式碼用來模擬 100 萬次井字遊戲的對奕，並統計出第一次選擇 9 個不同位置的勝率，具體的輸出結果如下 ：
```
Win probability for first move with random agents:
0.115 0.101 0.115
0.101 0.132 0.102
0.116 0.102 0.116
Player 1 won 585067 times
Player 2 won 288241 times
126692 ties
0.055166 seconds
18.127107 million games/sec
```
```c
static const uint32_t move_masks[9] = {
0x40040040, 0x20004000, 0x10000404, 0x04020000, 0x02002022,
0x01000200, 0x00410001, 0x00201000, 0x00100110,
};
```
至於程式碼的部份，首先可以看到位置是用 16 進位 32 bits 無號數來表示的，原因為 8 個 bits 正好可以對應到九宮格的 8 種連線方法，更具體一點的說，由最高位元開始依序往下會對應到的是九宮格中由上到下、由左到右的橫一、橫二、橫三、直一、直二、直三、右斜和左斜，而一組 nibble 中的前 3 個位元則代表了該位置是否已經被放置。以左上角第一個位置做說明，對於橫一、直一和右斜來講，就會是 100 ，其餘皆為 000 ，因此整體就會是 `0x40040040` 。
```c
static inline uint32_t is_win(uint32_t player_board)
{
return (player_board + BBBB) & 0x88888888;
}
```
此段程式碼是用來決定是否已經存在一方獲勝的狀態，也就是上述八種可能有其中一種被填滿，舉例來說，如果是橫二被填滿，那當下的狀態應為 `0x07022222` ，換句話說，被填滿的那條線必定會出現 `7(0b0111)` 這個數字，因此可以將當下的狀態加上 `0x11111111` 並和 `0x8888888` 做 `&` 運算，若有值則代表有存在著連線，否則會 `return 0` 。
`BBBB` 為 `0x11111111` 。
```c
static inline int mod7(uint32_t x)
{
x = (x >> CCCC) + (x & UINT32_C(0x7FFF));
/* Take reminder as (mod 8) by mul/shift. Since the multiplier
* was calculated using ceil() instead of floor(), it skips the
* value '7' properly.
* M <- ceil(ldexp(8/7, 29))
*/
return (int) ((x * UINT32_C(0x24924925)) >> 29);
}
```
此段程式碼想要得到 mod 7 的結果，在第一行利用了 Hacker's Delight 提出"若將 n 拆成兩部分並且將左半右移加上右半不會影響同餘的結果"的結論來進行範圍的縮小，可以看到 `x & UINT32_C(0x7FFF)` 是在計算右半邊 15 位，因此 `x >> CCCC` 就是在計算左半邊 17 位， `CCCC` 為 15 。
而這段程式碼也應用了另一定理 : $n\ \ (mod\ \ 7)\equiv\frac{8}{7}n\ \ (mod\ \ 8)$ ，因此在程式碼的最後出現了 `0x24924925` 這個數字，該數字即為 $\frac{2^{32}}{7}$ ，而再往右位移 29 位即可得到 $\frac{8}{7}$ 。
## [測驗3](https://hackmd.io/@sysprog/linux2024-quiz4#%E6%B8%AC%E9%A9%97-3)
[XTree](https://gist.github.com/jserv/7374b3e031ec3a81441b1b93b008c5d2) 為一種兼具 AVL tree 和 red-black tree 部分特性的結構，以下將說明節點刪除的部份 ：
```c
static void __xt_remove(struct xt_node **root, struct xt_node *del)
{
if (xt_right(del)) {
struct xt_node *least = xt_first(xt_right(del));
if (del == *root)
*root = least;
xt_replace_right(del, AAAA);
xt_update(root, BBBB);
return;
}
if (xt_left(del)) {
struct xt_node *most = xt_last(xt_left(del));
if (del == *root)
*root = most;
xt_replace_left(del, CCCC);
xt_update(root, DDDD);
return;
}
if (del == *root) {
*root = 0;
return;
}
/* 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(EEEE, FFFF);
}
```
在完整的程式碼中所定義的函式 ：
* `xt_first(n)` : 尋找由節點 n 所代表的子樹中最左側的節點，也就是該子樹中最小的節點。
* `xt_last(n)` : 尋找由節點 n 所代表的子樹中最右側的節點，也就是該子樹中最大的節點。
* `xt_replace_right(n, r)` : 用右子樹中的節點 r 來替換調節點 n ，並保持樹狀結構。
* `xt_replace_left(n, l)` : 用左子樹中的節點 l 來替換調節點 n ，並保持樹狀結構。
* `xt_update(n, m)` : 用來更新插入或刪除節點後的樹，利用一參數 `hint` 檢查樹是否平衡，以此決定是否需要 rotate 。
從程式碼的第一段可以看到，該程式碼會先判斷欲刪除的節點是否存在右子樹，若存在，則必須找到該子樹中最小的節點與 del 互換，並且在換完後重新調整以保持平衡，因此 `AAAA` 為 `least` ， `BBBB` 為 `xt_right(least)` 。而左子樹也是同樣道理，但由於 BST 的特性，若要交換則要找左子樹中最大的節點，因此使用了 `xt_last` 函式，其餘的部份就是將右子樹處理的方法稍做變換，因此 `CCCC` 為 `most` ， `DDDD` 為 `xt_left(most)` 。
在程式碼的最後，還需檢查整棵樹是否達到平衡，因此 `EEEE` 為 `root` ， `FFFF` 為 `parent` 。
---
## 問題記錄
* 第三周測驗二不懂為何要有 `(in | 1)` ，因為就算不強行轉為奇數即果還是正確的。
* 第四周測驗二的最後知道算出了 $\frac{8}{7}$ ，但為何這樣就同時做了 mod8 的動作。
* 在撰寫 mod9 和 mod5 的程式碼時，是由逼近的方法求解的，但數字一放大微小的差異也就跟著放大，因此我用了 if 來做調整，但好像可以不用這樣，應該從數值本身去探討其他寫法。
## 補充議題 : 浮點數乘法的實作探討
### float_mul2 探討
IEEE 754, float is single precision. Assume 32-bit
```c
float float_mul2(float x)
{
// using bitwise operation, no mul
int a = *((int*)&x);
int b = *((int*)&x);
a = (a & 0x7F800000) >> 23;
a++;
b = (b & 0x807FFFFF) | (a << 23);
return *((float*)&b);
}
```
考慮 overflow 後的做法 :
```c
float float_mul2(float x)
{
// using bitwise operation, no mul
int a = *((int*)&x);
int b = *((int*)&x);
a = (a & 0x7F800000) >> 23;
a++;
if(a == 0xFFFFFFFF)
return NAN;
else
b = (b & 0x807FFFFF) | (a << 23);
return *((float*)&b);
}
```
> https://godbolt.org/
* Instruction latency
從老師提供的 instruction tables 找到對應的處理器架構，並且將上方程式碼轉換成組合語言去逐行計算，得到的結果為 ：
```c
float_mul2: # @float_mul2
push rbp
mov rbp, rsp
movss dword ptr [rbp - 4], xmm0
mov eax, dword ptr [rbp - 4]
mov dword ptr [rbp - 8], eax
mov eax, dword ptr [rbp - 4]
mov dword ptr [rbp - 12], eax
mov eax, dword ptr [rbp - 8]
and eax, 2139095040
sar eax, 23
mov dword ptr [rbp - 8], eax
mov eax, dword ptr [rbp - 8]
add eax, 1
mov dword ptr [rbp - 8], eax
mov eax, dword ptr [rbp - 12]
and eax, -2139095041
mov ecx, dword ptr [rbp - 8]
shl ecx, 23
or eax, ecx
mov dword ptr [rbp - 12], eax
movss xmm0, dword ptr [rbp - 12] # xmm0 = mem[0],zero,zero,zero
pop rbp
ret
```
|Move instructions|Arithmetic instructions|Logic|Control transfer instructions|
|-|-|-|-|
|push * 1, mov * 12, movss * 2, pop * 1|add * 1|and * 2, sar * 1, shl * 1, or * 1|ret * 1|
計算過後可知該程式總共需要 29 (2+12+2+1+1+2+7+1+1)個 clock cycle 。
以下為稍微改進過後的版本，嘗試只用一個變數來節省記憶體空間 :
```c
float float_mul2(float x)
{
// using bitwise operation, no mul
int a = *((int*)&x);
a = ((a & 0x7F800000) + 0x00800000) | (a & 0x807FFFFF);
return *((float*)&a);
}
```
而這麼做還是太複雜，因為要先得到 exponent 的那 8 個 bits 加 1 後再把其他部分補回，因此更簡單的想法為直接在 exponent 那項加 1 ，修改過後的程式碼如下 :
```c
float float_mul2(float x)
{
// using bitwise operation, no mul
*(int*)&x += 0x00800000;
return x;
}
```
轉換為組合語言 :
```c
float_mul2: # @float_mul2
push rbp
mov rbp, rsp
movss dword ptr [rbp - 4], xmm0
mov eax, dword ptr [rbp - 4]
add eax, 8388608
mov dword ptr [rbp - 4], eax
movss xmm0, dword ptr [rbp - 4] # xmm0 = mem[0],zero,zero,zero
pop rbp
ret
```
|Move instructions|Arithmetic instructions|Logic|Control transfer instructions|
|-|-|-|-|
|push * 1, mov * 3, movss * 2, pop * 1|add * 1|none|ret * 1|
計算過後可知該程式總共需要 9 (2+3+2+1+1)個 clock cycle ，可以發現確實比第一個版本少了許多。
參考 [類神經網路的 ReLU 及其常數時間複雜度實作](https://hackmd.io/@sysprog/constant-time-relu?fbclid=IwZXh0bgNhZW0CMTAAAR0rA_gp9TtuqOBuqAJE2uutnvVZcFOSFlBPO9HspdhAA27KtGx6_rdC1aQ_aem_AbKPshOIUtxP5VUsyxPO6IG4EvaWPm2n-8Idh7E0hXijU7Y45cxOfZIWXO1McstmzjxNWThJiOkVZLV_kqQ-1Fby) 教材，裡面提到了由於對浮點數的操作成本太高，因此可以使用 union 的技巧來優化， union 可以讓浮點數和整數共用一段記憶體空間，並且在任一時刻只有一個成員有效，宣告的方法如下 ：
```c
union{
float f;
int i;
} out = {.f=x};
```
宣告一個名為 out 的 union ，裡面輸入需要操作的資料型態，而最後的 `{.f=x}` 是在做資料得初始化，將 input `x` 的值指派到 `f` ，後續則可以根據需求使用 `out.f` 及 `out.i` 來獲得兩種不同的資料型態。
以下給出一個簡單的例子 ：
```c
//code
float example(float x){
union{
float f;
int i;
} out = {.f=x};
printf("out.i = %d\n", out.i);
printf("out.f = %f\n", out.f);
}
int main()
{
float x = 2;
example(x);
return 0;
}
//output
out.i = 1073741824
out.f = 2.000000
```
原因為在這段記憶體裡存著 `0x40000000` ，以浮點數表示為 2 ，以整數表示則為 $2^{30}$ = 1073741824 。
最後，將原本的程式碼改寫 ：
```c
float float_mul2(float x)
{
// using bitwise operation, no mul
union{
float f;
int i;
} out = {.f=x};
out.i += 0x00800000;
return out.f;
}
```
轉換為組合語言 :
```c
float_mul2: # @float_mul2
push rbp
mov rbp, rsp
movss dword ptr [rbp - 4], xmm0
movss xmm0, dword ptr [rbp - 4] # xmm0 = mem[0],zero,zero,zero
movss dword ptr [rbp - 8], xmm0
mov eax, dword ptr [rbp - 8]
add eax, 8388608
mov dword ptr [rbp - 8], eax
movss xmm0, dword ptr [rbp - 8] # xmm0 = mem[0],zero,zero,zero
pop rbp
ret
```
|Move instructions|Arithmetic instructions|Logic|Control transfer instructions|
|-|-|-|-|
|push * 1, mov * 3, movss * 4, pop * 1|add * 1|none|ret * 1|
計算過後可知該程式總共需要 11 (2+3+4+1+1)個 clock cycle ，雖然相較於第二個版本多了 2 個 clock cycle ，但用 union 的寫法會比較安全，因為強制將資料型別轉換可能會造成 strict aliasing ，在 C 語言中會發生為定義的行為。
### float_mul_power_of_2 探討
```c
float float_mul_power_of_2(float x, int e)
{
// using bitwise operation, no mul
union{
float f;
int i;
} out = {.f=x};
out.i += (e << 23);
return out.f;
}
```
利用和上述 `float_mul2` 相同的概念，若要乘上 2 的冪，則把指數部分的 `+1` 改為 `+e` 即可，但若是直接用 `(e << 23)` 這樣的寫法會有問題 :
若該整數 `e` 夠大或夠小，則位移後會影響到最左方的 sign bit，因此我們需想辦法將其大小控制在 8 個 bits 以內，以符合 IEEE 754 的範圍。
解決方法 :
首先，由於浮點數的指數部分的範圍為 0 ~ 255 ，因此我們可以把指數部分單獨取出來看是否有在範圍內，對應的程式碼如下 :
```c
int exp = (out.i >> 23) & 0xFF;
exp += e;
if(exp < 0 || exp > 255)
return 0;
```
該程式碼的作用為先將 fraction bits 23 位過濾掉，然後再和 8 個 1 做 `&` 運算，以取得原數的指數部分，接著用整數型態的變數和 input `e` 相加，最後再判斷是否有在範圍內，完整的程式碼如下 :
```c
float float_mul_power_of_2(float x, int e)
{
// using bitwise operation, no mul
union{
float f;
int i;
} out = {.f=x};
int exp = (out.i >> 23) & 0xFF;
exp += e;
if(exp < 0 || exp > 255)
return 0;
else
out.i = (out.i & 0x807FFFFF) | (exp << 23);
return out.f;
}
```
其中 `return 0` 的部分若是引進 `math.h` 函式庫則可以再細分為不同的情況來表達，以下為修改過後的程式碼 :
```c
float float_mul_power_of_2(float x, int e)
{
// using bitwise operation, no mul
union{
float f;
int i;
} out = {.f=x};
int exp = (out.i >> 23) & 0xFF;
exp += e;
if(exp < 0)
return 0;
else if(exp >= 255)
return INFINITY;
else
out.i = (out.i & 0x807FFFFF) | (exp << 23);
return out.f;
}
```
Sign in with Wallet
Connect another wallet