# 2024q1 Homework4 (quiz3+4)
contributed by <[`hungyuhang`](https://github.com/hungyuhang)>
## 第 3 週測驗題
### [測驗 `1`](https://hackmd.io/@sysprog/linux2024-quiz3#%E6%B8%AC%E9%A9%97-1)
這個測驗內容是 [Digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) ，該演算法的目的是去找出某數字開平方，以下是簡單的運行方式：
將要開平方的數（也就是 $x$ ）拆成 2 的冪相加 (若是十進位則為 10 的冪)，也就是如下形式：$$ x = (000b_0b_1b_2\cdots b_{n-1}b_{n})^2 $$
其中 $$ 000b_0b_1b_2\cdots b_{n-1}b_{n} $$ 為 bits pattern，其中 $b_0$ 為最高位的 1 。
演算法的邏輯，就是從 $b_0$ 開始，一個一個的去判斷 $b_i$ 是 1 還是 0：
```
for i = [0 to n]:
determine if bi is 1 or 0
```
在版本一跟版本二的程式碼當中，每判斷一個位元的數值，就需要做一次相乘並且比較（判斷 $N^2 - P_m^2$ 大於還是小於零），而這些操作對效能會帶來負面影響。
但其實「每一個位元要是 0 還是 1 」這樣的資訊是可以透過上一次迭代的計算結果去算出來的，不需要在每一次迭代再重新做一次 $N^2 - P_m^2$ 這種計算成本很高的運算。所以在版本三的程式碼當中，就透過上一次迭代所產生的結果去計算每個 $b_i$ 到底是 1 還是 0 。
以下是版本三的程式碼：
```diff
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) {
+ for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
int b = z + m;
- z >>= BBBB;
+ z >>= 1;
if (x >= b)
x -= b, z += m;
}
return z;
}
```
其中各個變數代表的涵義如下（參考 [quiz3 測驗1](https://hackmd.io/@sysprog/linux2024-quiz3#%E6%B8%AC%E9%A9%97-1) 的敘述）：
- z 代表 $c_m$
- m 代表 $d_m$
- b 代表假設 $a_m = 2^m$ 的話， $Y_m$ 會等於的值
- x 代表 $X_{m+1}$
而 `if (x >= b)` 這行程式碼就是在判斷 $a_m$ 實際上是 $2^m$ 還是 $0$ ，如果 $a_m$ 真的是 $2^m$ 的話才會執行判斷式裡面的程式碼。
#### 使用 `ffs` / `fls` 取代 `__builtin_clz`
依照程式碼邏輯，應使用 `fls` 函式來取代原有的 `__builtin_clz` 函式。
首先我透過 `man` 命令去尋找 `fls` 函式的說明文件，但是沒有找到，只有找到 `ffs` 的函式說明。
實際將 `fls` 寫入程式碼裡面，並且引入 `<strings.h>` 標頭檔（這個標頭檔定義了 `ffs` 函式），程式碼也會因為找不到 `fls` 函式而無法成功編譯。所以在這裡我仍然使用 `__builtin_clz` 這個函式。
#### 實作無分支版本的程式碼
在這裡我嘗試將 `for` 迴圈內的 `if` 給去除，實作方式如下：
```diff
- if (x >= b)
- x -= b, z += m;
+ int r = x >= b;
+ z += m * r;
+ x -= b * r;
```
將兩個版本的 assembly code 做比較，無分支版本的程式碼的確少了一個 `jl` (jump if less) 指令。
使用 perf 對兩個版本的程式碼做測量，測量方式為對 1000000 開平方，並且重複測量一百次再取平均的測量數據作為結果：
- 有分支版本
```
Performance counter stats for './quiz3-1' (100 runs):
0.33 msec task-clock # 0.523 CPUs utilized ( +- 4.13% )
0 context-switches # 0.000 /sec
0 cpu-migrations # 0.000 /sec
59 page-faults # 178.352 K/sec ( +- 0.14% )
101,4122 cycles # 3.066 GHz ( +- 0.90% )
111,7387 instructions # 1.10 insn per cycle ( +- 0.09% )
20,1433 branches # 608.914 M/sec ( +- 0.08% )
5579 branch-misses # 2.77% of all branches ( +- 0.60% )
0.0006320 +- 0.0000245 seconds time elapsed ( +- 3.88% )
```
- 無分支版本
```
Performance counter stats for './quiz3-1' (100 runs):
0.32 msec task-clock # 0.520 CPUs utilized ( +- 2.55% )
0 context-switches # 0.000 /sec
0 cpu-migrations # 0.000 /sec
58 page-faults # 180.272 K/sec ( +- 0.15% )
99,9502 cycles # 3.107 GHz ( +- 0.81% )
111,1421 instructions # 1.11 insn per cycle ( +- 0.08% )
20,0526 branches # 623.263 M/sec ( +- 0.06% )
5485 branch-misses # 2.74% of all branches ( +- 0.68% )
0.0006186 +- 0.0000163 seconds time elapsed ( +- 2.63% )
```
從結果來看，兩個版本的程式碼的效能，甚至 branch 的數量幾乎一模一樣。目前的猜想是編譯器可能已經把有分支的版本做過一些優化了。
### 測驗 `2`
### 測驗 `3`
如果要計算以 2 為底的對數，且輸入跟輸出都是整數的話，那麼問題的解法就可以被簡化成下面兩個步驟。假設我們要求 18 以二為底的對數：
1. 將 18 轉成二進位表示式：`00010010`
2. 從右邊開始數（從零開始），找到最左邊的 `1` 是在第幾位，以 18 為例的話就是第 4 位，這個 4 就是答案。換成數學式的表示法就是 $\lfloor log_2(18) \rfloor = 4$ 。
依照這個邏輯，就可以對版本二的程式碼進行填空：
```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` 的版本：
```diff
int ilog32(uint32_t v)
{
- return (31 - __builtin_clz(DDDD));
+ return (31 - __builtin_clz(v));
}
```
### 測驗 `4`
以下是程式碼的填空：
```diff
struct ewma *ewma_add(struct ewma *avg, unsigned long val)
{
avg->internal = avg->internal
- ? (((avg->internal << EEEE) - avg->internal) +
+ ? (((avg->internal << avg->weight) - avg->internal) +
- (val << FFFF)) >> avg->weight
+ (val << avg->factor)) >> avg->weight
: (val << avg->factor);
return avg;
}
```
觀察程式碼，邏輯上我們大可將下面的式子（這邊稱作 `expr1` ）
```c
(((avg->internal << avg->weight) - avg->internal) +
(val << avg->factor)) >> avg->weight
```
改寫成下面這個樣子（這邊稱作 `expr2` ）
```c
(avg->internal - (avg->internal >> avg->weight)) +
((val << avg->factor) >> avg->weight)
```
其中， `expr2` 的邏輯也會比較接近 EWMA 的公式： $\alpha Y_t + (1 - \alpha)\cdot S_{t-1}$ 。
但是為什麼不使用 `expr2` 呢？
觀察 `expr2` 的組成，是先對下面兩個子算式進行運算完之後，才進行加總：
- `(avg->internal - (avg->internal >> avg->weight))`
- `((val << avg->factor) >> avg->weight)`
但是上面兩個算式都含有 `>>` 這個運算元，也就是說這兩個算式所算出的結果都有可能跟真實的答案是有所誤差的（比如說把 3 `>>` 一個位元的話會得到 1 ，而不是 1.5 ），而當我們把兩個有誤差的結果加起來之後，可能出現的誤差值只會變得更大。
現在我們來觀察 `expr1` 的算法，可以看到它先把 $(1 - \alpha)\cdot S_{t-1}$ 除以 $\alpha$ 倍，得到 $(S_{t-1} / \alpha) - S_{t-1}$ ，並且加上 $Y_t$ 之後 ，才把加總的結果乘上 $\alpha$ 倍。這麼一來我們在運算時，只有在最後乘以 $\alpha$ 的時候可能會產生誤差，相較於 `expr2` 來說，這個方式所造成的誤差會比較小。
### 測驗 `5`
```diff
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 > GGG) + 1;
+ return (r | shift | x > 1) + 1;
}
```
這個程式的運作原理其實跟測驗 3 的 `ilog2` 一樣，都是透過找到輸入的數字最左邊的 set bit 的位置來推出 $\lfloor log_2(x) \rfloor$ 是多少。但因為這個程式求的是 $\lceil log_2(x) \rceil$ ，所以在程式碼的最後會再把結果加上 1 來得到最後的輸出。
但是這個方法當遇到 $\lfloor log_2(x) \rfloor$ = $\lceil log_2(x) \rceil$ ，也就是 x 等於二的冪的時候就會出問題了，所以這個程式就在開頭把輸入的數字減一來繞過這個問題。
比如說當輸入 $2^5$ 的話，如果我們不先把輸入的數字減一的話，得到的輸出就會是 6 ，這顯然不是正確答案，但如果我們先把 $2^5$ 減去 1 再去做計算的話，就可以的到正確的答案，也就是 5 。
#### 程式碼改善
嘗試進一步改進程式碼，讓該函式能夠輸入 0 而不會出現錯誤。
可以觀察到在原來的程式碼中，如果輸入 0 的話，程式碼反而會變成在算 $log_2$(0xFFFFFFFF) ，這顯然不是正確的邏輯。
有一個簡單的方式可以改善程式碼，就是當輸入等於 0 的時候，就不要把輸入減一：
```c
if (x != 0)
x--;
```
但作業要求改善的方式也要是 branchless 的，所以進一步再把程式碼改成如下：
```c
x = x - !!x;
```
其中 `!!x` 只有在 x 不等於 0 的時候會等於 1 ，當 x 等於 0 的時候 `!!x` 也會等於 0 ，這麼一來就可以在輸入等於 0 的時候不做減一的動作了。
## 第 4 週測驗題
### 測驗 `1`
```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;
}
```
這個程式的目的，是去找出整數陣列 `nums` 中，所有整數對之間的 Hamming distance 的總和。
首先雙層 for 迴圈會計算陣列中兩兩之間的 Hamming distance 並且做加總。以下圖為例，假設 `numSize` 為 3 ，程式會在雙層迴圈中把下面九個 Hamming distance （ $d_{ij},\ 0 \le i \le 2,\ 0 \le j \le 2$ ）加總起來：
```
nums[0] nums[1] nums[2]
nums[0] d00 d01 d02
nums[1] d10 d11 d12
nums[2] d20 d21 d22
```
因為自己跟自己的 Hamming distance 必定為零，所以 $d_{00}$, $d_{11}$, $d_{22}$ 的值為零，而從剩下的值當中，我們可以發現 $d_{ij}$ 跟 $d_{ji}$ 其實都是在算同一個數對之間的 Hamming distance ，換句話說，這個雙層迴圈所算出來的值應該要除以 2 才會是最後的結果。所以最後回傳的值應該要往右移一位才會是正確的結果。
### 測驗 `2`
### 測驗 `3`
Sign in with Wallet
Connect another wallet