# 2024q1 Homework4 (quiz3+4)
contributed by < `yeh-sudo` >
## 第三周測驗題
### 測驗 `1`
#### 解釋程式碼運作原理
假設要求的平方根為 $x$ ，並且將 $x^2$ 表示為以下形式：
$$x^2=(000b_0b_1b_2...b_{n-1}b_n)^2$$
其中 $000b_0b_1b_2...b_{n-1}b_n$ 為 bit pattern ， $b_0$ 為最高位的 1 。所以 $x^2$ 可以寫成：
$$x^2=(b_02^n+b_12^{n-1}+b_22^{n-2}+...+b_{n-1}2^1+b_n2^0)^2$$
設 $a_0=b_02^n,a_1=b_12^{n-1},a_2=b_22^{n-2}...a_{n-1}=b_{n-1}2^1,a_n=b_n2^0$， $x^2$ 可以表示成：
$$x^2=(a_0+a_1+a_2+...+a_{n-1}+a_n)^2$$
根據 [Digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) ：
\begin{align*}
x^2&=(a_n+a_{n-1}+a_{n-2}+...+a_1+a_0)^2 \\
&=a_n^2+2a_na_{n-1}+a_{n-1}^2+2(a_n+a_{n-1})a_{n-2}+a_{n-2}^2+...+a_1^2+2(\sum_{i=1}^{n} a_i)a_0+a_0^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}^{n} a_i)+a_0]a_0
\end{align*}
設 $P_m=\sum_{i=m}^{n} a_i$ ，則 $x^2$ 可以表示為：
$$x^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_{n-2}+...+a_1+a_0$ 即為要求出的平方根。所以，從 $m=n$ 一路試到 $m=0$ ，測試 $P_m^2\le x^2$ 是否成立，分為兩種情況：
\begin{cases}
P_m=P_{m+1}+a_m & \quad \text{if } P_m^2\le x^2\\
P_m=P_{m+1} & \quad \text{otherwise}
\end{cases}
在檢驗時，如果每次都要計算 $P_m$ ，會需要很高的成本，所以可以利用重複利用每一次計算的差值 $X$ ，可以得出：
\begin{align*}
X_n&=x^2-P_n^2 \\
&=x^2-a_n^2
\end{align*}
\begin{align*}
X_{n-1}&=x^2-P_{n-1}^2 \\
&=x^2-(a_n+a_{n-1})^2 \\
&=x^2-a_n^2-a_{n-1}^2-2a_na_{n-1} \\
&=(x^2-a_n^2)-a_{n-1}^2-2a_na_{n-1} \\
&=X_n-a_{n-1}^2-2a_na_{n-1} \\
&=X_n-[2a_n+a_{n-1}]a_{n-1} \\
&=X_n-[2P_n+a_{n-1}]a_{n-1}
\end{align*}
\begin{align*}
X_{n-2}&=x^2-P_{n-2}^2 \\
&=x^2-(a_n+a_{n-1}+a_{n-2})^2 \\
&=x^2-a_n^2-a_{n-1}^2-a_{n-2}^2-2a_na_{n-1}-2a_na_{n-2}-2a_{n-1}a_{n-2} \\
&=X_n-[2P_n+a_{n-1}]a_{n-1}-2a_na_{n-2}-2a_{n-1}a_{n-2}-a_{n-2}^2 \\
&=X_{n-1}-[2(a_n+a_{n-1})+a_{n-2}]a_{n-2} \\
&=X_{n-1}-[2P_{n-1}+a_{n-2}]a_{n-2}
\end{align*}
\begin{align*}
...
\end{align*}
\begin{align*}
X_{m}&=X_{m+1}-[2P_{m+1}+a_m]a_m
\end{align*}
設 $Y_m=[2P_{m+1}+a_m]a_m$ ，則可以寫成以下遞迴關係：
$$X_{m}=X_{m+1}-Y_m$$
將 $Y_m$ 拆成 $c_m$ 與 $d_m$ 兩部份，
\begin{align*}
Y_m&=[2P_{m+1}+a_m]a_m \\
&=2P_{m+1}a_m+a_m^2 \\
&=2^{m+1}P_{m+1}+{(2^m)}^2 \\
&=c_m+d_m
\end{align*}
\begin{cases}
Y_m=c_m+d_m & \quad \text{if } a_m=2^m \\
Y_m=0 & \quad \text{otherwise}
\end{cases}
\begin{cases}
c_m=2^{m+1}P_{m+1} \\
d_m={(2^m)}^2
\end{cases}
可以從 $c_m$ 與 $d_m$ 推出下一輪的 $c_{m-1}$ 與 $d_{m-1}$ ，再用 $c_{m-1}$ 與 $d_{m-1}$ 推出 $c_{m-2}$ 與 $d_{m-2}$ ，直到推出 $c_{-1}=P_0$ 即為所求的平方根 $x$ ，而首項 $c_n=0, d_n={(2^n)}^2$ 。
\begin{align*}
c_{m-1}&=2^mP_m \\
d_{m-1}&={(2^{m-1})}^2={(2^m)}^22^{-2}=d_m/4
\end{align*}
此處的 $P_m=P_{m+1}+a_m$ 成立的條件為 $P_m^2\le x^2$ ，若不成立，則 $P_m=P_{m+1}$ ，所以當 $P_m^2\le x^2$ 時， $a_m=2^m$ ，反之則 $a_m=0$ 。所以最後可以將 $c_m$寫成以下關係：
\begin{gather*}
c_{m-1}=
\begin{cases}
c_m/2+d_m & \quad \text{if } a_m=2^m \\
c_m/2 & \quad \text{if } a_m=0
\end{cases}
\end{gather*}
在程式碼中， $c_m$ 對應到變數 `z` ， $d_m$ 對應到變數 `m` ， 差值 $X$ 對應到變數 `x` 。由於是紀錄差值，所以每當差值大於等於 $c_m+d_m$ ，也就是 `b = z + m` 時，就將差值扣掉 `b` ，並且，因為差值為正，代表 $P_m^2\le x^2$ ，也就是 $a_m=2^m$ ，所以 $c_{m-1}=c_m/2+d_m$ ，對應到的程式碼操作為 `z += m` ，持續進行迴圈直到 `m` 為 0 時，結束迴圈。
```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;
}
```
#### 以 ffs / fls 取代 __builtin_clz
程式碼 `(31 - __builtin_clz(x))` 的作用是找到最大的 bit 所在的位置，而使用第二週測驗題提供的 `__ffs` ，可以實作出找出對大 bit 的 `fls` 函式，實作方法是當 `x` 不為 0 時，持續尋找 `ffs` ，找到時就使用測驗題提供的 `__clear_bit` 將找到的 bit 清為 0 ，當 `x` 為 0 時，代表最後找到的 `ffs` 為最大的 bit 。
```c
static inline unsigned long fls(unsigned long x)
{
int result = 0;
while (x) {
unsigned int bit = __ffs(x);
if (bit > result)
result = bit;
__clear_bit(bit, &x);
}
return result;
}
```
另外還有一個方法可以找出 `fls` ，將 `x` 的逐位元反轉，再用 31 減掉 `ffs` 即可得到 `fls` 。
```c
static inline unsigned long fls(unsigned long x)
{
unsigned long result = x;
result = ((result & 0xffff0000) >> 16) | ((result & 0x0000ffff) << 16);
result = ((result & 0xff00ff00) >> 8) | ((result & 0x00ff00ff) << 8);
result = ((result & 0xf0f0f0f0) >> 4) | ((result & 0x0f0f0f0f) << 4);
result = ((result & 0xcccccccc) >> 2) | ((result & 0x33333333) << 2);
result = ((result & 0xaaaaaaaa) >> 1) | ((result & 0x55555555) << 1);
return 31 - __ffs(result);
}
```
#### 提供分支和無分支 (branchless) 的實作
在教材〈[解讀計算機編碼](https://hackmd.io/@sysprog/binary-representation#%E4%B8%8D%E9%9C%80%E8%A6%81%E5%88%86%E6%94%AF%E7%9A%84%E8%A8%AD%E8%A8%88)〉中，有描述不需要分支的設計，當 `INT32_MIN <= (b - a) <= INT32_MAX` 條件成立時，可以使用 bitwise 操作改寫分支。
- [ ] 改寫 `i_sqrt`
```diff
- if (x >= b)
- x -= b, z += m;
+ int tmp = x;
+ x -= (~(tmp - b) >> 31) & b;
+ z += (~(tmp - b) >> 31) & m;
```
- [ ] 改寫 `__ffs`
```diff
#if BITS_PER_LONG == 64
- if ((word & 0xffffffff) == 0) {
- num += 32;
- word >>= 32;
- }
+ num += !(word & 0xffffffff) & 32;
+ word >>= 32;
#endif
- 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;
+ num += !(word & 0xffff) & 16;
+ word >>= 16;
+ num += !(word & 0xff) & 8;
+ word >>= 8;
+ num += !(word & 0xf) & 4;
+ word >>= 4;
+ num += !(word & 0x3) & 2;
+ word >>= 2;
+ num += !(word & 0x1) & 1;
return num;
```
#### Linux 核心原始程式碼
```shell
~/linux/block master ❯ grep sqrt *.c
blk-wbt.c: * window is then shrunk to 100 / sqrt(scaling step + 1).
blk-wbt.c: int_sqrt((rqd->scale_step + 1) << 8));
```
在 Linux 核心中， block 這個目錄下的 `blk-wbt.c` 有使用到整數平方根的操作，使用的函式為 `rwb_arm_timer` 。
- [ ] `blk-wbt.c`
```c=397
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));
} else {
```
閱讀實作 `writeback throttling` 這個機制的作者 Jens Axboe 的 [commit messages](https://git.kernel.dk/cgit/linux-block/commit/?h=wb-buf-throttle) ， throttling 這個機制最主要是讓背景執行的 buffered writeback 不會影響到其他正在執行的程式，作者發現在執行大量的 buffered writeback 時，他想要開啟 Chrome 會打不開。
> But for as long as I can remember, heavy buffered writers have not behaved like that. For instance, if I do something like this:
>
> $ dd if=/dev/zero of=foo bs=1M count=10k
>
> on my laptop, and then try and start chrome, it basically won't start
before the buffered writeback is done. Or, for server oriented
workloads, where installation of a big RPM (or similar) adversely
impacts database reads or sync writes. When that happens, I get people
yelling at me.
這個機制有點像是 [CoDel](https://en.wikipedia.org/wiki/CoDel) networking scheduling algorithm ， `blk-wb` 監控一段時間中的最小延遲，如果這段時間中，有最小延遲超過設定的目標，就增加 `scaling step` 與壓縮佇列的深度，接著將監控的時間縮小 `100 / sqrt(scaling step + 1)` 。其中的整數平方根運算對應的程式碼在 `/lib/math/int_sqrt.c` ，其實作的演算法也是使用 digit-by-digit calculation 。
- [ ] `int_sqrt`
```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;
}
```
#### 閱讀〈Analysis of the lookup table size for square-rooting〉
#### 閱讀〈Timing Square Root〉
Sign in with Wallet
Connect another wallet