# 2022q1 Homework3 (quiz3)
contribute by < `kdnvt` >

### 測驗 `11`

unsigned long i_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;
}

[直式開方法](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation)

雖然教授問答的時候是從二進位的角度去想，但我第一次看到類似的原理其實是從十進位的直式開方法。不過直式開方法所有進位都適用，所以從十進位的想法轉成二進位也沒有很難。
直式開方法是把要開方的數拆成二的冪相加（如果是十進位就是十的冪），就可以拆成以下的形式：

$x = (000b_0b_1b_2\cdots b_{n-1}b_{n})^2$

上面的 $000b_0b_1b_2\cdots b_{n-1}b_{n}$ 為 bits pattern 。
$b_0$ 為最高位的 1 ，
寫成數學式如下：

$x = (2^{n}b_0+2^{n-1}b_1+\cdots+2^1b_{n-1}+2^0b_{n})^2$

$a_i = 2^{n-i}b_i$

$\begin{aligned}x
&=(a_0+a_1+\cdots+a_{n-1}+a_{n})^2 \\
&= a_0^{2}+2a_0a_1 + a_1^2 + 2(a_0+a_1)a_2 + a_2^2 + \cdots + a_n-1^2 + 2\left(\sum\limits_{i = 0}^{n-1}{a_i}\right)a_n + a_n^2 \\
&=a_0^2 + (2a_0+a_1)a_1 + [2(a_0+a_1)+a_2]a_2 + \cdots + [2\left(\sum\limits_{i = 0}^{n-1}{a_i}\right)+a_n]a_n \end{aligned}$

最後的形式為 $[2\left(\sum\limits_{i = 0}^{n-1}{a_i}\right)+a_n]a_n$ 的相加，而這每一項其實就是程式中每一輪迴圈

    if (x >= b) {
        x -= b;
        y += m;
    }

要判斷並可能減掉的 `b` 。

以下大致講解原理：
迴圈從 0 開始。
$y_i$ 為每個第 $i$ 個迴圈一開始的 `y` 值。
$m_i$ 為每個第 $i$ 個迴圈一開始的 `m` 值。
$l_i$ 為 $\sqrt m_i$ 。

$y_0 = 0$

$m_i = l_i^2$

$l_i = 2^{n-i}$

因為確定 $a_0$ 是最高位不為零，所以 $a_0 = l_0$

$y_1 = m_0 = l_0^2 = a_0^2 = l_0a_0$

而每經過一輪 $m$ 會右移兩位，他的平方根 $l$ 就等於右移一位：
$m_i = m_{i+1}\times4$
$l_i = l_{i+1}\times2$

$y_1 = l_0a_0$
$y_1 = 2l_1a_0$
$y_1+m_1 = 2l_1a_0 + l_1^2 = (2a_0+l_1)l_1$
如果此時 $y_1+m_1$ 的結果小於當前的 $x$ ，代表 $x$ 的展開式中 $(2a_0+a_1)a_1$ 這項不為零，也就是 $a_1 = 2^{n-1} = l_1$ 。將 $y_1$ 右移一位後加上 $m_1$ ，就會變成
$y_2 = \dfrac{y_1}{2} + m_1 = l_1a_0 + l_1^2 = l_1(a_0+a_1) = 2\times l_2(a_0+a_1)$

而如果此時 $y_1+m_1$ 的結果大於當前的 x ，代表 $x$ 的展開式中 $(2a_0+a_1)a_1$ 這項等於零，也就是 $a_1 = 0$ 。將 $y_1$ 右移一位後就會變成
$y_2 = \dfrac{y_1}{2} = l_1a_0 = l_1(a_0+a_1) = 2\times l_2(a_0+a_1)$

照同樣的原理推下去可以發現， $y_j$ 其實就等於 $2\left(\sum\limits_{i = 0}^{j-1}{a_i}\right)l_j$ 。因為不是所有的 $l_j 都會等於 a_j$ ，所以我用另一個變數 $z_j = 2\left(\sum\limits_{i = 0}^{j-1}{a_i}\right)a_j$ 來表達原式中的部份，$z_j = \dfrac{y_ja_j}{l_j}$ ，如果 $a_j = l_j$ 也就是 $a_j \neq 0$ ，則 $z_j$ 會等於 $y_j$ 而且下圖的 $a_i^2$ 會等於 $m_i$ ：

$\underbrace{a_0^2}_{z_0+a_0^2} + \underbrace{(2a_0+a_1)a_1}_{z_1+a_1^2} + \underbrace{[2(a_0+a_1)+a_2]a_2}_{z_2+a_2^2} + \cdots + \underbrace{[2\left(\sum\limits_{i = 0}^{n-1}{a_i}\right)+a_n]a_n}_{z_n+a_n^2}$

$z_i+a_i^2 = \left\{
\begin{array}{r}
0 , & \text{if }a_i=0 \\
y_i + m_i , & \text{if } a_i=l_i
\end{array}
\right.$

每一輪迴圈就是判斷 $x$ 有沒有 $z_i+a_i^2$ 並更新 $y_i$，

到了最後的 $y_n = 2\left(\sum\limits_{i = 0}^{n-1}{a_i}\right)l_n = 2\left(\sum\limits_{i = 0}^{n-1}{a_i}\right)2^0 = 2\left(\sum\limits_{i = 0}^{n-1}{a_i}\right)$ ， $y_n$ 會先往右移變成 $\sum\limits_{i = 0}^{n-1}{a_i}$ ，在加上最後的 $a_n$ ，變成 $\sum\limits_{i = 0}^{n}{a_i}$ ，就是我們要求的平方根。