# 整數開平方根 [Quiz2-2](https://hackmd.io/@sysprog/linux2025-quiz2#%E6%B8%AC%E9%A9%97-2)
contributed by <`MikazukiHikari`>
## sqrti 的原理
### count leading zero
首先是 `clz`，有關 `clz` 的詳細作法請參考 [count leading zero](https://hackmd.io/@sysprog/linux2025-quiz2#%E6%B8%AC%E9%A9%97-2) ，這邊不再贅述，我們只需要知道使用 `clz` 的目的是計算輸入值的前導零數量(位數)以縮小平方根的範圍，從而減少迭代的次數並加速 `sqrti` 的過程，也就是說其實不做 `clz` 也能得到答案，只是耗時會比較多。
#### 為什麼輸入值的位數和平方根有關?
這邊就需要舉例了，這也和我們要求 `shift` 是偶數，`m` 是 power of 4 有關：
$(4)_{10} = (100)_2$ 的平方是 $(16)_{10} = (10000)_2$ => 5 位數
$(8)_{10} = (1000)_2$ 的平方是 $(64)_{10} = (1000000)_2$ => 7 位數
可以發現當求一個數的平方根時，若它在二進位是 5~6 位數，則我們可以確定它的平方根必然在 4 到 7 之間。
也就是說透過 `clz` 可以計算輸入值的位數，再來我們就可以用 `m` 表達輸入值的平方根的範圍。
假設 $m = a^2$ ( $a = 2^n$ )，輸入值就在 $m$ ~ $4m$ ，或者說 $a^2$ ~ $(2a)^2$ 的範圍內，而平方根則會在 $a$ ~ $2a$ 的範圍內，這裡也能解釋為何 `m` 是 power of 4，得出可能的平方根的範圍也方便我們使用後面的 Digit-by-digit Method 去獲得最終平方根。
### Digit-by-digit Method
先來講解 Digit-by-digit 的概念，舉個例子：假設輸入為 226 ，我們也已透過 `clz` 得知它介於$8^2$ 到 $16^2-1$，而其平方根是在 8~15 之間；我們便會從 8 開始，逐步逼近。若把從 0 逼近 8 當成第一步的話，第二步就是$8>>1 = 4$，因此我們會計算它逼近 4 後的平方，也就是 $(8+4)^2 = 144$，因輸入 (226) >= 144，代表其平方根大於等於 $8+4=12$；接下來第三步是$4>>1 = 2$，計算它再逼近 2 後的平方，也就是 $(12+2)^2 = 196$，輸入 (226) >= 196，代表其平方根大於等於 $12+2=14$；最後第四步$2>>1 = 1$，計算它再逼近 1 後的平方，也就是 $(14+1)^2 = 225$，輸入 (226) >= 225，代表其平方根大於等於 $14+1=15$，最後我們得出了 226 的整數平方根是 15 。
#### 如何算出平方差
知道了 Digit-by-digit 的概念，我們就能透過同樣的方法算出任意 `uint64_t` 的數字的整數開平方根。然而，最後一個問題：如何算出下一步逼近的值的平方呢？
只要算出平方差，再加上原來的平方數便可得出下一步逼近的值的平方。
```c
while (m) {
uint64_t b = y + m;
y >>= 1;
if (x >= b) {
x -= b;
y += m;
}
m >>= 2;
}
return y;
```
參考部分 `sqrti` 的程式碼，已知 `m` 是每次逼近的值的平方，`x` 為輸入，而 `b` 其實就是平方差，每次迭代輸入`x` 減平方差 `b` 其實就等於反向的加平方差。
這邊就要用到一個公式：$(a+c)^2 - a^2 = a^2 + 2ac + c^2 - a^2$，平方差就是 $2ac + c^2$，到這邊就會發現它對應程式碼中的 `b = y + m`，`y` 是 $2ac$，`m` 是 $c^2$。
繼續使用剛才的 226 當例子：
```
第一次迭代 (逼近8):
(0+8)^2 - 0^2 = 0^2 + 2*0*8 + 8^2 -0^2 => m=8^2=64; y=2*0*8=0
第二次迭代 (逼近4):
(8+4)^2 - 8^2 = 8^2 + 2*8*4 + 4^2 -8^2 => m=4^2=64>>2; y=2*8*4=2*(0+8)*(8>>1)
第三次迭代 (逼近2):
(12+2)^2 - 12^2 = 12^2 + 2*12*2 + 2^2 -12^2 => m=2^2=16>>2; y=2*12*2=2*(8+4)*(4>>1)
第四次迭代 (逼近1):
(14+1)^2 - 14^2 = 14^2 + 2*14*1 + 1^2 -14^2 => m=1^2=4>>2; y=2*14*1=2*(12+2)*(2>>1)
```
`m` 是 $c^2$ 很好理解，相當於每次逼近的數的平方，每一輪的最後 `m >>= 2` 代表下次逼近的數為這次的一半，其平方則為四分之一(>>2)；至於 `y`為 $2ac$，舉 226 的第二次和第三次迭代，前者的 `y` 為 $2*8*4$ 之後執行 `y >>= 1` ，變成 $2*8*(4>>1)$，再執行 `y += m` ，即加上`m` 也就是第二次的 $c^2=4^2=4*4=2*4*(4>>1)$，等於 $2*(8+4)*(4>>1)$，正好就是第三次迭代的 $2ac=2*12*2$，因為 $a$ 每次都會多了上一次逼近的數，$c$ 又等於上一次逼近的數的一半，因此 $2ac$ 就會剛好等於上一次的 `y>>1` 加上 `m`。
而最後一次迭代的 `y` 為 $2ac$，`y >>= 1` 後為 $ac$，又最後一次的 `c` 必為 1 ，`y` 為 $a$，最後再根據比較 `x` 和 `b` 判斷是否要 +1 ，最終 `y` 為整數平方根。
整段程式最酷的是，我們可以只用二進位的 shift 和加減來表達整個過程，包含裡面的所有乘法和平方。
## 改寫程式碼，不使用 clz / clz64，不能用乘法或除法
前面說到 `clz` 僅用於求輸入之位數，換句話說所有能求輸入之位數的方法都可替代 `clz` 。
### 方法 1：逐位檢查（最簡單但較慢）
```c
int count_leading_zeros(uint64_t x) {
if (x == 0) return 64;
int count = 0;
for (uint64_t mask = 1ULL << 63; mask > 0; mask >>= 1) {
if (x & mask) break;
count++;
}
return count;
}
```
從最高位（第63位）開始逐位檢查，直到遇到第一個 1。
### 方法 2：二分搜尋（高效位操作）
```c
int count_leading_zeros(uint64_t x) {
if (x == 0) return 64;
int count = 0;
if (x <= 0x00000000FFFFFFFF) { count += 32; x <<= 32; }
if (x <= 0x0000FFFFFFFFFFFF) { count += 16; x <<= 16; }
if (x <= 0x00FFFFFFFFFFFFFF) { count += 8; x <<= 8; }
if (x <= 0x0FFFFFFFFFFFFFFF) { count += 4; x <<= 4; }
if (x <= 0x3FFFFFFFFFFFFFFF) { count += 2; x <<= 2; }
if (x <= 0x7FFFFFFFFFFFFFFF) { count += 1; }
return count;
}
```
通過二分法逐步縮小範圍，每次判斷高位是否全零。
### 方法 3：查表法（空間換時間），可和二分搜尋合併
```c
static const uint8_t leading_zero_table[256] = {
8,7,6,6,5,5,5,5,4,4,4,4,4,4,4,4, // 0x00-0x0F
3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3, // 0x10-0x1F
// ... and so on
};
int count_leading_zeros(uint64_t x) {
if (x == 0) return 64;
int count = 0;
if (x <= 0x00000000FFFFFFFF) { count += 32; x <<= 32; }
if (x <= 0x0000FFFFFFFFFFFF) { count += 16; x <<= 16; }
if (x <= 0x00FFFFFFFFFFFFFF) { count += 8; x <<= 8; }
count += leading_zero_table[(x >> 56)];
return count;
}
```
將 64 位元分為 8 位元一組，先用二分搜尋判定前導1的範圍，再通過查表快速定位前導零數。
## 如何加速 sqrti
使用 `clz` + [牛頓迭代法](https://hackmd.io/@sysprog/sqrt#%E7%89%9B%E9%A0%93%E6%B3%95)
還是一樣先用 `clz` 計算輸入值的前導零數量(位數)以縮小平方根的範圍，從而減少迭代的次數並加速 `sqrti` 的過程，然而這次是用於牛頓迭代法的初始值；再使用牛頓迭代法，利用誤差的平方遞減特性，指數級地使平方根的有效位數翻倍，雖然每次迭代的計算量稍大（需一次除法），但總次數極少，整體速度更快，數學分析請參考[如何檢驗「加速」版本的 sqrti](https://hackmd.io/@DboOgKS6RtOmMLCltFsBig/SkS4ZDUeeg#%E5%A6%82%E4%BD%95%E6%AA%A2%E9%A9%97%E3%80%8C%E5%8A%A0%E9%80%9F%E3%80%8D%E7%89%88%E6%9C%AC%E7%9A%84-sqrti)。
```c
uint64_t sqrt_newton(uint64_t a) {
if (a == 0 || a == 1) return a;
uint64_t m = 0;
int shift = (63 - clz64(a) + 1) >> 1; // Correct bit-width calculation
m = 1ULL << shift;
uint64_t x = m;
uint64_t y = (x + a / x) >> 1; // Newton's method formula
do{
x = y;
y = (x + a / x) >> 1;
// Check: exit early if the next step does not change
if (y >= x && (x + a / x) >> 1 == y) break;
} while (y < x);
return (y > x) ? x : y;
}
```
### 相比純粹的牛頓法的關鍵改進
* 改善初始猜測值：
基於輸入值的位數直接估算平方根位數(`clz`)，大幅減少迭代次數。
* 殘差檢查：
防止因整數除法地板效應導致的無效迭代。
* 判斷輸入是否為完全平方數：
由於牛頓法正常來說得到的平方根會是 `ceil_sqrti`，可將結果減 1 得到 `floor_sqrti` ，然而需先透過 `x` 和 `y` 來判斷輸入是否為完全平方數後，才能決定是否能減 1。
### 為何牛頓法能比逐位法還快呢?
逐位法雖然能不使用乘法或除法，然而由於其需要一位一位的逼近，在 `uint64_t` 的輸入的條件之下，最壞可能會遇到例如 $2^{64}-1$ 這種高達 64 位數的輸入，導致其必須迭代 32 次才能將平方根從第 31 位逼近到第 0 位，因此在這邊會消耗大量的時間。
而牛頓法雖然每次迭代的計算量稍大（需一次除法），但其理論上最多需要6次迭代(可參考下一段)，因此整體速度更快。
## 如何檢驗「加速」版本的 sqrti
牛頓迭代法在計算整數平方根時，理論上最多需要6次迭代的依據來自其二次收斂性和 64位元整數的位數限制。以下是詳細分析：
1. 誤差遞推公式
設真實平方根為 $s=⌊a⌋$，當前估算值為 $x_k$，誤差定義為：
$ϵ_k=x_k-s$
在整數牛頓迭代中：
$x_{k+1}=⌊\frac{x_k+⌊\frac{a}{x_k}⌋}{2}⌋$
誤差遞推滿足不等式：
$ϵ_{k+1}\le \frac {ϵ_k^2}{2s}$
2. 64位整數的參數設定
最大輸入值：$a=2^{64}-1$
真實平方根：$s=2^{32}-1$（32位二進制數）
初始誤差：若基於CLZ的初始猜測 $x_0=2^{32}$，則誤差：
$ϵ_0=x_0-s=1$
3. 收斂次數推導
代入誤差公式，迭代次數 $k$ 滿足：
$ϵ_{k}\le \frac {ϵ_0^{2^k}}{(2s)^{2^k-1}}$
對於64位整數，代入 $ϵ_0=1$ 和 $s=2^{32}$：
$ϵ_{k}\le {\frac {1^{2^k}}{(2^{33})^{2^k-1}}}=\frac {1}{(2^{33})^{2^k-1}}$
要求誤差小於 1 （即精確到整數）：
$2^{33(2^k-1)}>1$ => $33(2^k-1)>0$ => $2^k>1$ => $k>0$
但由於整數除法的地板效應，考慮到整數除法的截斷效應，增加 1 次迭代作為安全邊界。考慮最壞情況（初始誤差為 $2^{32}$ ），推導修正為：
$2^k \ge log_2(32)+1=6$ => $k=6$
幾何級數收斂驗證：
| 迭代次數 | 誤差上界 | 有效二進制位數 |
| -------- | -------- | -------- |
| 0 | $2^{32}$ | 0位 |
| 1 | $2^{31}$ | 1位 |
| 2 | $2^{29}$ | 3位 |
| 3 | $2^{25}$ | 7位 |
| 4 | $2^{17}$ | 15位 |
| 5 | $2^{1}$ | 31位 |
| 6 | $2^{-31}$ | 63位 （誤差<1） |
每次迭代誤差的有效位數幾乎翻倍。
舉例：若初始誤差 $ϵ_0=2^{32}$，$s≈2^{32}$（對64位整數），根據誤差遞推不等式 $ϵ_{k+1}\le \frac {ϵ_k^2}{2s}$：
* 第1次迭代：$ϵ_{1}\le 2^{32}$
* 第2次迭代：$ϵ_{2}\le \frac{(2^{32})^2}{2*2^{32}}=2^{31}$
* 第3次迭代：$ϵ_{3}\le \frac{(2^{31})^2}{2*2^{32}}=2^{29}$
* ...
* 第k次迭代：$ϵ_{k}\le 2^{32-2^{k-1}}$
當 $ϵ_k<1$ 時，已經精確到整數，這通常只需6次（因為$2^{32-2^{6-1}}=1$)。
然而由於搭配 `clz` ，得出較為精確的牛頓迭代法初始值，且整數除法導致實際迭代次數少於理論值（因地板效應提前終止），因此大部分的實際迭代次數可能都小於6次。
下面的測試程式輸出的牛頓法和逐位法的迭代數也可證明牛頓法的迭代數遠小於逐位法。
worse case 應為 64 位元，也就是 `clz` 為 0 的數字，逐位法需迭代 32 次(`n_digit`)；而牛頓法只需迭代最多 6 次
測試原始碼如下，使用隨機生成的 1000000 個值測試：
```c
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <time.h>
#define N 1000000
uint64_t data[N];
int clz64(uint64_t x) {
if (x == 0) return 64;
return __builtin_clzll(x);
}
uint64_t sqrt_newton(uint64_t a) {
if (a == 0 || a == 1) return a;
int shift = (63 - clz64(a) + 1) >> 1;
uint64_t x = 1ULL << shift;
uint64_t y = (x + a / x) >> 1;
do {
x = y;
y = (x + a / x) >> 1;
} while (y < x);
return (y > x) ? x : y;
}
uint64_t sqrt_digit(uint64_t x) {
if (x == 0 || x == 1) return x;
uint64_t m, y = 0;
int shift = (63 - 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;
}
void prepare_data() {
srand(0);
for (int i = 0; i < N; i++) {
data[i] = ((uint64_t)rand() << 32) | rand();
}
}
void test_newton_all() {
volatile uint64_t dummy = 0;
for (int i = 0; i < N; i++) {
dummy = sqrt_newton(data[i]);
}
}
void test_digit_all() {
volatile uint64_t dummy = 0;
for (int i = 0; i < N; i++) {
dummy = sqrt_digit(data[i]);
}
}
int main(void) {
prepare_data();
#ifdef ONLY_NEWTON
test_newton_all();
#endif
#ifdef ONLY_DIGIT
test_digit_all();
#endif
return 0;
}
```
以下是使用 perf 個別分析後的結果
clz + 牛頓法：
```
$ gcc -DONLY_NEWTON -O2 -o sqrt_test sqrt_test.c
$ perf stat ./sqrt_test
Performance counter stats for './sqrt_test':
76.76 msec task-clock # 0.995 CPUs utilized
0 context-switches # 0.000 /sec
0 cpu-migrations # 0.000 /sec
2,004 page-faults # 26.109 K/sec
286,569,753 cycles # 3.734 GHz
196,498,060 instructions # 0.69 insn per cycle
39,006,486 branches # 508.192 M/sec
445,043 branch-misses # 1.14% of all branches
0.077175204 seconds time elapsed
0.075169000 seconds user
0.002004000 seconds sys
```
clz + 逐位法：
```
$ gcc -DONLY_DIGIT -O2 -o sqrt_test sqrt_test.c
$ perf stat ./sqrt_test
Performance counter stats for './sqrt_test':
145.69 msec task-clock # 0.993 CPUs utilized
4 context-switches # 27.456 /sec
1 cpu-migrations # 6.864 /sec
2,004 page-faults # 13.755 K/sec
591,639,505 cycles # 4.061 GHz
379,322,939 instructions # 0.64 insn per cycle
98,156,441 branches # 673.746 M/sec
16,212,165 branch-misses # 16.52% of all branches
0.146752388 seconds time elapsed
0.139861000 seconds user
0.006993000 seconds sys
```
將兩者都執行，查看各函式 CPU 佔用比例：
```
$ perf report
61.70% sqrt_test sqrt_test [.] test_digit_all
22.02% sqrt_test sqrt_test [.] test_newton_all
11.77% sqrt_test libc.so.6 [.] __random
1.50% sqrt_test libc.so.6 [.] __random_r
0.69% sqrt_test libc.so.6 [.] rand
0.32% sqrt_test sqrt_test [.] prepare_data
0.26% sqrt_test [unknown] [k] 0xffffffff8c058511
0.25% sqrt_test sqrt_test [.] rand@plt
0.16% sqrt_test [unknown] [k] 0xffffffff8cea9a0c
0.14% sqrt_test [unknown] [k] 0xffffffff8c046c82
0.13% sqrt_test [unknown] [k] 0xffffffff8c07812f
0.13% sqrt_test [unknown] [k] 0xffffffff8c05015f
0.13% sqrt_test [unknown] [k] 0xffffffff8d001192
0.13% sqrt_test [unknown] [k] 0xffffffff8c080f20
0.13% sqrt_test [unknown] [k] 0xffffffff8c05011a
0.13% sqrt_test [unknown] [k] 0xffffffff8d001248
0.13% sqrt_test [unknown] [k] 0xffffffff8c0edbac
0.13% sqrt_test [unknown] [k] 0xffffffff8c07c6e9
0.12% sqrt_test [unknown] [k] 0xffffffff8c080e39
0.01% perf-exec [unknown] [k] 0xffffffff8bfb8e80
0.00% perf-exec [unknown] [k] 0xffffffff8bcbe0e3
```
可以發現：
| 函式名稱 | CPU 佔用比例 | 解釋 |
| ----------------- | ---------- | ------------------------------- |
| `test_digit_all` | **61.70%** | 逐位法總共消耗最多 CPU 時間，是主要瓶頸 |
| `test_newton_all` | 22.02% | 牛頓法只佔不到 1/3 的時間，明顯較快 |
| `__random` 等 | 14%+ | 隨機數生成也占了可觀比例，建議固定輸入以排除這部分對實驗的干擾 |
### 牛頓法會迭代幾次？
```c
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#define N 1000000
uint64_t data[N];
// variables for counting
int iteration_histogram[100] = {0}; // Assume the number of iterations does not exceed 100.
int max_iter = 0;
int total_iter = 0;
int clz64(uint64_t x) {
if (x == 0) return 64;
return __builtin_clzll(x);
}
uint64_t sqrt_newton(uint64_t a) {
if (a == 0 || a == 1) return a;
int shift = (63 - clz64(a) + 1) >> 1;
uint64_t x = 1ULL << shift;
uint64_t y = (x + a / x) >> 1;
int n = 1;
do {
x = y;
y = (x + a / x) >> 1;
n++;
if (y >= x && (x + a / x) >> 1 == y) break;
} while (y < x);
if (y > x) y--;
// Count the number of iterations.
if (n < 100) iteration_histogram[n]++;
total_iter += n;
if (n > max_iter) max_iter = n;
return y;
}
void prepare_data() {
srand(0);
for (int i = 0; i < N; i++) {
data[i] = ((uint64_t)rand() << 32) | rand();
if (data[i] == 0) data[i] = 1;
}
}
void test_newton_all() {
for (int i = 0; i < N; i++) {
sqrt_newton(data[i]);
}
printf("The number of iterations（number : count）\n");
for (int i = 1; i <= max_iter; i++) {
if (iteration_histogram[i] > 0)
printf("%2d : %d\n", i, iteration_histogram[i]);
}
printf("The number of max iterations：%d\n", max_iter);
printf("The number of average iterations：%.4f\n", (double)total_iter / N);
}
int main() {
prepare_data();
test_newton_all();
return 0;
}
```
輸出如下：
```
The number of iterations（number : count）
2 : 60
3 : 17789
4 : 289946
5 : 692025
6 : 180
The number of max iterations：6
The number of average iterations：4.6745
```
將統計結果作圖：
這個結果非常清楚地展示了牛頓法在 1000000 筆隨機資料上：
* 平均只需要約 4.67 次迭代
* 最多只需 6 次迭代
* 大多數輸入集中在 第 4～5 次迭代 完成
這證明了牛頓法在實務上非常快速、收斂次數穩定，且比逐位法效率高很多（後者通常迭代固定次數，例如 32 次）。
### 牛頓法與逐位法的 worst case 測試比較
我們使用上一節找到的 180 個對牛頓法與逐位法皆是 worst case 的例子執行 `180*1000` 次去比較兩者的差距：
```c
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#define N 180*1000
int clz64(uint64_t x) {
if (x == 0) return 64;
return __builtin_clzll(x);
}
uint64_t sqrt_newton(uint64_t a) {
if (a == 0 || a == 1) return a;
int shift = (63 - clz64(a) + 1) >> 1;
uint64_t x = 1ULL << shift;
uint64_t y = (x + a / x) >> 1;
do {
x = y;
y = (x + a / x) >> 1;
} while (y < x);
return (y > x) ? x : y;
}
uint64_t sqrt_digit(uint64_t x) {
if (x == 0 || x == 1) return x;
uint64_t m, y = 0;
int shift = (63 - 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;
}
uint64_t input_data[180] = {
0x7ed344d32f378c0f, 0x7f31d0aa5c5733bf, 0x7c8b6314430c9035, 0x7cc9b36d25619e58,
0x7fda725056f807d8, 0x7f8ab4610ea18d0a, 0x7f5ff0c91fb79424, 0x7d6249c37885d25e,
0x7c9a527113f2046b, 0x0791b5461e256b57, 0x7ed7202a0a102827, 0x7dec0fcc27debde0,
0x7f2ca4d014d241d7, 0x70b468d10dec9ddd, 0x7cef688d2aa916d7, 0x1fb678681528d887,
0x77e99a761929a85c, 0x7fc0afef0c2cb445, 0x6f8c30cf10f74355, 0x7fef05044db09608,
0x7bb1a19522fc9dd1, 0x7ba195e658d941a7, 0x7f867baf477df979, 0x20a27ac9287289d2,
0x2233c5bd41a76c35, 0x7e7ce3990d446713, 0x7c932d7f40bdef41, 0x7fef7b5f272351ea,
0x7751403b5c3fad56, 0x7f50fdb6786f1b73, 0x711b8bb97df77bb2, 0x78efb4cb15b0a325,
0x77973450659c3259, 0x790e592c7dcaf3f2, 0x7e49a0d2206e405a, 0x7d4559256a9e5230,
0x7db6b881322bc30b, 0x2159d4233fa70a7f, 0x7931cc587beb2875, 0x2274ad2a3a1dafea,
0x1c254d52538241d1, 0x1f12488f24cff92b, 0x7af36a1334ece5c9, 0x7d490bcc7d0229bf,
0x7b3a25cf4053e9ff, 0x7890a3cc1186a59c, 0x1fc9a44170629652, 0x769f7519406cf3d9,
0x7c6729476d013ea6, 0x7e7f94fb6efcfc76, 0x7cd4d50b0c517ade, 0x7fb864cc12ae49ae,
0x75048b9948c8b314, 0x7f1f3a082a9874a0, 0x7619e71910e46dcc, 0x022797a52ba69961,
0x7d816ff713bab46f, 0x7c4f42d52101c04f, 0x7fade98a39b861f0, 0x7818b1655c43ba0a,
0x7e75ac65173fb5fd, 0x209afe85111f6c94, 0x7ca5482e76255f39, 0x7a9c476e5da66a8c,
0x7cd771e871b43be2, 0x7dad557771e57c04, 0x20380b9779ea368e, 0x7ab1cb0c3c2c08cc,
0x7e7f1f1e4015f7d9, 0x7fcdd6e0471f3324, 0x7ec982bc69b0ad49, 0x20d2010e0229a26b,
0x7c0a71551b2e7ff1, 0x7e18d75e34d215e4, 0x1fe2bcca44300eac, 0x7eec53af0703bc85,
0x7cffbf484174dab4, 0x20b96fbe35e9b67d, 0x7d7f2fcd1926a01d, 0x78ebbd387226e3b7,
0x7c622041461977bb, 0x7c3f473c7caeeba4, 0x22d082d3568cb903, 0x7f9a9c6b70917e01,
0x7eae3baf2ed254bd, 0x7f95281e5ff4980d, 0x081160cf0d1d6ee3, 0x7f4b5e1936275899,
0x7af0718672dbc4a6, 0x7f97628f7cbbd54d, 0x7d665da4063b2cd6, 0x75b3d4a25d192261,
0x203d8c5722dadc31, 0x79b406320a91e571, 0x776fc903237ffd67, 0x7ce9d7ff4a5eef86,
0x7f70a3af44706ee2, 0x7691872e716aad2f, 0x22714f5b7a640ca0, 0x1fd9e54229f2520f,
0x75e3babc3062169c, 0x20eaccf83b4dcd78, 0x7c9bffa506ebf76f, 0x7c5b68443260a1fb,
0x7f31c9875111090d, 0x7e1785290612ac9f, 0x1fba498362afcd09, 0x1f9d5e076934693a,
0x2433302f0d282bbc, 0x76178e555e40077e, 0x7e5b5fdc0e5a1fb7, 0x7d047aca233135d5,
0x7f7e397656cd8491, 0x7a9966a71a371600, 0x1e6e0e8160dc8a9a, 0x7d633dbb5fe6cc33,
0x7fc98ce76adc1798, 0x71dc10666b357fdd, 0x73f4647220ff5d8d, 0x7ee5e57f10055957,
0x7c45e401540cf907, 0x7e9556a8185ac28e, 0x1ea483351a77758e, 0x79fb7c8f2cabc813,
0x7cb3186f0718a1a3, 0x7e953c033a769222, 0x7957226f30248999, 0x7749d85614f74b24,
0x6edbf71b39067296, 0x72c668495d08d113, 0x209a390823634bbc, 0x7f60f29a692ea4cd,
0x7f238de60f30abbd, 0x7f8056582ba9d0ce, 0x7f0bfc6d3328a5a9, 0x7ffee7f47d1521e4,
0x7a5636466b6e105c, 0x7e04d20057c9bb90, 0x744784813b162e4e, 0x20f807d32fd84771,
0x7ca823c13d6d4940, 0x7a31c3253f7a98f4, 0x7e35c2d02e8216bf, 0x216d8f7b08b28737,
0x7c332f3279e411b6, 0x7fd5c6e90182c53a, 0x1f5cba124b97dda3, 0x7ef5dbf64f509e96,
0x1fd39f370b1fafcc, 0x1fb22a4043f51c4b, 0x7e1879d803b61424, 0x7abf2cf142c81ef3,
0x7beba50811c09258, 0x77dcab3a2cdcb213, 0x79f08c8b44b4f626, 0x7f5ccd5a2ad46a65,
0x7da2c8be3c17c34b, 0x1e49467677b645e9, 0x212dd7622681ba9c, 0x7ccaf2f47f034ed3,
0x76512316378a2af0, 0x7db6424e57b2b84d, 0x7ea1666245b2f7c0, 0x7992b4fc273d9214,
0x20fc812f473a5e69, 0x22294b092123a111, 0x6e51c1b6694583a6, 0x7e16b62b67ed3b13,
0x7ff6d68c194fed38, 0x7f7afed21eafa7b7, 0x7e8c81954c4b98a1, 0x7d560e164d1c7181,
0x20c21ffc1e252cc6, 0x6ec112c523dbd404, 0x79a43dbd46f19efd, 0x7cf2884c5dfac2b2,
0x70c3b91f3604fefc, 0x22749fdf27cd4a69, 0x7e3ccc576d442bd1, 0x207e5d523e47c33d
};
void test_newton_all() {
volatile uint64_t dummy = 0;
for (int i = 0; i < N; i++) {
dummy = sqrt_newton(input_data[i%180]);
}
}
void test_digit_all() {
volatile uint64_t dummy = 0;
for (int i = 0; i < N; i++) {
dummy = sqrt_digit(input_data[i%180]);
}
}
int main(void) {
test_newton_all();
test_digit_all();
return 0;
}
```
`input_data[0]` 輸出：
```
n_digit:32 //迭代數理論最高 32 次
n_newton:6 //迭代數理論最高 6 次
answer:3023032209,3023032209
```
查看各函式 CPU 佔用比例：
```
$ perf report
61.02% sqrt_test1 sqrt_test1 [.] test_digit_all
30.17% sqrt_test1 sqrt_test1 [.] test_newton_all
0.95% sqrt_test1 [unknown] [k] 0xffffffff8c047e73
0.81% sqrt_test1 [unknown] [k] 0xffffffff8d0015d0
0.35% sqrt_test1 ld-linux-x86-64.so.2 [.] __GI___tunables_init
0.33% sqrt_test1 [unknown] [k] 0xffffffff8bfd3d98
0.33% sqrt_test1 [unknown] [k] 0xffffffff8d001192
0.32% sqrt_test1 [unknown] [k] 0xffffffff8bda1937
0.32% sqrt_test1 [unknown] [k] 0xffffffff8be32e50
0.32% sqrt_test1 [unknown] [k] 0xffffffff8bd77d16
0.32% sqrt_test1 [unknown] [k] 0xffffffff8bd6ffb5
0.32% sqrt_test1 [unknown] [k] 0xffffffff8bd7f08c
0.32% sqrt_test1 [unknown] [k] 0xffffffff8bd918f9
0.32% sqrt_test1 [unknown] [k] 0xffffffff8bccdc52
0.32% sqrt_test1 [unknown] [k] 0xffffffff8bd7f084
0.32% sqrt_test1 [unknown] [k] 0xffffffff8d001168
0.32% sqrt_test1 [unknown] [k] 0xffffffff8be15065
0.32% sqrt_test1 [unknown] [k] 0xffffffff8d0011d9
0.32% sqrt_test1 [unknown] [k] 0xffffffff8d001217
0.32% sqrt_test1 [unknown] [k] 0xffffffff8be14ba7
0.32% sqrt_test1 [unknown] [k] 0xffffffff8cea9b2f
```
可知在雙方皆為 worst case 下，牛頓法的效率還是比逐位法要高。
## 那為何 Linux 核心不採用呢？
Linux 核心未採用牛頓法的主要原因，是因為某些處理器缺乏高效的除法指令，在這種情況下，牛頓法反而可能因需使用除法而導致效能下降。
Sign in with Wallet
Connect another wallet