# 牛頓法
contributed by `TonyLinX`

## 定義
牛頓法是用來求根的方法，主要的目的是求近似值，所以通常的用法是給你一個函數 $f(x)$，要你求這個函數的根，它會透過不斷地迭代逐漸收斂到一個近似值，這個近似值會非常靠近你的根。

:::info
初始值 ($x_0$) 有幾種常見的選法：
* $x_0$ = a
* $x_0$ = a/2
:::

## 公式推導
假設給你一個函數

$$
f(x) = y = x^2 - a
$$

那這個 a 可以是任何數字，例如：2, 3, 5，所以這個函數就是我們常見的要你求 $\sqrt{2}$, $\sqrt{3}$, $\sqrt{5}$ 的近似值。但一開始我們不知道它的根是多少，只能透過牛頓迭代公式推導，牛頓迭代公式的邏輯：
- 先給定一個 $x_n$ (a)，若 $n=0$，就使用常見的初始值設定
- 在 $(x_n, f(x_n))$ 做切線，與 x 軸交於 $x_{n+1}$ (b) 點
- 在 $(x_{n+1}, f(x_{n+1}))$ 做切線，與 x 軸交於 $x_{n+2}$ (d) 點
- 一直重複上述操作，就會得到 $\sqrt{N}$

公式推導： 假設 $x_0$ = a , $f(x) = x^{2} - a$

求 $(x_0, f(x_0))$ 的切線
$$
m = f'(x_n) = \frac{f(x_n) - 0}{x_n - x_{n+1}}
=> x_n - x_{n+1} = \frac{f(x_n)}{f'(x_n)}
=> x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}
$$

通用的牛頓法迭代公式：
$$
x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}
$$

帶入 $f(x) = y = x^2 - a$ 以及 $f'(x) = 2x$ :
$$
x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}
=> x_{n+1} = x_n - \frac{x^2 - a}{2x}
=> x_{n+1} = \frac{1}{2}(x_n + \frac{a}{x_n})
$$

開平方根的牛頓法迭代公式：
$$
x_{n+1} = \frac{1}{2}(x_n + \frac{a}{x_n})
$$

#include <stdio.h>
#include <stdint.h>
#include <math.h>

float sqrt_w_fpu_newton(int a) {
    float x = a;  
    float epsilon = 1e-6f;
    
    while (fabsf(x * x - a) > epsilon)
        x = 0.5f * (x + a / x);

    return x;
}

int sqrt_no_fpu_newton(int a) {
    int x = a;  
    while (x * x > a  || (x + 1) * (x + 1) <= a)
        x = (x + a / x) / 2;
    
    return x;
}

float sqrt_fixed_q8_8(int a_q8_8) {
    
    int x = a_q8_8; //初始值
    
    for (int i = 0; i < 5; i++) {
        x = (x + (((int32_t)a_q8_8 <<8) / x)) >> 1;
    }
    
    return x /256.0f;  
}

int main() {
    int a[] = {1, 2, 3, 4, 5, 9};
    int len = (int)(sizeof(a) / sizeof(a[0]));
    
    for (int i = 0; i < len; i++) {
        float fpu_result = sqrt_w_fpu_newton(a[i]);
        int no_fpu_result = sqrt_no_fpu_newton(a[i]);  
        float fixed_result = sqrt_fixed_q8_8(a[i]<<8);
        
        float error = fabsf(fpu_result - (float)no_fpu_result)/fpu_result * 100;
        float error_fixed = fabsf(fpu_result - fixed_result)/fpu_result * 100;
        
        printf("sqrt(%d) ≈ %.6f (FPU), %d (no FPU), error = %12.8f %%, fixed = %.8f, error_fixed = %12.8f %%\n",
               a[i], fpu_result, no_fpu_result, error, fixed_result, error_fixed);
    }
    
    return 0;
}

結果：

sqrt(1) ≈ 1.000000 (FPU), 1 (no FPU), error =   0.00000000 %, fixed = 1.00000000, error_fixed =   0.00000000 %
sqrt(2) ≈ 1.414214 (FPU), 1 (no FPU), error =  29.28931999 %, fixed = 1.41406250, error_fixed =   0.01068001 %
sqrt(3) ≈ 1.732051 (FPU), 1 (no FPU), error =  42.26497650 %, fixed = 1.73046875, error_fixed =   0.09134522 %
sqrt(4) ≈ 2.000000 (FPU), 2 (no FPU), error =   0.00000000 %, fixed = 2.00000000, error_fixed =   0.00000000 %
sqrt(5) ≈ 2.236068 (FPU), 2 (no FPU), error =  10.55728149 %, fixed = 2.23437500, error_fixed =   0.07571372 %
sqrt(9) ≈ 3.000000 (FPU), 3 (no FPU), error =   0.00000000 %, fixed = 3.00000000, error_fixed =   0.00000000 %

## 開 k 次方根牛頓法公式推導

通用的牛頓法迭代公式：
$$
x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}
$$

套進我們的 $f(x) = x^k -a$ 得：

1. 函數:
$$
f(x) = x^k - a
$$
2. 導數:
$$
f'(x) = kx^{k-1}
$$

3. 帶入牛頓公式：
$$
x_{n+1} = x_n - \frac{x^{k}_n - a}{kx^{k-1}}
=> x_{n+1} = \frac{1}{k}((k-1)x_n + \frac{a}{x^{k-1}_n})
$$

開 k 次方根的牛頓法迭代公式：
$$
x_{n+1} = \frac{1}{k}((k-1)x_n + \frac{a}{x^{k-1}_n})
$$