How to use digital-by-digital implement sqrt() in C

二進位逼近

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

先看以下範例，我們利用此法求出

$\sqrt{13}$ 的整數部分

1. $13=2^3+2^2+2^0$, 因此我們依以下假設使用
   $(S_m{)}^2$ 來逼近
   
   $$13={(a_3{2}^3+a_2{2}^2+a_1{2}^1+a_0{2}^0)}^2={\left(\sum _{i=0}^3{a}_i{2}^i\right)}^2\text{和}$$

$$S_m=(a_3{2}^3+\cdots +a_m{2}^m)=\sum _{k=m}^3{a}_k{2}^k$$

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

• $(S_3{)}^2=(a_3{2}^3{)}^2,(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,(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,(2^1{)}^2\le 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,(2^1+2^0{)}^2\le 13,\text{ therefore }{a}_0=1$

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

換個角度想

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

我們都知道簡單的平方和公式：

${(a+b)}^2=a^2+2ab+b^2$

若套用在二進位並多加一位：

$\text{let }n1>n2>n3\ge 0$

$\begin{array}{rl}(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{array}$

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

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;
}

$N$的整數部分一定能使用二進位表示，因此假設

$N$的整數部分為n，且：

$n=\sum _{i=1}^k{2}^{n_i},\text{ for }{n}_i>n_{i+1}\mathrm{\forall }1\le i<k$

$N$ (i.e.x)

$\approx \sum _{i=1}^k{2}^{n_i}=(2^{n_1}+2^{n_2}+\cdots +2^{n_k}{)}^2$

看完後可以推測出m為最高位的平方，y"可能"為中間項，會說可能就是因為y長得非常神奇，無法看出是如何得到中間項。可以想像從第一項開始，若y為

$2(2^{n_1})(2^{n_2}+\cdots +2^{n_k})$應該會是個遠大於此算出的值。

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

但並不是

$2(2^{n_i})(2^{n_2}+\cdots +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}$

• 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}+\cdots +2^{n_k}$

More details

同樣假設

$\sqrt{N}$ 的整數部分為

$n=\sum _{i=1}^k{2}^{n_i}$
且

$N\approx (\sum _{i=1}^k{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}+\cdots +2^{n_k}{)}^2=(2^{n_1}{)}^2+2(2^{n_1})(2^{n_2}+\cdots +2^{n_k})+(2^{n_2}+\cdots +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}+\cdots +2^{n_k})$ 來保證

$N$ 在扣除後剩餘的數為

$2^{n_2}+2^{n_3}+\cdots +2^{n_k}$，那會不會發生存在某數

$t,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 在右移成相對應的數時，過程中仍是很大的數，所以過程中不會出現此情形進入條件式，此部分就不再詳加證明。

這就是個屌方法。