# 開平方根運算

## 二分逼近法

當要找出某個數 $n$ 的平方根時，可以先定義上界和下界，當然 $n$ 是介於上下界的區間的，接著找出上下界的中間點，接著利用這個中間點去和 $n$ 比較，若是大於 $n$ 則代表 $n$ 介於下界和中間點的區間，此時將新的上界變為中間點，反之小於也一樣，下界變為中間點，並且遞迴下去，直到達成中止條件(和目標的差值小於某個數值等等)。

#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$ 的誤差小於某個數值)。

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` 做比對。


在執行時間上面, 兩者在前面遇到較小的數字時的所花費的時間明顯要比較長(待研究)

而我寫的牛頓法遇在遇到一些數字在收斂(找切線)的過程會非常非常久，仍須觀察問題出在哪裡。

再來觀察誤差，由於 floating number 的 mantissa 能夠表達足夠精準的位數大致為小數點後七位

> $\log_{10} {2^{24}}$≈7.224719

所以就擷取到後七位來做誤差觀察，從圖中可以發現兩條線是貼合的


先將測試資料改為 1 到 1000，再來將兩種方法的結果做誤差曲線

<!--  -->


可以看到在輸入值很小的時候，誤差非常大，推測是因為我本來的牛頓法寫的誤差容忍就不夠小，再來因為 IEEE 754 的 single precision float 只有 24 位元的精度，在小數區域，鄰近的兩個數字差距會變非常小，用 float 表示會喪失一定的精度，因為能表示的「有效位數」始終是 24 位元，當數值很小時，能精確表示的十進位數字變少，誤差比例相對就會變大。


## 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}$。

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;
}

接下來要把它擴充成可以對浮點數去做開根號

<!--

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include<linux/time.h>


float newton_sqrt(float target, float n)
{
    float tolerance = 0.000001f;

    while (fabs(target*target - n) > tolerance){
        float x = (target * target - n) / (2*target);
        target -= x;
    }

    return target;
}

int main(int argc, char *argv[]){
    struct timespec start, end;
    double cpu_time_used;
    
    float x = atof(argv[1]);
    clock_gettime(CLOCK_MONOTONIC, &start);
    float num = newton_sqrt(x,x);
    clock_gettime(CLOCK_MONOTONIC, &end);

    cpu_time_used = (end.tv_sec - start.tv_sec) + (end.tv_nsec - start.tv_nsec)/1e9;
    
    FILE *fp = fopen("result.txt", "a");
    if (fp != NULL){
        fprintf(fp, "%.12f %.12f\n" ,x ,cpu_time_used);
        fclose(fp);
    } else{
        printf("File opened failed\n");
    }

    return 0;

}

<!--
Max error: 2.14569999990033e-06
Mean error: 1.3244187690582708e-07
Max related error: 1.5172390462556032e-06
Mean related error: 2.6878791535642885e-08 -->