# 2024q1 Homework4 (quiz3+4)
contributed by < [youjiaw](https://github.com/youjiaw/lab0-c) >

## 第 3 週測驗題
### 測驗 `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_0})^2, a_m=2^m\ or\ a_m=0
$$
並根據 $(x+y）^{2}=x^{2}+2xy+y^{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}^n a_{i})+a_{0}]a_{0}
$$

以 n=2 舉例，$N^2=(a_{2}+a_{1}+a_{0})^2$ 可以用下方的矩陣表示
$$
\left[
\begin{array}{ccc}
a_0a_0 & a_1a_0 & a_2a_0 \\
a_0a_1 & a_1a_1 & a_2a_1 \\
a_0a_2 & a_1a_2 & a_2a_2 \\
\end{array}
\right]
$$
將算式調整後可得
$$
N^2=a_{2}^2+[2a_{2}+a_{1}]a_{1}+[2(a_{2}+a_{1})+a_{0}]a_{0}
$$
分別代表以下三個部分
$$
\left[
\begin{array}{ccc}
& & \\
& &\vdots \\
& ... & a_2a_2 \\
\end{array}
\right]
$$

$$
\left[
\begin{array}{ccc}
\ddots & \vdots & \vdots \\
... & a_1a_1 & a_2a_1 \\
... & a_1a_2 & \ddots \\
\end{array}
\right]
$$

$$
\left[
\begin{array}{ccc}
a_0a_0 & a_1a_0 & a_2a_0 \\
a_0a_1 & \ddots & \vdots \\
a_0a_2 & \vdots & \\
\end{array}
\right]
$$

接著令 $P_{m}=\sum_{i=m}^n a_{i}$，則 $P_{0}=a_{n}+a_{n-1}+...+a_{0}$ 就會是我們要求的平方根 N，因此我們需要找出每個 $a_{m}$ 的值。

想檢測 $a_{m}$ 的值是 $2^m$ 還是 0，可以計算 $P_{m}^2$ 與 $N^2$ 的差值，若 $P_{m}^2 \leq N^2$ 就代表 $a_{m}$ 為 $2^m$，否則為0，但是此方法每輪的運算成本太高了。

若是令
* $X_{m}=N^2-P_{m}^2=X_{m+1}-Y_{m}$
* $Y_{m}=P_{m}^2-P_{m+1}^2=(P_{m+1}+a_{m})^2-P_{m+1}^2=2P_{m+1}a_{m}+a_{m}^2$

就可以透過紀錄 $P_{m+1}$ 來取代 $P_{m}^2$ 的計算。

最後將 $Y_{m}$ 拆為 $c_{m}$、$d_{m}$，可得
* $c_{m}=2P_{m+1}a_{m}$
* $d_{m}=a_{m}^2$
* $Y_m=\left.
\begin{cases}
c_m+d_m & \text{if } a_m=2^m \\
0 & \text{if } a_m=0
\end{cases}
\right.$

並藉由位元運算，由 $c_{m}$、$d_{m}$ 推出下一輪的 $c_{m-1}$、$d_{m-1}$。
* $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}=\dfrac{d_m}{4}$

最後求得的 $c_{-1}$ 即為所求的 $P_{0}=N$。

程式碼如下：

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

在確認 `x` 的範圍後，會使用三個變數，`z`、`m` 與 `b`，分別對應到先前的 $c_m$、$d_m$ 與 $Y_{m}$。
在每輪迴圈中，無論 $a_{m}$ 值為何，$c_{m}$ 都會除以 2，$d_m$ 則是除以 4，因此 `AAAA` 為 2，`BBBB` 為 1。

但是我還沒理解為什麼要做 `& ~1UL`，相當於限制左移位元數為偶數，目前猜測是因為 $d_{m}=a_{m}^2=2^{2m}$，所以 `m` 僅可能為 2 的偶數次方。

#### 嘗試用 ffs / fls 取代 __builtin_clz
ffs 會回傳最低有效位元的索引，從 1 開始計算，若 `x` 為 0 則會回傳 0。因此可以不斷將 `x` 右移 `ffs(x)` 個位元，直到 `x` 為 0 時就停止，最後累計的位移次數減一就會是最高有效位元的索引。

因此加入以下程式碼

int tmp = x, msb = 0;
while (tmp) {
    int index = ffs(tmp);
    tmp >>= index;
    msb += index;
}

並更改迴圈條件為

-    for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
+    for (int m = 1UL << ((msb - 1) & ~1UL); m; m >>= 2) {

而 fls 是直接回傳最高有效位元的索引

，因此更改如下

-    for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
+    for (int m = 1UL << ((fls(x) - 1) & ~1UL); m; m >>= 2) {

### 測驗 `3`
此程式碼目的是取得輸入值 `i` 在二進位表示中最高有效位元的索引。

版本一是將 `i` （$2^n$）不斷右移一個位元直到為 0，此時右移次數即為 n，迴圈執行次數也為 n。
版本二則是依照 `i` 的大小，提供四種右移的方式，只要 `i` 大於 $2^k$，就一次右移 k 個位元，讓迴圈的執行次數可以小於 n。

static size_t ilog2(size_t i)
{
    size_t result = 0;
    while (i >= 2^16) {
        result += 16;
        i >>= 16;
    }
    while (i >= 2^8) {
        result += 8;
        i >>= 8;
    }
    while (i >= 2^4) {
        result += 4;
        i >>= 4;
    }
    while (i >= 2) {
        result += 1;
        i >>= 1;
    }
    return result;
}

版本三使用了 `__builtin_clz` 改寫，但因為輸入值不能為 0，所以先與 1 做 `|` 運算。

return (31 - __builtin_clz(v | 1));

### 測驗 `5`
此程式碼以測驗三做延伸，計算$\lceil log_2(x) \rceil$。

<!-- 一開始的 `x-1` 是為了讓 -->

若 (x > 0xFFFF) 回傳 1，則 `r` 的值為 `10000`，因此會將 `x` 右移 16 個位元。

r = (x > 0xFFFF) << 4;                                   
x >>= r;

對應到測驗三的

while (i >= 2^16) {
    result += 16;
    i >>= 16;
}

而因為每次左移的位元數都不同，代表 `r` 和 `shift` 不會有相同位元都為 1 的情況，因此計算 result 的方式可以用 `r |= shift` 代替。

shift = (x > 0xFF) << 3;
x >>= shift;
r |= shift;
shift = (x > 0xF) << 2;
x >>= shift;
r |= shift;

最後再判斷 `x` 是否大於 1。

shift = (x > 0x3) << 1;
x >>= shift;
return (r | shift | x > 1) + 1;

#### 改進程式碼
最直覺的想法就是當 `x` 為 0 時減 0，其他情況減 1。目前採用下方的方法。

-    x--;
+    x -= (x > 0)

## 第 4 週測驗題
### 測驗 `1`