# 2024q1 Homework4 (quiz3+4)
**contributed by <`csotaku0926`>**
## 第三週
### 測驗一
- [ ] 原理
根據 [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+...+a_n)^2,a_m=2^m (or 2^m=0), m \in \{0,1,...,n\}
$$
已知 $(a+b)^2=a^2+2ab+b^2$，則有以下遞迴關係
(單純看式子很難懂，可以用平方矩陣去想比較好理解)
$$
\begin{bmatrix}
a_0a_0 & a_0a_1 & a_0a_2\\
a_1a_0 & a_1a_1 & a_1a_2\\
a_2a_0 & a_2a_1 & a_2a_2\\
\end{bmatrix}
$$
$N^2=((a_0+a_1+...+a_{n-1})+a_n)^2=a_n^2+2a_n\sum_{k=0}^{n-1}a_k+(a_0+a_1+...+a_{n-1})^2$
$=a_n^2+(2P_n+a_{n-1})a_{n-1}+...+(2P_1+a_0)a_0$
令平方根逼近值$P_m=\sum_{i=m}^{n}a_i$，$m\in \{0,1,...,n\}$
$P_0=\sum_{i=0}^{n}a_i$ 即是欲求平方根
顯然 $P_{m-1}=P_m+a_{m-1}$
因此對於每個 $m\in\{n,...,1,0\}$，於$P_{m-1}$ 加上 $a_m$，由於是 2 的冪，因此只有以下兩種可能:
$P_{m-1}=P_m+a_{m-1},(P_m^2<N^2)$,
$P_{m-1}=P_m,(otherwise)$
但我們的程式並不是用這個較直白的方法判別是否加 $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$
想了很久，原來只是單純將加法的兩項拆開來看，即:
$c_m=P_{m+1}2^{m+1}$, $d_m=(a_{m})^2$
$Y_m=c_m+d_m$
每次迴圈更新:
$\begin{gather*}
c_{m-1}=P_m2^m=(P_{m+1}+a_m)2^m=\begin{cases}
c_m/2+d_m \ \ if \ \ a_m = 2^m\\
c_m/2 \ \ if \ \ a_m=0
\end{cases}
\end{gather*}$
$d_{m-1}=d_m/4$
計算到最後 $c_{-1}=P_0=N$, $d_0=0$
太難理解了，單單根號運算為了追求表現，竟然會變得這麼複雜 (我資質駑鈍)
- [ ] 以 [ffs](https://man7.org/linux/man-pages/man3/ffs.3.html) 改寫
> [commit e413059](https://github.com/csotaku0926/Linux_Hw4_Lab)
原先 `(31 - __builtin_clz(x))` 的部份是為了算 MSB，現在用 `strings.h` 提供的 `ffs` 達成
- [ ] Linux kernel 的平方根運算
要理解 linux 核心程式碼並不是容易的事情，光是說明提供參考的 [block](https://github.com/torvalds/linux/tree/master/block) 目錄就有許多檔案，要如何在其中找到平方根實作？
透過 google 關鍵字搜尋，總算是在 `block/blk-wbt.c` 這個檔案中找到有關平方根的動作:
```
- If the minimum latency in the above window exceeds some target, increment
* scaling step and scale down queue depth by a factor of 2x. The monitoring
* window is then shrunk to 100 / sqrt(scaling step + 1).
```
由註解看出這是一段有關調整 monitor window size 的程式碼，有點像網路程式設計學到的滑動視窗 (sliding window) ，可以根據 scaling step 逐漸縮小
實際呼叫平方根的函式為 `rwb_arm_timer` :
```c
static void rwb_arm_timer(struct rq_wb *rwb)
{
struct rq_depth *rqd = &rwb->rq_depth;
if (rqd->scale_step > 0) {
/*
* We should speed this up, using some variant of a fast
* integer inverse square root calculation. Since we only do
* this for every window expiration, it's not a huge deal,
* though.
*/
rwb->cur_win_nsec = div_u64(rwb->win_nsec << 4,
int_sqrt((rqd->scale_step + 1) << 8));
```
而定義 `int_sqrt` 的函式位於 `/lib/math/int_sqrt.c`
```c
unsigned long int_sqrt(unsigned long x)
{
unsigned long b, m, y = 0;
if (x <= 1)
return x;
m = 1UL << (__fls(x) & ~1UL);
while (m != 0) {
b = y + m;
y >>= 1;
if (x >= b) {
x -= b;
y += m;
}
m >>= 2;
}
return y;
}
EXPORT_SYMBOL(int_sqrt);
```
除了獲取 MBS 的方式是使用 `_fls`, 將迴圈改為 while loop 之外，整體演算法是相同的
- [ ] 閱讀論文與改進
- [論文1](https://web.ece.ucsb.edu/~parhami/pubs_folder/parh99-asilo-table-sq-root.pdf)
- [論文2](https://web.archive.org/web/20210208132927/http://assemblyrequired.crashworks.org/timing-square-root/)
### 測驗二
- [ ] 程式解讀
> 根據[教師所言](https://www.facebook.com/groups/system.software2024/?multi_permalinks=1572555173522375¬if_id=1711546821881676¬if_t=feedback_reaction_generic&ref=notif)，這題目是希望引導學員思考針對 RV32I/RV64I 這類沒有整數乘法和除法指令的指令集，該如何撰寫有效的程式
希望以 bitwise operation 實作除法，但若除數 (這個例子中, 10) 包含非 2 的冪的因數，會有不精確的問題
[項目說明](https://hackmd.io/@sysprog/linux2024-quiz3) 中透過推導證明，若存在 $n\in N$，選定一 $x$ 使 $\frac{n}{x}$ 直到小數第一位都是精準的，即 $a.b\le \frac{n}{x} \le a.b9$
則對於任何小於 n 的非負整數 l ，$\frac{l}{x}$ 直到小數第一位也是精準的
考慮輸入不大於 19 的狀況下，找出能夠維持小數第一位精確度的除數 x :
$1.9\le \frac{19}{x}\le1.99$ --> $9.55\le x\le10$
由於是 bitwise operation，除數的形式應該要是 $\frac{a}{2^n}$，a 是任意配對的正整數
接下來在測試程式處，有個匪夷所思的地方:
```c
q = ((((n >> 3) + (n >> 1) + n) << 3) + d0 + d1 + d2) >> 7;
```
這段程式碼想要求得 `13 * tmp` ，但他的步驟是先求得 `13 * tmp / 8` 後透過左移變回 `13 * tmp`，又此時因經過右移導致部份位元遭捨棄，因此又用 d0, d1, d2 補回。
為什麼不直接求 `13 * tmp` 就好了？ 可以進一步寫成以下形式，節省變數運用:
```c
q = ((n << 3) + (n << 2) + n) >> 7;
```
可見測試題中的程式碼只是其中一種可行方法，不代表最優解
而包裝程式又換了一種寫法:
```c
uint32_t x = (in | 1) - (in >> 2); /* div = in/10 ==> div = 0.75*in/8 */
uint32_t q = (x >> 4) + x;
```
這一行初步認為是將除數 x 換成 $0.75\times in$，而 q 的值指派為 $\frac{51\times in}{64}=0.796\times in$
接下來又是一個匪夷所思的部份:
```c
q = (q >> 8) + x;
```
這行程式碼重複了四次，也就是說，in 不斷加入 `in >> 8` 數值，若 in 是小於 256 的正整數，這些程式對 in 毫無作用。也許要結合後段程式碼才會知道這樣寫的目的
最後，q 需要再右移一個數，才是最終的商。考慮
$1.9\le \frac{19}{x}\le1.99$ --> $9.55\le x\le10$
需要將 $\frac{64}{51}$ 乘上適當的 2 的冪，使其符合 x 的條件
而 $\frac{64\times 8}{51}=10.03$，最符合上述條件，可知 `*div = (q >> 3);`
最後一行是在算餘式，已知 $div=q>>3$，$mod=n-div\times10$
`q & ~0x7` 相當於 $div \times 8$
所以: `*mod = in - ((q & ~0x7) + (*div << 1));`
- [ ] instruction cycle 計算
> [項目參考](https://www.agner.org/optimize/instruction_tables.pdf)
> 一行行看指令對照表有些費時，如果可以做個程式自動計算 cycle 應該很不錯
想要計算這部份，第一步是要確認電腦的處理器版本，我使用的型號為 Intel 12th gen Core
```
Model name: 12th Gen Intel(R) Core(TM) i5-12450H
```
對照圖表，並沒有發現直接對應的微處理器，先以型號最接近的 "Nehalem" 探討
先以 `gcc div10.c -o div10` 產出執行檔後，以 `objdump -d div10 > div10.asm` 觀察反組譯的結果 (-O 是 gcc 最佳化選項，<s>預設是 -O0</s> )
> 預設最佳化等級取決於 GNU toolchain 編譯時的選項 :notes: jserv
首先是以除法和模數寫成的直覺版本:
```c
void div10_naive(uint32_t n, uint32_t *q, uint32_t *r)
{
*q = n / 10;
*r = n % 10;
}
```
實際查看指令表格，有許多欄位。觀看敘述得知 "latency" 以 cpu cycle 為時間單位測量指令延遲，因此要對照這個欄位，以下是出現在反組譯檔案中的指令之對應 cycle 數:
| 指令 | cycle | count |
| --------- | ----- | ----- |
| push r | 3 | 1 |
| mov r r/i | 1 | 18 |
| imul | 3 | 2 |
| add/sub | 1 | 3 |
| shr/shl | 1 | 4 |
| pop r | 2 | 1 |
佔用 CPU 週期總計 36 個
至於只使用位元運算的除法:
```c
void div10(uint32_t n, uint32_t *q, uint32_t *r)
{
*q = ((n << 3) + (n << 2) + n) >> 7; /* n / (128/13) */
*r = n - (((*q << 2) + *q) << 1);
}
```
| 指令 | cycle | count |
| --------- | ----- | ----- |
| push r | 3 | 1 |
| mov r r/i | 1 | 18 |
| lea | 1 | 3 |
| add/sub | 1 | 4 |
| shr/shl | 1 | 2 |
| pop r | 2 | 1 |
佔用 CPU 週期總計 32 個，較先前版本略為減少，發現 `imul` 被換成 `lea` 指令
參考 [var-x 的 writeup](https://hackmd.io/@vax-r/linux2024-homework4#%E6%B8%AC%E9%A9%97%E4%BA%8C) ，可以利用 `perf stat <command>` 這個命令測量 cpu cycle
- 註: 第一次使用 `perf stat` 會遇到權限問題，需要修改 `/proc/sys/kernel/perf_event_paranoid`
- 可以利用 `sudo sysctl <parameter>=<value>` 快速修改，[ref](https://askubuntu.com/questions/1471162/how-to-change-kernel-perf-event-paranoid-settings)
於同個 c 檔案定義好兩個除法函式後，以 `main()` 呼叫之。編譯出執行檔後測試 `perf stat ./div10` :
```c
int main(void)
{
uint32_t r, q;
div10(123, &q, &r);
}
```
初始版本的除法耗費 977,489 個 cycle，而不用除法模數的版本只花費 880,667 個，這是很大的進步
- [ ] 撰寫 modulo 5, modulo 9 (不依賴除法指令)
注意先前的 div10 指令，預期輸入是不比 19 大的正整數，但現在有效範圍必須是 `uint32_t` 以內
[Modulus without Division, a tutorial](https://homepage.cs.uiowa.edu/~dwjones/bcd/mod.shtml#exmod3) 一文針對模數處理，提出了不需除法指令的實作
> [commit fd9cce](https://github.com/csotaku0926/Linux_Hw4_Lab/commit/fd9cce61339ec323019661786e6e7bf8ce70d659)
考慮 modulo 5，已知 $8\mod{5}=3$，因此符合以下條件:
$a\mod{5}=(\frac{3a}{8}+a\mod{8})\mod{5}$
可以用這個式子，加上迴圈判斷數值是否還需要做減去的動作。
這個原則可以運用於任一正整數模數，即:
$$
a\mod{m}=((2^n\mod{m})\frac{a}{2^n}+a\mod{2^n})\mod{m}
$$
其中 m 為任一正整數，n 為非負整數
但是這麼做的時間成本還是太高，其中一種作法是用 mod 15 推算出 mod 5
$$
a\mod{5}=(a\mod{15})\mod{5}
$$
因為 mod 15 存在快速用位元運算計算的方式
> [commit b408910](https://github.com/csotaku0926/Linux_Hw4_Lab/commit/b408910659ce8ffcfa86df06c7682e6be668c231)
至於 mod 9，目前想不到較為直接的方式計算
我想到的一種方法為先計算 mod 63，再推算出 mod 9
### 測驗三
程式碼回傳 log2 數值
顯然，只須回傳最高位數 (MSB) 是第幾位即可
對於 32 位元的無號整數，則可以透過 GNU extension 的 `__builtin_clz()` 達成
#### linux 核心 log2 案例
透過 google 搜尋，找到 [bootlin](https://elixir.bootlin.com/linux/latest/source/include/linux/log2.h) 這個網站，裡面記載了許多核心的程式碼
其中可以找到位於核心目錄的 `/include/linux/log2.h`
針對 arch 架構，允許透過在 `asm/bitops.h` 撰寫比使用 `fls` 更有效率的 log2 實作
這也解釋為什麼許多標頭檔都會加上 `#ifndef...endif` 這類的巨集
```c
#ifndef CONFIG_ARCH_HAS_ILOG2_U32
static __always_inline __attribute__((const))
int __ilog2_u32(u32 n)
{
return fls(n) - 1;
}
#endif
```
可以發現具有常數 log2 的實作方式
```c
#define const_ilog2(n) \
( \
__builtin_constant_p(n) ? ( \
(n) < 2 ? 0 : \
(n) & (1ULL << 63) ? 63 : \
(n) & (1ULL << 62) ? 62 : \
// ...
(n) & (1ULL << 2) ? 2 : \
1) : \
-1)
```
根據 [gnu 官網](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html)，`__bulitin_constant_p()` 可用於檢驗表達式在編譯時間是否為常數
參考 [編譯器和最佳化原理篇](https://hackmd.io/@sysprog/c-compiler-optimization#From-Source-to-Binary-How-A-Compiler-Works-GNU-Toolchain)，在編譯器最佳化的策略上，經常可以看到將迴圈內容展開的步驟，稱作 loop unrolling
也有提供非常數的實作:
```c
#define ilog2(n) \
( \
__builtin_constant_p(n) ? \
((n) < 2 ? 0 : \
63 - __builtin_clzll(n)) : \
(sizeof(n) <= 4) ? \
__ilog2_u32(n) : \
__ilog2_u64(n) \
)
```
原來在 linux 核心追求的程式設計，即使是像 `log2` 這樣基本的函式，不只有正確性的考量，也需要考量硬體配置的差異 (如 32 位元及 64 位元)
### 測驗四
根據 EWMA 公式:
$\begin{gather*}
S_t=\begin{cases}
Y_0 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ t=0\\
\alpha Y_t+(1-\alpha)S_{t-1} \ \ t>0
\end{cases}
\end{gather*}$
- $\alpha$ 為資料遞減的程度 ($0\le\alpha<1$)
- $Y_t$ 為在時間為 t 時的資料數值
- $S_t$ 為在時間為 t 時的 EWMA
對應到程式碼，有以下架構:
```c
struct ewma {
unsigned long internal;
unsigned long factor;
unsigned long weight;
};
```
`internal` 對應到 $S_t$，
由於平均數涉及到小數點，這裡使用 `factor` 表示定點數，`weight` 則對應 $\alpha$
於是 `ewma_add` 對應到上方公式
```c
avg->internal = avg->internal
? (((avg->internal << avg->weight) - avg->internal) +
(val << avg->factor)) >> avg->weight
: (val << avg->factor);
```
#### linux 核心案例
### 測驗五
TODO
## 第四週
Sign in with Wallet
Connect another wallet