2024q1 Homework4 (quiz3+4)

# 2024q1 Homework4 (quiz3+4) contributed by < [`jouae`](https://github.com/jouae) > Github [`jouae/sqrt`](https://github.com/jouae/sqrt) 這個頁面談及到關於開根號的演算法跟分析，其中涵蓋 Linux int_sqrt.c Python sqrt 等整數開根號部分 This page discusses algorithms and analysis related to the numerical method for integer square root, including integer square root operations in Linux int_sqrt.c and Python's sqrt. 關鍵字：根號、平方根、快速平方根、牛頓法、查表平方根、整數平方根 keywords: integer square root, Newton's method, fast integer square root, Heron's method ## 平方根 ### Digit-by-digit Method 在例題中使用平方和公式來理解，如何從乘法到單純加減與位移運算來開平方根。先觀察其規律 $$ (a+b)^2 = a^2+(2a+b)b $$ 三元平方和 $\begin{align} (a+b+c)^2 &= ((a+b)+c)^2 \\ &= (a+b)^2 + (a+b)c +c(a+b) + c^2 \\ &= (a+b)^2 + (2(a+b)+c)c \\ &= a^2 + (2a+b)b + (2(a+b)+c)c \end{align}$ 四元平方和 $\begin{align} (a+b+c+d)^2 &= ((a+b+c)+d)^2 \\ &= (a+b+c)^2 + (a+b+c)d + d(a+b+c) + d^2 \\ &= (a+b+c)^2 + (2(a+b+c)+d)d \\ &= a^2 + (2a+b)b + (2(a+b)+c)c + (2(a+b+c)+d)d \end{align}$ 考慮一般化的平方和 $$ N^2 = (a_n+a_{n-1}+\dots+a_{1}+a_0)^2 $$ 觀察規律可以發現 $m$ 元平方和可以藉由 $1$ 到 $m-1$ 元平方和總和，同時加上一項 $$ Y_{m}=\left[\left(2\sum^{n}_{i=m+1} a_{i}\right) + a_{m}\right]a_{m} = \left( 2P_{m+1}+ a_{m}\right)a_{m} $$ 其中 $P_{m+1}=a_n+a_{n-1}+\dots+a_{m+2}+a_{m+1}$ 而 $P_0=N$ 且 $$ P_{m} = P_{m+1} + a_{m},\quad \text{for } 1\leq m<n $$ 總結上述兩個觀察的結果： 1. $n$ 元平方和可以藉由 $1$ 到 $n-1$ 元平方和總和，再加一項 $Y_{n}$ 2. $Y_{n}=\left( 2P_{n-1}+ a_n\right)a_n$ 其中 $P_{n-1} = P_{n-2} + a_{n-1}$ 數值系統在二進位表示下為： $$ x = (00b_0b_2\dots b_{n-1}b_n)_2^2 $$ $b_0$ 為最高位為 $1$ 之位元，二進位轉十進位如下： $$ \sum_{i=0}^n b_i\times 2^{n-i} $$ 假設 $N^2$ 可以用二進位系統逼近： $$ \begin{align} N^2 &\approx (b_n\times 2^0 + b_{n-1}\times 2^{1}+\dots b_ 1\times2^{n-1} +b_0\times2^n)^2 \\ &= (a_n + a_{n-1} + \dots a_1 + a_0)^2 \end{align} $$ 其中 $b_i\in\lbrace 0,1\rbrace,\text{ for }i=0,1\dots,n$ 且 $P_m=P_{m+1}+a_{m}=P_{m+1}+b_{m}\times 2^{m}$。故逼近過程只要從 $m=n$ 也就是最高位元為 $1$ 開始計算到最低位元 $m=0$，確認 $a_m=2^m$ 或是 $a_m=0$ 對 $0 \leq m < n$ ，藉由上述 $P_m=P_{m+1}+b_{m}\times 2^{m}$ 可以確認有兩種可能： $$ P_m= \begin{cases} P_{m+1}+2^{m}, &\text{if }P_m^2 \leq N^2,\\ P_{m+1}, &\text{otherwise.} \end{cases} $$ 舉例來說，設 $N^2=10$ 則 $N^2=(10)_{10}=(1010)_2=(b_3\times2^3+b_1\times 2^1)$ ，最高位元 $b_3$ 從 $m=3$ 開始計算。 * $m=3$ 時 $P_{3}^2 = a_3^2 = (2^3)^2 = 64> 10 = N^2$ ，因此 $P_3=P_4$ 且 $a_3=0$。 * $m=2$ 時 $P_2^2 = (a_3 + a_2)^2 = (2^2)^2 = 16 > 10 = N^2$ ，因此 $P_2 = P_3$ 且 $a_2=0$ * $m=1$ 時 $P_1^2 = (a_3+a_2+a_1)^2 = (2^1)^2 = 4 < 10 =N^2$ ，因此 $P_1=P_2+2^1$ 且 $a_1=2^1$ * $m=0$ 時 $P_0^2 = (a_3+a_2+a_1+a_0)^2 = (2^1 + 2^0)^2 = 3^2= 9 \leq 10 = N^2$ ，因此 $P_0=P_1+ 2^0$ 且 $a_0=2^0$ 故 $N=P_0= a_1 + a_0 = 2^1+2^0= 3$ ，因此 $10$ 的整數平方根為 $3$。令 $X_m=N^2-Y_m$ 為差值表實際值與逼近值的誤差，其中 $Y_m=\left( 2P_{m+1}+ a_{m}\right)a_{m}$ 運用以下性質可以縮減上述計算過程： 1. $X_m=N^2-P_m^2=X_{m+1}-Y_m$ 2. $Y_m=P_m^2-P_{m+1}^2=(2P_{m+1}+a_m)a_m$ 由於 $Y_m=2a_mP_{m+1}+a_m^2$ 與 $a_m=2^m$ 則可以將 $Y_m$ 表示為 $$ Y_m = P_{m+1}2^{m+1} + (2^{m})^2 $$ 令 $c_m =2a_mP_{m+1}= P_{m+1}2^{m+1}$ 與 $d_m = a_m^2=(2^{m})^2$ 。重新將 $Y_m$ 表示為 $$ Y_m = \begin{cases} c_m+d_m, &\text{ if }a_m =2^m, \\ 0, &\text{ if } a_m = 0. \end{cases} $$ 同時藉由 $P_{m} = P_{m+1}+a_{m}$ 則 $$ \begin{align} c_m &= P_{m+1}2^{m+1} \\ &= (P_m-a_m)2^{m}2 \\ &= 2P_m2^{m} - 2a_m2^{m} \\ &= 2c_{m-1} - 2a_m2^m. \end{align} $$ 如果 $a_m\neq 0$ 則 $$ c_{m-1} = c_m/2 +d_m. $$ 同時 $$ d_m = (2^m)^2 = (2\times2^{m-1})^2 = 4(2^{m-1})^2 $$ 則 $d_{m-1} = d_m/4$。總結迭代關係為 $$ c_{m-1} = \begin{cases} c_m/2 +d_m, &\text{ if } a_m\neq 0, \\ c_m/2, &\text{ if } a_m = 0. \end{cases} $$ 與 $$ d_{m-1} = d_m/4 $$ 則可以藉由 $c_{-1}$ 得知平方根為 $$ c_{-1} = P_02^0 = P_0 = N. $$ ### Japanese Algorithm 英國 University of Sheffield 數學與統計系任教的 Frazer Jarvis 教授在他的[著作](http://www.afjarvis.org.uk/maths/index.html)中有一篇〈[Square roots by subtraction](http://www.afjarvis.org.uk/maths/jarvisspec02.pdf)〉介紹了一個計算平方根的方法。據其所言： > When I was at school, my mathematics teacher showed me the following very strange method to work out square roots, using only subtraction, which is apparently an old Japanese method. Mayer Goldberg 在〈[A Spigot-Algorithm for Square-Roots: Explained and Extended](https://arxiv.org/pdf/2312.15338.pdf)〉稱其為 Japanese algorithm ，並提出一改進算法，在內文 3.2 Reconstructing the Original “Japanese Algorithm” 中提及，所謂的 Japanese algorithm 似乎是針對「算盤」計算而演化來。 > The “Japanese algorithm” is essentially the optimized algorithm of the previous section, with one additional, non-mathematical optimization, that is probably intended to simplify computation on an abacus. :::success 要會算盤計算的人驗證。 ::: 以下敘述 Frazer Jarvis 的〈[Square roots by subtraction](http://www.afjarvis.org.uk/maths/jarvisspec02.pdf)〉一文中 Japanese algorithm 作法跟證明過程。此算法在十進位系統中，針對所求 $n$ 的平方根逼近先初始化數對 $(a,b)$ 為 $$ (a,b) = (5n, 5). $$ 再藉由兩個迭代規則更新數對 $(a,b)$ * 規則一(R1)：當 $a\geq b$ 時，令 $a=a-b$ 且 $b=b+10$ * 規則二(R2)：當 $a<b$ 時，在 $a$ 的右側插入兩個零且在 $b$ 最右側數來第二位插入一個零。虛擬碼如下： ``` /* Pesudo code of finding square root of n */ Initial a = 5*n and b = 5 as a pair (a, b) MAX_DIGITS - The maximum number of digits represented by a square root. do loop until MAX_DIGITS /* R1 */ if a >= b then a = a - b b = b + 10 /* R2 */ elif a < b then a = 100 * a b = 10 * b - 45 end do return b as the approximation of n ``` 顯而易見的這個算法的缺陷之一在於 `int a` 滿足 R2 時，假如 `a = 100*a` 超過 `INT_MAX` 會有 overflow 的情況。 ### Spigot-Algorithm Mayer Goldberg 基於以下目標改進 Japanese algorithm ： >1. To uncover the layered representation that is at the heart of this algorithm, in order to understand what insights it captured about square-roots. >2. To use this layered representation to understand how the algorithm works. >3. To generalize the algorithm to other roots. 原文中使用十進位系統，此處使用二進位系統表示 $b_i\in \lbrace 0,1 \rbrace$ 且專注於整數型態根號 $$ [b_0;b_1;\dots;b_{n-1};b_n] = 2^nb_0+2^{n-1}b_1 + \dots +2^1 b_{n-1} + 2^0b_n=D_n. $$ 此處 $b_0$ 為最高位為 1 的位元。同時藉由上述表示法有以下關係： $$ \begin{align} [b_0;b_1;\dots;b_{n};b_{n+1}] &= 2^{n+1}b_0+2^{n}b_1 + \dots +2^1 b_{n} + 2^0b_{n+1} \\ &= 2(2^nb_0+2^{n-1}b_1 + \dots +2^0 b_{n}) + b_{n+1} \\ &= 2[b_0;b_1;\dots;b_{n}] + b_{n+1} \\ &= 2D_n +b_{n+1}. \end{align} $$ 手法如 digit-by-digit 一樣，令所求平方根為 $$ M=(2^0b_n + 2^1b_{n-1} + \dots +2^{n-1}b_1 +2^nb_0)^2 $$ 藉由 Nested Parentheses 算平方和的方法，每次可以提取出一項 $A_{i+1}=(b_{i+1}^2+4D_ib_{i+1}+4A_{i})$ 且 $A_1 = b_1^2 + 4D_0b_1 + 4b_0^2$ 。舉例來說 $n=3$ 時，可以寫成 $b_3$ 與 $2^1b_2+2^2b_1+2^3b_0$ 相加的平方 $$ \begin{align} A_3=(2^0b_3 + 2^1b_2 + 2^2b_1 + 2^3b_0)^2 &= (b_3^2 + (2^1b_2 + 2^2b_1 + 2^3b_0))^2 \\ &= b_3^2 + 2(2^1b_2 + 2^2b_1 + 2^3b_0)b_3 + (2^1b_2 + 2^2b_1 + 2^3b_0)^2 \\ &= b_3^2 + 4(2^0b_2 + 2^1b_1 + 2^2b_0)b_3 + [2(2^0b_2 + 2^1b_1 + 2^3b_0)]^2 \\ &= b_3^2 + 4D_2b_3 + 4(2^0b_2 + 2^1b_1 + 2^2b_0)^2 \\ \end{align} $$ 經由這一輪的拆解後就有另一個平方和 $(2^0b_2 + 2^1b_1 + 2^2b_0)^2$ 包含在第一個括號裏頭，一樣用同樣的手法可以寫成 $$ \begin{align} A_2=(2^0b_2 + 2^1b_1 + 2^2b_0)^2 &= (b_2 + (2^1b_1 + 2^2b_0))^2 \\ &= b_2^2 + 2(2^1b_1 + 2^2b_0)b_2 + (2^1b_1 + 2^2b_0)^2 \\ &= b_2^2 + 4(2^0b_1 +2^1b_0)b_2 + [2(2^0b_1 + 2^1b_0)]^2 \\ &= b_2^2 + 4D_1b_2 + 4(2^0b_1 + 2^1b_0)^2 \\ &= b_2^2 + 4D_1b_2 + 4A_1 \end{align} $$ 再來一次 $$ A_1 = (2^0b_1 + 2^1b_0)^2= b_1^2 + 4D_0b_1 + 4b_0^2 $$ 在這個算法重要的地方是我們想要專注於使用減法來找出平方根，所以藉由改寫一些形式我們可以對 $b_i^2+4D_ib_{i+1}$ 以減法的方式找出 $b_{i+1}$ 原文中由於是十進位系統，所以當 $b_{i+1} \neq 0$ 時需要知道起碼 $b_{i+1} - 1$ 個可能才能推出 $b_{i+1}$ 。但在二進位系統中，只要確認 $b_{i+1}$ 是不是 $0$ 即可  那 $b_i^2+4D_ib_{i+1}$ 可以寫成 $$ \begin{align} b_{i+1}^2 + 4D_ib_{i+1} &= \begin{cases} 0, &\text{if } b_{i+1}=0, \\ 4D_i &\text{if }b_{i+1}=1. \end{cases} \\ &= \begin{cases} 0, &\text{if } b_{i+1}=0, \\ 2^2D_i &\text{if }b_{i+1}=1. \end{cases} \\ &= \begin{cases} 0, &\text{if } b_{i+1}=0, \\ [b_0;\dots;b_i;0;0] &\text{if }b_{i+1}=1. \end{cases} \end{align} $$ ### Linux 中的整數開根號 `int_sqrt` ```c unsigned long int_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; } ``` ## 藉由查表的平方根運算 > The algorithm is unusual in that it is **almost branch-free**, **requiring only a single final correction step**, and in that the number of required iterations is easily predicted in advance. The algorithm is efficient both at small scales and asymptotically, and represents an attractive compromise between speed and simplicity. 由 Mark Dickinson 整理的查表演算法，並在 Python 3.8 中採用。 * Heron's method 令 $f(x)=x^2-t$ ，由於 $f$ 為多項式函數，所以 $f$ 為光滑函數(無窮可微連續)，藉由牛頓法求解 $x^2-t=0$ 則一次微分 $f'(x)=2x$ 且牛頓法迭代表示為 $$ x_{n+1}=x_n-\dfrac{x^2_n-t}{2x_n}=\dfrac{1}{2}\left(x_n+\dfrac{t}{x_n}\right) $$ 在 [從 √2 的存在談開平方根的快速運算：牛頓法](https://hackmd.io/@sysprog/sqrt#%E7%89%9B%E9%A0%93%E6%B3%95) 一文中提及截斷誤差為 $$ \dfrac{(t-t_0)^2}{t_0} $$ 上述並非對整數平方根逼近，故為了處理整數取值逼近，上述牛頓法的迭代過程以下取整數函數改寫成 $$ g(a) = \left\lfloor \dfrac{1}{2}\left(a+\left\lfloor\dfrac{n}{a} \right\rfloor\right) \right\rfloor $$ 再來 Mark 照著以下思路至證明存在整數根號 1. 迭代生成的數列 $\lbrace a_i \rbrace$ 有下界 $\sqrt{t}$ 2. 數列有序 $a_i \leq a_{i+1}$ 3. 藉由良序原則(well-ordering principle)聲明最小 $a_i$ 存在 - [ ] Lemma 2. 任意正整數 $a$ 滿足 $$ \left\lfloor\sqrt n \right\rfloor \leq \left\lfloor \dfrac{1}{2}\left(a+\left\lfloor\dfrac{n}{a} \right\rfloor\right) \right\rfloor $$ * Proof of Lemma 1 一樣使用平方和公式。對任意正整數 $a$ 以及正數 $n$ 有 $$ 0 \leq (a-\sqrt n)^2=a^2-2a\sqrt n + n $$ 將根號部分移至不等式右方後，由於 $a$ 為正整數，兩邊同除 $2a$ 則 $$ \sqrt n \leq \dfrac{1}{2}\left(a+\dfrac{n}{a}\right) $$ 實際上就是對 $a$ 和 $\dfrac{n}{a}$ 的算術-幾何不等式，最後對不等式兩側做下取整數函數 $$ \left\lfloor\sqrt n \right\rfloor \leq \left\lfloor \dfrac{1}{2}\left(a+\dfrac{n}{a} \right) \right\rfloor = \left\lfloor \dfrac{1}{2}\left\lfloor\left(a+\dfrac{n}{a} \right) \right\rfloor \right\rfloor =\left\lfloor \dfrac{1}{2}\left(a+\left\lfloor\dfrac{n}{a} \right\rfloor\right) \right\rfloor $$ - [ ] Lemma 4. 講了初始 $a_0$ 滿足 $\left\lfloor\sqrt n \right\rfloor \leq a_0$ 時，以迭代式生成的數列 $a_i$ 有下界存在且為 $\left\lfloor\sqrt n \right\rfloor$ - [ ] Lemma 5. 證明少了 $\left\lfloor\sqrt n \right\rfloor=a_i$ only if $a_{i}\leq a_{i+1}$ 部分 > 非空正整數子集才能應用良序原則，證明中的細節少聲明這部分。參考 [PR #14: Modify the proof of Therorem 6](https://github.com/mdickinson/snippets/pull/14) 在文中 Mark 提及該算法計算出的該數列不一定需要收斂至 $\left\lfloor\sqrt n \right\rfloor$ ，文中使用 "**in the sense**" 一詞（[a general feeling or understanding](https://dictionary.cambridge.org/zht/%E8%A9%9E%E5%85%B8/%E8%8B%B1%E8%AA%9E-%E6%BC%A2%E8%AA%9E-%E7%B9%81%E9%AB%94/sense)）描述該算法最終會是一定數，但缺少相關證明與佐證，以下為原文： >Note that while the sequence is guaranteed to encounter $\left\lfloor\sqrt n \right\rfloor$ eventually, it's **not** necessarily true that the sequence **converges** to $\left\lfloor\sqrt n \right\rfloor$, in the sense that it's eventually constant. 綜上所述怎麼樣選取初始值是改進這個迭代法效率的要點。實作上對初始值的選取如下： 1. 令正整數初始值 $a_0$ 為 $n$ 2. 假設初始值 $a_0$ 要滿足 $a_0 \geq \left\lfloor\sqrt n \right\rfloor$ 3. 如不滿足 2. 則使用 Lemma 2 的結果 $\left\lfloor\sqrt n \right\rfloor \leq g(a_0)$ ，取代 $a_0$ 將 $g(a_0)$ 設置為新初始值  > 流程圖的繪製要依照國際標準 ISO-5807-1985 ```mermaid graph LR A("Let initial a0 be n") --> B{"Is a0 => floor(n) ?"} B -- yes --> C([do find_isqrt function]) B -- no --> D["replace a0 with g(a0)"] --> B ``` 想要不使用邏輯判斷分支， Mark 提出要從初始值的選取下手。 ### 選擇初始值的策略一個直觀的初始值選擇為超過 $\left\lfloor\sqrt n \right\rfloor$ 最小的 2 的冪(power of 2) 從 bit 的角度以 $2$ 的冪討論這個問題，先給定義： - [ ] Definition 7. For a nonnegative integer $n$, the *bit length* of $n$, which we'll write simply as $\text{length}(n)$, is the least nonnegative integer $e$ for which $n < 2^e$. 藉由不等式的推論過程 Mark 給出 $e$ 的最小下界：假設有一非負整數 $n$ 和 $e$ 滿足 $$ \left\lfloor\sqrt n \right\rfloor < 2^e $$ 根據該下取整函數不等式關係得到 $$ \left\lfloor\sqrt n \right\rfloor \leq \sqrt{n} < 2^e $$ 不等式兩側由於都是正數，故兩側各別平方不改變不等式方向 $$ n < 2^{2e} $$ 對 $n$ 取位元長度則不等式可以表示為 $$ \text{length}(n) \leq 2e $$ 由於該不等式關係 $\text{length}(n)$ 小於或等於 $2e$ ，當 $2e$ 減去 $1$ 後則只滿足小於關係 $$ \text{length}(n) < 2e + 1 $$ 不等式兩側同減 $1$ 後同除以 $2$ 則 $$ \dfrac{\text{length}(n) - 1}{2} < e $$ 再次從取整函數不等式關係得到 $$ \left\lfloor \dfrac{\text{length}(n) - 1}{2} \right\rfloor \leq \dfrac{\text{length}(n) - 1}{2} < e $$ 再次從 $\text{length}(n) \leq 2e$ 不等式關係出發，對兩側同時加 $1$ 後在同時除以 $2$ $$ \dfrac{\text{length}(n) + 1}{2} \leq e+\dfrac{1}{2} $$ 考慮不等式兩側同時做下取整函數，不等式左側 $e+\frac{1}{2}$ 由於 $e$ 為整數，故 $e+\frac{1}{2}$ 帶入下取整函數後會得到 $e$ 則該不等式關係為變成 $$ \left\lfloor \dfrac{\text{length}(n) + 1}{2} \right\rfloor \leq e $$ 以符號表式可以寫成 $$ \begin{align} \left\lfloor\sqrt n \right\rfloor < 2^e &\Rightarrow \sqrt{n} < 2^e \\ &\Rightarrow n < 2^{2e} \\ &\Rightarrow \text{length}(n) \leq 2e \\ &\Rightarrow \text{length}(n) < 2e + 1 \\ &\Rightarrow \left\lfloor \dfrac{\text{length}(n) - 1}{2} \right\rfloor < e \end{align} $$ 其中 $$ \text{length}(n) \leq 2e \Rightarrow \left\lfloor \dfrac{\text{length}(n) + 1}{2} \right\rfloor \leq e $$ 具體的來說，該等不等式之間是邏輯等價的 $$ \begin{align} \left\lfloor\sqrt n \right\rfloor < 2^e &\iff \sqrt{n} < 2^e \\ &\iff n < 2^{2e} \\ &\iff \text{length}(n) \leq 2e \\ &\iff \text{length}(n) < 2e + 1 \\ &\iff \left\lfloor \dfrac{\text{length}(n) - 1}{2} \right\rfloor < e \\ &\iff \left\lfloor \dfrac{\text{length}(n) + 1}{2} \right\rfloor \leq e \end{align} $$ 僅需要證明 $\left\lfloor \dfrac{\text{length}(n) + 1}{2} \right\rfloor \leq e$ 和 $\left\lfloor\sqrt n \right\rfloor < 2^e$ 邏輯等價，以及 $\left\lfloor \dfrac{\text{length}(n) - 1}{2} \right\rfloor < e$ 和 $\left\lfloor\sqrt n \right\rfloor < 2^e$ 邏輯等價。只要確認該兩命題成立，便可以由 $\left\lfloor\sqrt n \right\rfloor < 2^e$ 推至其他不等式關係。由於我們已經證明完該不等式 $\left\lfloor \dfrac{\text{length}(n) + 1}{2} \right\rfloor \leq e$ 為 $\left\lfloor\sqrt n \right\rfloor < 2^e$ 的充分條件，故只要證明該不等式同時也為 $\left\lfloor\sqrt n \right\rfloor < 2^e$ 的必要條件即可。由於 $e$ 為整數， $e=\left\lfloor e + \frac{1}{2} \right\rfloor$ ，故 $\left\lfloor \dfrac{\text{length}(n) + 1}{2} \right\rfloor \leq e$ 可以寫成 $$ \begin{align} \left\lfloor \dfrac{\text{length}(n) + 1}{2} \right\rfloor \leq \left\lfloor e + \frac{1}{2} \right\rfloor &\Rightarrow \dfrac{\text{length}(n) + 1}{2} \leq e + \frac{1}{2} \\ &\Rightarrow \text{length}(n) + 1 \leq 2e +1 \\ &\Rightarrow \text{length}(n) \leq 2e \\ &\Rightarrow n \leq 2^{2e} \\ &\Rightarrow \sqrt{n} \leq 2^{e} \\ &\Rightarrow \left\lfloor\sqrt n \right\rfloor < 2^e \end{align} $$ 證明 $\left\lfloor \dfrac{\text{length}(n) - 1}{2} \right\rfloor < e$ 和 $\left\lfloor\sqrt n \right\rfloor < 2^e$ 邏輯等價方式類似，在此不贅述。所以該等不等式之間彼此間邏輯等價。該等邏輯等價關係告訴我們：該初始值選擇條件可以從 $n$ 的位元長度計算 $$ \left\lfloor \dfrac{\text{length}(n) + 1}{2} \right\rfloor $$ ### 鄰近平方根 near square root - [ ] 定義9：假設 $n$ 為一正整數。則當 $(a-1)^2<n<(a+1)^2$ 稱正整數 $a$ 為 $n$ 鄰近平方根。此定義可以等價表示：如果 $a=\left\lfloor\sqrt x \right\rfloor$ 或 $\left\lceil\sqrt x \right\rceil$ 則 $a$ 為一鄰近平方根。假如有一 $n$ 的鄰近平方根 $a$，則 $a^2 \leq n$ 時，$n$ 的整數根號為 $a$；則 $a^2 > n$ 時， $n$ 的整數根號為 $a-1$ 。總得來說，一個找出鄰近平方根的演算法，提供了一個計算整數平方根的方式。此演算法的想法： > Our algorithm is based on the idea of “lifting” a near square root of a quotient of $n$ to a near square root of $n$.  此想法藉由[帶餘數除法（division with remainder）](https://en.wikipedia.org/wiki/Euclidean_division)，將被開根號數 $n$ 寫成一個除數與商數的乘積，加上一個餘數，該除數以某個正整數平方表示且除數小於或等於該被開根號數 $n$。同時，當 $n$ 為完全平方時，除數會於 $n$；否則，藉由除數與商數的乘積近似該被開根號數 $n$，除法結果近似於該商數。對此被除數與除數開根號後相除，商數開根號會近似該除式。再次從帶餘數除法得知，商數開根號乘上除數會近似 $n$ 開根號。假設 $j$ 為一正整數，滿足 $j^2\leq n$ 和 $b$ 為 $\left\lfloor \frac{n}{j^2} \right\rfloor$ 一鄰近平方根。換言之，對 $\left\lfloor \frac{n}{j^2} \right\rfloor$ 開平方根會得到 $b$ ，也就是說可以將 $b$ 近似表示為 $$ b\approx \frac{\sqrt{n}}{j} $$ 則可以將 $j,b,n$ 的關係寫成 $$ jb \approx \sqrt{n} $$ 也就是說，對於求解 $jb=\sqrt{n}$ 可以得到一近似平方根 $\sqrt{n}$，則照 Heron's method ，寫出一改進的近似方法 $$ g(jb) = \left\lfloor \dfrac{1}{2}\left(jb+\left\lfloor\dfrac{n}{jb} \right\rfloor\right) \right\rfloor \approx \sqrt{n} $$ 原文中給出一個限制在於 $j^2 \leq n$ ，假設 $j$ **不要太大**的情況下，使用該近似方法可以求得 $n$ 的鄰近平方根。具體的作法是令 $j=2k$ 為一偶數，則根據假設 $b$ 為 $\left\lfloor \frac{n}{j^2} \right\rfloor=\left\lfloor \frac{n}{4k^2} \right\rfloor$ 的鄰近平方根，其中 $2kb$ 為 $\sqrt{n}$ 的近似值。則替換掉 $j$ 後近似方法為 $$ \begin{align} g(2kb) &= \left\lfloor \dfrac{1}{2}\left(2kb+\left\lfloor\dfrac{n}{2kb} \right\rfloor\right) \right\rfloor \\ &= \left\lfloor kb + \dfrac{1}{2}\left\lfloor\dfrac{n}{2kb} \right\rfloor \right\rfloor \\ &?= kb + \left\lfloor\dfrac{n}{4kb} \right\rfloor \end{align} $$ $$ \left\lfloor \dfrac{1}{2}\left\lfloor\dfrac{n}{2kb} \right\rfloor\right\rfloor = \left\lfloor\dfrac{n}{4kb} \right\rfloor ? $$ - [ ] 定理10：假設 $n,k$ 為正整數，且滿足 $4k^2\leq n$。令 $b$ 表示為 $\left\lfloor\dfrac{n}{4kb} \right\rfloor$ 的鄰近平方根。則 $n$ 的鄰近平方根表示為 $$ a = kb + \left\lfloor\dfrac{n}{4kb} \right\rfloor $$ * 證明根據鄰近平方根的定義，將 $b$ 與 $\left\lfloor\dfrac{n}{4kb} \right\rfloor$ 有一不等式關係 $$ (b-1)^2 \leq \left\lfloor\dfrac{n}{4kb} \right\rfloor \leq (b+1)^2 $$ 需修正 example 11 typo floor(n/k^2) --> floor(n/4k^2) ## 測試使用 Python script 先生成出 $1$ 至 $2^{32}-1$ 的平方，儲存成一 `perfect_square32.txt`。再使用 C 讀取該 `perfect_square32.txt` 執行十分逼近法開根號。 ### 完全平方根 ``` Performance counter stats for './test' (40 runs): 29.30 msec task-clock # 0.979 CPUs utilized ( +- 0.54% ) 1 context-switches # 34.128 /sec ( +- 8.01% ) 0 cpu-migrations # 0.000 /sec 187 page-faults # 6.382 K/sec ( +- 0.07% ) 49,244,519 cycles # 1.681 GHz ( +- 0.10% ) (72.81%) 452,034 stalled-cycles-frontend # 0.92% frontend cycles idle ( +- 0.91% ) (86.22%) 10,980,482 stalled-cycles-backend # 22.30% backend cycles idle ( +- 0.67% ) (86.40%) 147,406,087 instructions # 2.99 insn per cycle # 0.07 stalled cycles per insn ( +- 0.04% ) (86.50%) 34,442,449 branches # 1.175 G/sec ( +- 0.71% ) (86.41%) 170,372 branch-misses # 0.49% of all branches ( +- 1.36% ) (81.65%) 0.029942 +- 0.000165 seconds time elapsed ( +- 0.55% ) ``` ### 不同整數數值輸入 ## 參考資料 * [2024q1 第 3 週測驗題](https://hackmd.io/@sysprog/linux2024-quiz3) * [2024q1 第 4 週測驗題](https://hackmd.io/@sysprog/linux2024-quiz4) * [Linux, int_sqrt](https://github.com/torvalds/linux/blob/8f2c057754b25075aa3da132cd4fd4478cdab854/lib/math/int_sqrt.c#L49) * [從 √2 的存在談開平方根的快速運算](https://hackmd.io/@sysprog/sqrt) * [Analysis of the lookup table size for square-rooting](https://web.ece.ucsb.edu/~parhami/pubs_folder/parh99-asilo-table-sq-root.pdf) * [Timing square root](https://web.archive.org/web/20210208132927/http://assemblyrequired.crashworks.org/timing-square-root/) * [Wikipedia: Digit-by-digit calculation](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit-by-digit_calculation) * [Square roots by subtraction by Frazer Jarvis](http://www.afjarvis.org.uk/maths/jarvisspec02.pdf) * [查表 isqrt](https://github.com/mdickinson/snippets/tree/main/papers/isqrt) * [perf.wiki](https://perf.wiki.kernel.org/index.php/Main_Page)