# 2024q1 Homework4 (quiz3+4)
contributed by < `w96086123` >

## [第三週測驗](https://hackmd.io/@sysprog/linux2024-quiz3)
### 測驗一

#### 公式解釋

為求 $N^2$ 的開根號，可利用 [Digit-by-digit calculation ](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) 的方式取得。

可將 $N^2$ 拆分為由 n 個整數之和： $N^2 = (a_n + a_{n-1} + a_{n-2} + ... + a_0)^2,a_m=2^m\ or\ a_m=0$
將此展開之後可得：\begin{split} \\
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(\displaystyle\sum_{i=1}^{n}a_i)+a_0]a_0 \\
=&\ 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_m =&\ a_n+a_{n-1}+...+a_m \\
=&\ p_{m+1}+a_m\\
\end{split}
因此可知 $P_0$ 即為所求。

目的是從試著找出所有 $a_m$ ，因此由 n 試到 0 並且因為 $a_m$ 只有兩種可能 $a_m=2^m or\ 0$ ，因此求 $a_m$ 時只須試驗 $P_m^2 \leq N^2$ ，即可知道 $a_m$ 的值。
\begin{cases}
P_m = P_{m+1} + 2^m, & \text{if $P_m^2 \leq N^2$}\\
P_m = P_{m+1}, & \text{otherwise}
\end{cases}

由於每輪計算 $N^2 - P_m^2$ 的成本過高，因此改為利用上輪的差值 $X_{m+1}$ 減去 $Y_m$ 可得 $X_{m}$
$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$
也就是紀錄上一輪的 $P_{m+1}$ 來調整。
\begin{split}
Y_m=&
\begin{cases}
c_m+d_m & \text{if } a_m=2^m \\
0 & \text{if } a_m=0
\end{cases}\\
c_m =& P_{m+1}2^{m+1} \\
d_m =& (2^m)^2 \\
\end{split}
拆分 ${Y_m}$ 可得上述的 $c_m$ 和 $d_m$ ，並且藉由位元運算推出下一輪
\begin{split}
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}=&\ (2^{m-1})^2 = (\dfrac{2^m}{2})^2 = \dfrac{(2^m)^2}{4}= \dfrac{d_m}{4}
\end{split}

由上述可知 `AAAA` 為 2 ，`BBBB` 為 1 。

#### 利用第二週測驗的 ffs 替代 `__builtin_clz`

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 >>= 1) {
        int b = z + m;
        z >>= 2;
        if (x >= b)
            x -= b, z += m;               
    }
    return z;
}

#### 在 Linux 核心原始程式碼找出對整數進行平方根運算的程式碼，並解說其應用案例
此函數在 `\lib\math\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;
}

找到 `int_sqrt` 被用在 `rwb_arm_timer` 中，而 `rwb_arm_timer` 的路徑是 `\block\blk-wbt.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));
	} else {
		/*
		 * For step < 0, we don't want to increase/decrease the
		 * window size.
		 */
		rwb->cur_win_nsec = rwb->win_nsec;
	}

	blk_stat_activate_nsecs(rwb->cb, rwb->cur_win_nsec);
}

此函數 `rwb` 代表的是 `request write back` ，因此可以知道它是一個寫入的函數。
此函數主要想要判斷是否有需要提昇寫入速度，若有需要則使用調整 window size 的方法進行提昇，且調整方法為透過 [Fast inverse square root](https://en.wikipedia.org/wiki/Fast_inverse_square_root) 。
但為何要使用 8 與 4 這樣狀況無法得知，只能猜測可能是基於實驗過後的結果所得出的結論。

### 測驗二

#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 >> CCCC);
    *mod = in - ((q & ~0x7) + (*div << DDDD));   
}

主要想法為想要使用 bitwise 的方式進行除法運算，但 10 不是 2 的冪次，因此無法直接使用，所以此方法想把 in 改變成 $\frac{8}{10} in$ ，後使用右移三位取得 $\frac{1}{10}in$ 。
但 in 無法使用直接轉成 $\frac{8}{10}$ ，因此使用逼近法取得數值 `q = (q >> 8) + x` 。
根據上述的計算過程，可得知要取得商數只需要將 `q` 右移三位即可，因此 `CCCC` 為 3 。
則想要取得餘數可利用餘數定理取得，`((q & ~0x7) + (*div << DDDD))` 則為 $g*Q$ ，因此 $mod = in - g * Q$ 即為餘數。

### 測驗三
主要作法與 `__builtin_clz` 相似，尋找 MSB 的數值以取得以 2 為底的對數整數。
以長度為 32 位元的為例，若 `i` 大於 $2^{16}$ 代表前 16 位元中含有 bit 為 1 ，以此類推可得以下程式：

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;
}

由於此方法可以直接使用 `__builtin_clz` 代替，因此也可得以下程式：

int ilog32(uint32_t v)
{
    return (31 - __builtin_clz(v));
}

### 測驗四
### 測驗五

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;       
}

`shift` 用來記錄當下 msb 的值，並且使用 `r` 進行 `shift` 的累加，以取得 `ilog2` 的數值。在這過程中，需要注意到我們要取得的是上限值，因此在一開始進行 `-1` ，主要為了處理 `x` 剛好為 2 的指數次方狀況。如果沒有做此動作，在輸出時會因為回傳的 `+1` 而比正確答案多一。

#### 改進程式碼，使其得以處理 x = 0 的狀況，並仍是 branchless

可以在進行運算前，先進行判斷是否 x 為 0 ，若是則改為 1 。後續處理方式一樣，程式碼如下。

int ceil_ilog2(uint32_t x)
{
    uint32_t r, shift;

    x |= (x == 0);

    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;
}