# [2024q1 Homework3 (quiz3+4)](https://hackmd.io/@sysprog/HkatSCZCT)
contributed by < [devarajabc](https://github.com/devarajabc) >

## 第三週測驗 `1-1`
### 1. 解釋上述程式碼運作原理，並嘗試用第 2 週測驗題提到的 ffs / fls 取代 `__builtin_clz`，使程式不依賴 GNU extension。

>參考 [wikipedia](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) 、[quiz3](https://hackmd.io/@sysprog/linux2024-quiz3) 、[浮點數運算和定點數操作](https://hackmd.io/@NwaynhhKTK6joWHmoUxc7Q/H1SVhbETQ?type=view)

#### 解釋程式碼運作原理
在解釋程式碼之前我們要先思考如何計算 $x$ 的平方根?
可以分成以下幾個步驟：
##### 1.
假設 $x$ 的平方根是 $N$ ，則 $x=N^2$
而其又可再拆解成：

$N^2=(a_{n}+a_{n-1}+a_{n-2}+...+{a_0})^2, a_m=2^m\ or\ a_m=0$

利用乘法公式可將 $(a_{n}+a_{n-1}+a_{n-2}+...+{a_0})^2$ 拆解成：
$\begin{split} \\
& 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 \\
\end{split}$

##### 2.
為了方便計算，使用 $P_m$ 來表示 $(\displaystyle\sum_{i=m}^{n}a_i)$ ，也就是 $P_m=a_n+a_{n-1}+...+a_m$ ，不難發現 $P_0=a_n+a_{n-1}+...+a_0$ 即為平方根 $N$

至於要如何得到 $P_0$ 呢？ 可以將 $P_m$ 寫成以下遞迴式：
$\begin{split} \\
P_m = P_{m+1} + a_m\
\end{split}$

接著利用 $P_m^2 \leq N^2$ 的特性來決定 $a_m=2^m$ 還是 $a_m=0$

將遞迴式改寫成：

\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}
遞迴的初始條件如下：

$P_n=a_n= 2^n$ ， $n$ 為能使 $a_n^2=(2^n)^2＝ 4^n \leq N^2$ 的最大值。

##### 3.
若在每次遞迴的過程中都要計算 $N^2 - P_m^2$ 的值會消耗大量的計算資源，因改以變數 $X_m$ 來記錄$N^2 - P_m^2$ 的值，並以遞迴的方式計算差值 $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$
為了方便計算，將 $Y_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 = P_{m+1}2^{m+1}$
* $d_m = (2^m)^2$

而 $c_{m},d_{m}$ 可以透過下列遞迴式求得：
* $c_{m-1}=(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}$

##### 4.
這個時候 **有 趣 的 事 情** 發生了，當 $m=-1$ 時：

$c_{-1} = P_02^0=P_0=N$

竟然透過計算 $c_m$ 就可以得到所求的平方根 $N$ ，連 $Y_m, X_m$ 之類的都不用管了嗎？ 當然不是，遞迴式需要透過 $X_m$ 的值來判斷 $a_m = 0$ 或是 $a_m = 2^m$ ，而 $X_m$ 的值需要透過 $Y_m$ 來計算。

##### 5.
接著就是如何透過程式碼執行以下遞迴式：
* $c_{m-1}=\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}$

初始條件改為：

$\begin{cases}
c_n=0 \\
d_n=2^n
\end{cases}$

$n$ 為能使 $a_n^2=(2^n)^2= d_n^2 = 4^n \leq N^2$ 的最大值，
此時 $c_n$ 為 $0$ ， 為什麼呢？
由於 $P_n=a_n= 2^n$ 且 $P_m = P_{m+1} + 2^m$
故 $P_{n+1}=0$, $c_n = 0$

##### 6.
程式運作的流程如下：

1. 利用變數 `x` 來儲存並計算 $X_{m+1} = N^2 - P_{m+1}^2 = N^2$
2. 利用變數 `m = 1UL << ((31 - __builtin_clz(x)) & ~1UL)` 來記錄 $d_n=2^n$
3. 利用 $d_{m-1}=\dfrac{d_m}{4}$ 來更新 `m`
4. 利用變數 `b` 來計算 $Y_m$ ，也就是 $c_m+d_m$
5. 計算 $c_{m-1}=\dfrac{c_m}{2}$
6. 如果 $X_m$ $>$ $Y_m$ ，代表 $a_m$ 不為零
* 利用 $X_m = X_{m+1} - Y_m$ 更新 $X_m$ 的值
* 將 $c_{m-1}$ 再加上 $d_m$

int i_sqrt(int x)
{
    if (x <= 1) /* Assume x is always positive */
        return x;

    int z = 0;//c_n
    for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
        int b = z + m;//b = y_m, m = d_m
        z >>= 1;
        if (x >= b)// x = x_n
            x -= b, z += m;               
    }
    return z;
}

sqrt_(26) = 5 // 正確答案應為 5.09901951359.....

測試完之後發現問題很大，在計算非完全平方數的平方根時會出現錯誤答案，
錯誤的原因不難理解，因為在程式的運算過程皆以整數 (int) 形態來儲存變數，為了讓輸出更接近真實答案，我打算使用[定點數](https://en.wikipedia.org/wiki/Fixed-point_arithmetic)的技巧來增加精準度。

#### 修改

#### 試用 `ffs / fls` 取代 `__builtin_clz`
> 參考 [ffs](https://man7.org/linux/man-pages/man3/ffs.3.html) :
> The ffs() function returns the position of the first (least
significant) bit set in the word i. The least significant bit is
position 1 and the most significant position is, for example, 32
or 64.

### 2. 在 Linux 核心原始程式碼找出對整數進行平方根運算的程式碼，並解說其應用案例，至少該包含 block 目錄的程式碼。
### 3. 閱讀〈Analysis of the lookup table size for square-rooting〉和 Timing Square Root，嘗試撰寫更快的整數開平分根運算程式。

[quiz-4](https://hackmd.io/@sysprog/linux2024-quiz4)