2024q1 Homework3 (quiz3+4)

contributed by < devarajabc >

第三週測驗 1-1

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

參考 wikipedia 、quiz3 、浮點數運算和定點數操作

解釋程式碼運作原理

在解釋程式碼之前我們要先思考如何計算 $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{array}{r}\\ & 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{array}$

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{array}{r}\\ P_m=P_{m+1}+a_m\end{array}$

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

將遞迴式改寫成：

$$\{\begin{array}{ll}{P}_m=P_{m+1}+2^m,& \text{if }{P}_m^2\le N^2\\ P_m=P_{m+1},& \text{otherwise}\end{array}$$
遞迴的初始條件如下：

$P_n=a_n=2^n$ ， $n$ 為能使 $a_n^2=(2^n{)}^2＝4^n\le 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}\ast a_m+a_m^2$
為了方便計算，將 $Y_m$ 分別以 $c_m=2P_{m+1}\ast a_m$ 和 $d_m=a_m^2$ 來表示：

$Y_m=\{\begin{array}{ll}{c}_m+d_m& \text{if }{a}_m=2^m\\ 0& \text{if }{a}_m=0\end{array}$

將 $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_m{2}^m=\{\begin{array}{ll}{c}_m/2+d_m& \text{if }{a}_m=2^m\\ c_m/2& \text{if }{a}_m=0\end{array}$
• $d_{m-1}={\displaystyle \frac{d_m}{4}}$

4.

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

$c_{-1}=P_0{2}^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{array}{ll}{c}_m/2+d_m& \text{if }{a}_m=2^m\\ c_m/2& \text{if }{a}_m=0\end{array}$
• $d_{m-1}={\displaystyle \frac{d_m}{4}}$

初始條件改為：

$\{\begin{array}{l}{c}_n=0\\ d_n=2^n\end{array}$

$n$ 為能使 $a_n^2=(2^n{)}^2=d_n^2=4^n\le 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}={\displaystyle \frac{d_m}{4}}$ 來更新 m
4. 利用變數 b 來計算 $Y_m$ ，也就是 $c_m+d_m$
5. 計算 $c_{m-1}={\displaystyle \frac{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) 形態來儲存變數，為了讓輸出更接近真實答案，我打算使用定點數的技巧來增加精準度。

修改

試用 ffs / fls 取代 __builtin_clz

參考 ffs :
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