開平方根運算

# 開平方根運算 ## 二分逼近法當要找出某個數 $n$ 的平方根時，可以先定義上界和下界，當然 $n$ 是介於上下界的區間的，接著找出上下界的中間點，接著利用這個中間點去和 $n$ 比較，若是大於 $n$ 則代表 $n$ 介於下界和中間點的區間，此時將新的上界變為中間點，反之小於也一樣，下界變為中間點，並且遞迴下去，直到達成中止條件(和目標的差值小於某個數值等等)。 ```c #include<stdio.h> #include<stdlib.h> float bisqrt(float n) { float top = n; float bottom = 1; float mid = (top+bottom)/2; while((top-bottom) >= 0.000001f){ if (mid * mid > n){ top = mid; } else if (mid * mid < n) { bottom = mid; } else return mid; mid = (top+bottom)/2; } return mid; } ``` ## 牛頓法用切線的 $0$ 去逼近非線性函數的 $0$ 根，也就是我們要求的根平方根， ### 固定點定理對某個函數 $f$，若存在一個點 $x$，使得 $f(x) = x$，那這個 $x$ 就叫做固定點。 #### 存在性證明要怎麼確認這個固定點保證存在呢? 先假設 $g(x)$ 是一個 $[a,b]$ 閉區間中連續的函數，並且 $g(x)\in[a,b]$ 對所有 $x\in[a,b]$，代表該函數值都落在其中。接著定義一個新函數 $h(x) = g(x) - x$，也是連續函數，這時就有兩種情況: 1. 第一種是固定點在邊界上，也就是 $g(a) = a; g(b) = b$ 2. 第二種就是固定點在區間內，所以 $g(a) - a > 0$ 而且 $g(b) - b < 0$ 而如果是第二種情況，這時候 $h(a) = g(a) - a > 0$ 且 $h(b) = g(b) - b < 0$，所以根據中間值定理: > 設函數 $f$ 在閉區間 $[a,b]$ 上連續，且 $f(a)\neq f(b)$，那麼對於任何介於 $f(a)$ 和 $f(b)$之間的值 $y$，一定存在某個 $c \in (a,b)$ 使得 $f(c) = y$。在 $a$ 到 $b$ 之間必定存在一個數 $p$ 能讓 $h(p) = 0$，而因為$h(p) = 0$，所以 $g(p) = p$，即存在固定點。 #### 唯一性證明利用均值定理: > 假設函數 $f$在閉區間 $[a,b]$ 連續且在開區間 $(a,b)$ 可微，則存在一點 $c$，$a < c < b$，使得 $f'(c) = \frac{f(b)-f(a)}{b-a}$> 但因為 $g(p) = p ; g(q) = q$，所以上面式子變成 $g'(c) = \frac{p-q}{p-q}=1$，違反了假設 $\vert g'(x) \vert\leq M <1$，故 $p \neq q$ 不成立，固定點唯一。 #### 近似求解牛頓法的原理，是透過在函數圖形上某一點作切線，並以該切線與 $x$軸的交點作為下一次的估計值，如此反覆進行，以逼近方程式實根的近似解。在求平方根的情境中，我們考慮的方程是： $f(x) = x^2 - N = 0$ 一開始我們給予初始值 $a_0$，計算函數值 $f(a_0)$，並透過倒數 $f'(x) = 2x$ 得出在該點的切線斜率。該切線與$x$軸的焦點，根據牛頓法公式如下: $$x_{n+1}=x_n - \frac{f(x_n)}{f'(x_n)} = x_n - \frac{(x_n)^2 - N}{2x_n}=\frac12(x_n + \frac{N}{x_n})$$ 再來取切線和$x$軸的交點作為下一次的輸入值，以此遞迴下去，直到達成設定條件(例如和 $N$ 的誤差小於某個數值)。 ```c float nt_sqrt(float target, float n) { float tolerance = 0.0001f; while (fabs(target*target - n) > tolerance){ float x = (target * target - n) / (2*target); target -= x; } return target; } ``` 以此演算法下去做實驗，觀察一下執行時間，測資從 $1$ 到 $100$，step 為 $0.01$，並且跟 `sqrtf` 做比對。 ![result](https://hackmd.io/_uploads/B1_lAFGVgl.png) 在執行時間上面, 兩者在前面遇到較小的數字時的所花費的時間明顯要比較長(待研究) 而我寫的牛頓法遇在遇到一些數字在收斂(找切線)的過程會非常非常久，仍須觀察問題出在哪裡。再來觀察誤差，由於 floating number 的 mantissa 能夠表達足夠精準的位數大致為小數點後七位 > $\log_{10} {2^{24}}$≈7.224719 所以就擷取到後七位來做誤差觀察，從圖中可以發現兩條線是貼合的 ![result](https://hackmd.io/_uploads/rk8wvAGNgl.png) 先將測試資料改為 1 到 1000，再來將兩種方法的結果做誤差曲線  ![result](https://hackmd.io/_uploads/SkcKiZENgg.png) 可以看到在輸入值很小的時候，誤差非常大，推測是因為我本來的牛頓法寫的誤差容忍就不夠小，再來因為 IEEE 754 的 single precision float 只有 24 位元的精度，在小數區域，鄰近的兩個數字差距會變非常小，用 float 表示會喪失一定的精度，因為能表示的「有效位數」始終是 24 位元，當數值很小時，能精確表示的十進位數字變少，誤差比例相對就會變大。 ![image](https://hackmd.io/_uploads/HJdTEzN4eg.png) ## Digit-by-digit (base 2) $N^2 = (a_n+.....+a_0)$，其中每個 $a_m$ 因為是位元的關係在二進位中僅會是$0$或$1$的表示，在十進位中就是 $2^m$ 或是 $0$ $a_m = \begin{cases} 2^m \\ 0 \end{cases}$，而目標就是去迭代這些位元，從 $n$ 到 $0$，我們可以先建立近似解 $$P_m = a_n + a_{n-1} + .... + a_{m+1} + a_m$$ 而為了要確定 $a_m$ 究竟是 $2^m$ 或是 $0$，我們使 $P_m = P_{m+1} + a_m$，如果 ${P_m}^2 \leq N^2$，那就代表說我們近似解的平方包括了 $2^m$ 也不會超過目標值，所以 $a_m = 2^m$，否則 $a_m = 0$，$P_m = P_{m+1}$。 ```c int digit(int n) { unsigned int bit = 1<<31; unsigned int res = 0; while (!(bit&n)){ bit >>= 1; } int ans = 0; for (;bit > 0;bit >>= 1){ int tmp = bit | ans; if (tmp * tmp <= n){ ans = tmp; } } return ans; } ``` 接下來要把它擴充成可以對浮點數去做開根號