C語言牛頓迭代法求平方根：原理、初值選擇與快速逆平方根算法

題目回顧

迭代公式：X_{n+1} = (X_n + a/X_n) / 2，收斂條件 |X_{n+1} - X_n| < 1e-5。程式中初值 `x1 = a/2`，問為何選此初值？

牛頓法數學原理

求解 f(x) = x² - a = 0 的根。牛頓迭代公式：

x_{n+1} = x_n - f(x_n)/f'(x_n) = x_n - (x_n² - a)/(2x_n) = (x_n + a/x_n)/2

初值選擇分析

• 收斂條件：對凸函數 `f(x)=x²-a`，牛頓法在 **`x₀ > √a`** 時單調遞減收斂；在 `0 < x₀ < √a` 時第一步會跳到 `> √a` 後收斂。
• a/2 ≥ √a 當 a ≥ 4。對 `a ≥ 4`，`a/2` 滿足 `x₀ ≥ √a`，保證收斂。
• 對 `0 < a < 4`，`a/2 < √a`，但第一次迭代後會超過 `√a`，同樣收斂（只是多一步）。
• 通用性：`a/2` 簡單、無需條件判斷、對所有 `a > 0` 皆可收斂，是工程上常用的啟發式初值。

更嚴謹的初值選擇可用 `frexp` 取得指數部分估算，或查表法，但 `a/2` 在雙精度下通常 5-6 次迭代即達機器精度，足夠實用。

標準庫最佳實踐

**實務直接用 `` 的 `sqrt()`**：經過高度優化、硬體指令支援（`FSQRT`）、正確處理邊界（0、負數、次正規數、NaN、Inf）。

#include 
#include 

int main() {
    double a;
    printf("輸入非負數: ");
    scanf("%lf", &a);
    if (a < 0) { printf("輸入錯誤\n"); return 1; }
    printf("平方根 = %.10f\n", sqrt(a));
    return 0;
}

傳奇算法：快速逆平方根（Quake III 0x5f3759df）

原始回答提供的 `SquareRootFloat` 即為此算法變體。核心思想：

1. 利用 IEEE 754 單精度浮點數位元結構：`value = (1+m) * 2^(e-127)`
2. `log2(value) ≈ e + m - 127`（線性近似）
3. `log2(1/√value) = -0.5 * log2(value)`
4. 對應位元操作：`i = 0x5f3759df - (i >> 1)` 得到 `1/√x` 的極佳初估
5. 再執行 1-2 次牛頓法修正：`y = y * (1.5 - 0.5 * x * y * y)`

此算法在無硬體平方根指令的早期 GPU/CPU 上極快，現代 CPU 有 `RSQRTSS`/`SQRTSS` 指令，編譯器優化通常更優，但仍是位元操作經典教案。

完整可執行版本（含邊界處理）

#include 
#include 

float fast_inv_sqrt(float x) {
    if (x == 0) return INFINITY;
    float xhalf = 0.5f * x;
    int i = *(int*)&x;           // 位元重新解讀
    i = 0x5f3759df - (i >> 1);   // 魔術常數
    x = *(float*)&i;
    x = x * (1.5f - xhalf * x * x); // 牛頓法 1 次
    x = x * (1.5f - xhalf * x * x); // 可選第 2 次
    return x;
}

float fast_sqrt(float x) {
    if (x <= 0) return 0;
    return x * fast_inv_sqrt(x);
}

int main() {
    float a;
    printf("輸入非負數: ");
    scanf("%f", &a);
    printf("sqrt = %.8f\n", fast_sqrt(a));
    printf("std  = %.8f\n", sqrtf(a));
    return 0;
}