# How to use [digital-by-digital](https://hackmd.io/@sysprog/sqrt#Digit-by-digit-Method) implement [sqrt()](https://hackmd.io/@sysprog/linux2025-quiz2#測驗-2) in C

## 二進位逼近

在實作前先介紹何謂二進位逼近。

先看以下範例，我們利用此法求出 $\sqrt{13}$ 的整數部分

1. $13 = 2^3 + 2^2 + 2^0$, 因此我們依以下假設使用 $(S_{m})^2$ 來逼近
$$13 = \left( a_{3}2^3 + a_{2}2^2 + a_{1}2^1 + a_{0}2^0 \right)^2 = \left(\sum^{3}_{i=0} a_{i}2^{i} \right)^2 \text{和}$$

$$S_{m} = (a_{3}2^3 + \dots + a_{m}2^m) = \sum^{3}_{k=m} a_{k}2^k$$

2. 由最大項 $m=3$ 開始逼近

- $(S_{3})^2 = (a_{3}2^3)^2, \quad (2^3)^2 > 13, \text{ therefore }a_{3}=0$

- $(S_{2})^2 = (a_{3}2^3 + a_{2}2^2)^2=(a_{2}2^2)^2, \quad (2^2)^2 > 13, \text{ therefore } a_{2}=0$

- $(S_{1})^2 = (a_{3}2^3 + a_{2}2^2 + a_{1}1^1)^2=(a_{1}2^1)^2, \quad (2^1)^2 \leq 13, \text{ therefore } a_{1}=1$

- $(S_{0})^2 = (a_{3}2^3 + a_{2}2^2 + a_{1}1^1 + a_{0}2^0)^2=(2^1 + a_{0}2^0)^2, \quad (2^1+2^0)^2 \leq 13, \text{ therefore } a_{0}=1$

3. 所以 $\sqrt{13}$ 的整數部分為 $(2^1 + 2^0) = 3$

## 換個角度想

以上方法是從最高位元逐步向最低位元**遞增**逼近目標數值。
同樣地，我們也可以透過扣除的方式來逼近。接下來將說明如何使用此方法。

我們都知道簡單的平方和公式：
$\left(a+b \right)^2 = a^2+2ab+b^2$

若套用在二進位並多加一位：
> $\text{let } n1 > n2 > n3 \geq 0$

$\begin{align}
(2^{n1}+2^{n2}+2^{n3})^2
&= (2^{n1}+(2^{n2}+2^{n3}))^2 \\
&= (2^{n1})^2 + 2(2^{n1})(2^{n2}+2^{n3}) + (2^{n2}+2^{n3})^2
\end{align}$

因此可以透過類似遞迴逐漸減去最高位及中間項，再重複利用同樣的方法得到下一位，這便是程式實作的核心想法。

## Impliment in C

程式實作如下

void sqrti(uint64_t x)
{
    uint64_t m, y = 0;

    int total_bits = 64;

    /* clz64(x) returns the count of leading zeros in x.
     * (total_bits - 1 - clz64(x)) gives the index of the highest set bit.
     * Rounding that index down to an even number ensures our starting m is a
     * power of 4.
     */
    int shift = (total_bits - 1 - clz64(x)) & ~1;
    m = 1ULL << shift;
        while (m) {
            uint64_t b = y + m;
            y >>= 1;
            if (x >= b) {
                x -= b;
                y += m;
            }
            m >>= 2;
        }
    // return y;
    return y;
}

:::info
$N$的整數部分一定能使用二進位表示，因此假設$N$的整數部分為n，且：
$n = \sum^{k}_{i=1}2^{n_i}, \text{ for }n_i>n_{i+1} \quad \forall 1 \leq i < k$
$N$ (i.e.`x`) $\approx \sum^{k}_{i=1}2^{n_i} = ( 2^{n_1} + 2^{n_2} + \dots + 2^{n_k})^2$
:::

看完後可以推測出`m`為最高位的平方，`y`"可能"為中間項，會說可能就是因為`y`長得非常神奇，無法看出是如何得到中間項。可以想像從第一項開始，若`y`為 $2(2^{n_1})(2^{n_2} + \dots + 2^{n_k})$應該會是個遠大於此算出的值。

然而，他確實是平方和中間項的值。

但並不是 $2(2^{n_i})(2^{n_2} + \dots + 2^{n_k})$ ，讓我們觀察每次扣除時`y`的變化：
> m 的部分明顯為 $(2^{n_i})^2$，因此只針對 y 來做分析

1. `m` = $(2^{n_1})^2$ 第一次進入`if (x >= b)`，扣除的 y 為 0
2. `m` = $(2^{n_2})^2$ 第二次進入`if (x >= b)`，扣除的 y 為
$$ y = 2^{2n_1 - (n_1-n_2-1)}=2^{n_1+n_2+1} = 2(2^{n_1})(2^{n_2})$$
並且 y 更新為 $2^{n_1+n_2}+2^{2n_2}$

:::warning
- m 每次右移 2 bit，所以 m 從 $2^{2n_1}$ 右移成 $2^{2n_2}$ 時，迴圈總共跑了 $n_1 - n_2$ 次
- 並且，因為是更新 b 才右移 y，所以 b 中的 y 會比同個迴圈的 y 少更新一次，所以扣除的 y 從 $2n_1$ 右移了 $n_1 - n_2 -1$ 次
- 後續也是相同的概念
:::

3. `m` = $(2^{n_3})^2$ 第二次進入`if (x >= b)`，扣除的 y 為
$$y = \frac{2^{n_1+n_2}+2^{2n_2}}{2^{n_2-n_3-1}} = 2^{n_1+n_3+1}+2^{n_2+n_3+1} = 2(2^{n_3})(2^{n_1}+2^{n_2})$$
並且 y 更新為 $2^{n_1+n_3}+2^{n_2+n_3}+2^{2n_3}$

最後就能得到 $y = 2^{n_1} + 2^{n_2} + \dots + 2^{n_k}$

## More details

> 同樣假設 $\sqrt{N}$ 的整數部分為 $n = \sum^{k}_{i=1}2^{n_i}$
且 $N \approx (\sum^{k}_{i=1}2^{n_i})^2$

在此實作中，找到 $2^{n_i}$ 後的算法為：

- 扣除 $(2^{n_i})^2$
- 在找到 $2^{n_{i+1}}$ 時扣除 $2(2^{n_i})(2^{n_{i+1}})$

也就是說這不是根據

$$( 2^{n_1} + 2^{n_2} + \dots + 2^{n_k})^2 = (2^{n_1})^2 + 2(2^{n_1})(2^{n_2} + \dots + 2^{n_k}) +( 2^{n_2} + \dots + 2^{n_k})^2$$

來實現，而是用了類似像累積的方式，若將每次扣除的值條列出來如下：

| m | b = y + m | y (after shifted) |
| -------- | -------- | -------- |
| $(2^{n_1})^2$ | $0+(2^{n_1})^2$ | 0 |
| $(2^{n_2})^2$ | $2(2^{n_1})(2^{n_2})+(2^{n_2})^2$ | $2^{n_1+n_2}$ |
| $(2^{n_3})^2$ | $2(2^{n_1}+2^{n_2})(2^{n_3})+(2^{n_3})^2$ | $2^{n_1+n_3}+2^{n_2+n_3}$ |
| $(2^{n_4})^2$ | $2(2^{n_1}+2^{n_2}+2^{n_3})(2^{n_4})+(2^{n_4})^2$ | $2^{n_1+n_4} + 2^{n_2+n_4} + 2^{n_3+n_4}$ |

乍看之下，扣除的數就像在累積成下一個平方和一樣，只不過這是用減的
事實上也是如此，每階段總扣除如下：

| m | $\sum{b}$ |
| -------- | -------- |
| $(2^{n_1})^2$ | $(2^{n_1})^2$ |
| $(2^{n_2})^2$ | $(2^{n_1} + 2^{n_2})^2$ |
| $(2^{n_3})^2$ | $(2^{n_1} + 2^{n_2} + 2^{n_3})^2$ |
| $(2^{n_4})^2$ | $(2^{n_1} + 2^{n_2} + 2^{n_3}+ 2^{n_4})^2$ |

---

最後，腦中出現了個問題：

若在迴圈中每次扣除的數不是像 $(2^{n_1})^2 +$ ==$2(2^{n_1})(2^{n_2} + \dots + 2^{n_k})$== 來保證 $N$ 在扣除後剩餘的數為 $2^{n_2}+2^{n_3} + \dots + 2^{n_k}$，那會不會發生**存在某數 $t, \quad n_1 > t > n_2$ 且在`m`還未右移成 $(2^{n_2})^2$ 時進入 if 條件式？** 也就是說 $m$ 在 $(2^{n_1})^2$ 右移成 $(2^{n_2})^2$ 前，在 $m = (2^{t})^2$ 時，進入了 if 條件式？

不會。

因為 b = y + m，儘管過程中可能出現 x > m 的情形，但還有 y 這個變數在把關，y 在右移成相對應的數時，過程中仍是很大的數，所以過程中不會出現此情形進入條件式，此部分就不再詳加證明。

這就是個屌方法。