---
tags: sysprog
---

# 從 √2 的存在談開平方根的快速運算
> 貢獻者: jouae, jserv

## $\sqrt{2}$ 的緣起

邊長為 1 的正方形對角線長度是多少？中學數學提到，直角邊長為 1 的等腰**直角** 三角形的斜邊長就是 $\sqrt{2}$，而這個數是個無理數，也就是說它無法寫成二個整數相除的形式。
> 延伸閱讀: [第一次數學危機，天才是怎麼被扼殺的？](https://youtu.be/nAOVQEcqjSM)

畢氏定理的畢氏，指古希臘學者[畢達哥拉斯](https://en.wikipedia.org/wiki/Pythagoras)。相傳此定理在當時是由畢達哥拉斯所發現，雖然古巴比倫人在他出生前的一千多年前就已發現該道理，但至少可知道，畢達哥拉斯把這個定理又帶回來給西方世界，所以西方世界多用畢氏定理來稱這個常見的定理。

在數學之外，畢達哥拉斯也是集數學、天文、音樂理論於一身的高手，他願意將部分的知識傳播給大眾，而且他也是當時少數讓女性也能來聽講的學者 (當時禁止女性聽公開演講)。不過，這些知識也有分等級，分享給大眾聽的，是比較一般的知識，若你想提問，或者你想知道更高深的知識，你需要加入畢達哥拉斯學派。

畢達哥拉斯學派的正統信徒必須宣誓，不會將學派的秘密以及較高深的學說外流出去，同時也必須經過嚴格的訓練，才會被承認為正式的教徒。除此之外，還必須嚴守各種的清規戒律，某種程度就是在修行；同時也會將自己所研究出來的新成果，成就歸於導師畢達哥拉斯。

畢達哥拉斯等人對數學是有一定的狂熱，他們不是為了實際的應用而去鑽研數學。他們相信，透過數學，即可以探究這世界，上段所提的音樂理論，也是從數學推演出來的。為了抽象地去形容宇宙，或者單純只是為了美，畢達哥拉斯的學派在數學上也有很多的特殊想法，例如：萬物皆數。

萬物皆數，照字面上的意思就是萬物的本質皆是數字，這也是為什麼他們相信能利用數學來尋找真理。這個數是有嚴謹定義的，或者說，他們認為所有的數都有一個性質：[可公度性](https://en.wikipedia.org/wiki/Commensurability_(mathematics))。

要如何理解可公度性呢？若用線段長度來看，兩條線段之間必定能用一個特殊長度的整數倍來形容，例如 2 公分和 6 公分，可寫成 $1 \times 2$ 公分和 $3 \times 2$ 公分，皆可用 2 公分的整數倍來形容；7 公分和 5.5 公分，可寫成 $14 \times \dfrac{1}{2}$ 公分和 $11 \times \dfrac{1}{2}$ 公分。要找這個特殊長度，可以先把長的扣掉短的，得到的新線段，再拿第二短的去扣掉最短的又得到新線段，接著重複上述過程，直到能剛好減到 0，此時最短線段的長度就是這個特殊長度。

例如 7, 5.5 $\to$ 7, 5.5, 1.5 $\to$ 7, 5.5, 4, 1.5 $\to$ 7, 5.5, 4, 2.5, 1.5 $\to$
　7, 5.5, 4, 2.5, 1.5, 1 $\to$ 7, 5.5, 4, 2.5, 1.5 ,1 , 0.5，最後得到 7, 5.5, 4, 2.5, 1.5,1, 0.5, 0.5，這個特殊長度即是 `0.5`，類似輾轉相除法找最大公因數的手法。

可公度性的另一種說法是：畢達哥拉斯認為所有的數，皆是整數或者二個非零整數的比值，也就是有理數。等等，這觀點和我們今日所知的數學原理不一致，發展的轉折就要提及[希帕索斯](https://en.wikipedia.org/wiki/Hippasus) (Hippasus)。

希帕索斯是畢達哥拉斯學派的信徒，傳聞他非常認真，亦愛鑽研可公度性這個性質。有次他鑽研正五邊形的性質時，發現有兩條線段的長度，不論怎樣也無法找到可公度性所需要的特殊長度，只能無窮無盡的切割下去，沒有停止的那天。換言之，這代表這兩條線段並沒有可公度性這個性質，那麼，他們當初所相信的一個性質 —— 所有數都是可公度的，便是錯誤的，或者說，這世界上存在著不是整數、亦不能寫成兩整數比值的數。

後來希帕索斯改良使用更為直覺的方式，利用正方形的對角線與邊長，來說明存在兩數是不可公度的。
> 延伸閱讀：〈[不可公度性的發現與克服——Hippasus 和 Eudoxus 對于人類理性文明的重大貢獻](http://episte.math.ntu.edu.tw/articles/ar/ar_wy_geo_02/page4.html)〉

在畢達哥拉斯的指使或者漠視下，希帕索斯被門徒淹死在河裡。有人說畢達哥拉斯也搞政治，但其地位後來已岌岌可危，不容許再出現「學說錯誤」這樣的汙點。有的故事則說，希帕索斯提出不可公度這件事時，畢達哥拉斯已死去一段日子，被處死這件事情與畢達哥拉斯沒有太大關係。也有故事說，希帕索斯是與同門搭船時，與大家爭論不可公度這件事情時，被惱怒的同門丟下船溺死的。

無論是哪個版本，都指向諷刺的結果：希帕索斯運用理性找出的真相，被當時最能代表理性、邏輯和數學的一群人否定，而他的生命也被這群人不理性的抹殺。

仔細想想，畢達哥拉斯或者其學派的人，不可能不理解希帕索斯提出的疑慮，或者說不定，以畢達哥拉斯這樣的大師來說，他早就發現這件事情也不稀奇，但何以這些信徒仍選擇抹殺這個發現問題的學生呢？

因為信仰，不是宗教上的信仰，而是內心認同的價值觀。

畢達哥拉斯一生都在追求真理，同時他也不斷的宣揚萬物皆數這樣的想法。試著這麼想，如果有一天，你發現「真理」與你所想像的不同，跟你這一生所信奉的不一樣，也與你這大半輩子所宣揚的有所違背，你會做甚麼選擇？承認錯誤，讓自己變得更接近真理？

只可惜多數的人，可能也包含畢達哥拉斯與其信徒，沒有那麼大的勇氣去否定他們的信仰，去承認過去的努力只是個錯誤。因為這也是在否定他自己。如果他們是自己發現了這個事實，那麼能採取的，只能欺騙自己。

如果一件事情變成信仰，而信仰者卻沒有透過自省或質疑來去強化信念，那麼他的信仰即使看起來很堅定，也就只是虛有其表。在這種情況之下，你會發現他們為鞏固那並不是很堅固的信仰，會編撰出很多理由，即使這些理由從邏輯或知識上來看根本過於荒謬。

那麼，是否因他們邏輯推理能力不足、學識不夠豐富，才無法藉由理性與科學證據來修正自己的信仰呢？

不是，邏輯和知識當然是承認錯誤的必要條件，但即使他們有這些能力，若不曾自省與質疑，那麼這些能力也只會變成他們用來找尋藉口或荒謬理由的工具而已。看前面的故事就知道，即使當代最有邏輯與學識的一群人，一旦把教條當成不可質疑的信仰時，就算有人提出毫無瑕疵的證據質疑其說法，也只會像希帕索斯一樣，被這樣想要堅定自己搖搖欲墜的信仰之人推向深淵抹除，而在這種情形下，恐怕還有很多不同名字的希帕索斯這樣的消失在這個世界上。

我們藉由中學的數學，已知無理數的存在 (或者不可公度的存在)，可大概推得，希帕索斯的發現並沒有隨著他的生命而消失。這個理論後來似乎隨著畢達哥拉斯學派的式微或者教條約束力的減弱，也漸漸流傳出去。這理論沒有沉到水裡，至於希帕索斯的發現是否對他的信仰 (數學) 造成了傷害呢？的確這個發現造成了第一次數學危機，但是也讓數學發展向前邁進了一大步。

別讓信仰殺死希帕索斯，因為他會帶領我們更走向真理。

> 以上改寫自 [因信仰而死的希帕索斯](https://m.facebook.com/notes/%E5%AE%8B%E6%94%BF%E8%92%B2/%E5%9B%A0%E4%BF%A1%E4%BB%B0%E8%80%8C%E6%AD%BB%E7%9A%84%E5%B8%8C%E5%B8%95%E7%B4%A2%E6%96%AF/10207331849244095/)

## 兩河流域古文明時期

針對這個數學問題，考古學家發現在西元前 1800 年的古巴比倫泥板上已記載答案。

兩河流域文明時代最早的居民是蘇美人，在公元前 4000 年就抵達幼發拉底河和底格里斯河之間的美索不達米亞平原，大體位於現今的伊拉克。屬于塞姆語系的阿卡德人、巴比倫人、亞述人以及迦勒底人，繼承和發展蘇美人的成就，使兩河流域的文明成為人類文明史上重要的一頁。其中，巴比倫人的成就最大，因此兩河流域的文明又被稱為巴比倫文明。

兩河流域的居民主要用使用牛、驢拉著木犁耕地，古代兩河地區的金屬制造工藝達到了相當純熟的水平。古巴比倫人使用楔形文字表示數字，以[六十進位](https://en.wikipedia.org/wiki/Sexagesimal)進行計算。考古學家在中東地區發掘了數以萬計的古巴比倫泥板碎片，其中一片現存於耶魯大學的泥板依稀刻印著一個正方形的對角線的圖案，以及一排以楔形文字的數字。

後來考古學家解譯出上圖這組楔形數字是: `1`, `24`, `51`, `10`。
由於古巴比倫人用的是六十位系統，將這組數字轉換為十進位的數值得到以下：
$$
1 + \dfrac{24}{60} + \dfrac{51}{60^2} + \dfrac{10}{60^3} = \dfrac{305470}{216000} = 1.41421296...
$$

這個數字極為近似 $\sqrt{2}$ 的值 $1.41421356$。

泥板上還給出了一個計算的例子，在正方形邊線上標註數字 $30$，也就是十進位 $0.5$，代表正方形邊長。這個時候對角線長度的楔形數字為 $42$, $25$, $35$，轉換十進位得到：
$$
0 + \dfrac{42}{60} + \dfrac{25}{60^2} + \dfrac{35}{60^3} = \dfrac{152735}{216000} = 0.707106481...
$$

這個值正是 $\sqrt{2}$ 的一半，也就是邊長 $0.5$ 正方形的對角線長度。

至於東亞採用一日十二時辰的制度可能源於漢代曆法，採取十二支作為劃分是基於天文觀測的結果。後來為了更精確的計時，改用一刻六十分的制度，考量點如下：
* 由十時制改為十二時制,是為了更精確地反映一天中不同時段的差異。觀察天文現象發現，一天的時間不能平均分為直觀的十份，十二等分更能貼近實際情況。
* 採用六十分作為一刻的進位制，是因地球繞太陽公轉的周期約為 365.25天，與 360 度的圓周相吻合。六十進位也更利於分秒的計算。

整個時間計量系統的發展反映人類持續強化對自然界週期性運行規律的認識。從十時制到十二時制、從六十分制到公制，都是為了更適應天文觀測和農耕生活的需要。

### 十分逼近法算 $\sqrt{2}$

古巴比倫人的時代距今已約 3800 年之久，我們無法確實知道古巴比倫人是如何計算出如此精準的 $\sqrt{2}$ 數值。

在估算 $\sqrt{2}$ 的數值前，先回想平方根的定義：一個數 n 的平方根會滿足
$$
n = (\sqrt n)^2
$$

現在考慮三個數 `1`, `2`, `4`，由於 `1 < 2 < 4`，如果對這三個數開平方根也會有下面的關係：
$$
\sqrt{1} < \sqrt{2} < \sqrt{4}
$$

這裡 `1` 和 `4` 都是平方數，可以很容易地算出 $\sqrt{1}$ = 1 和 $\sqrt{4}$ = 2，帶入上式可得到：
$$
1 < \sqrt{2} < 2
$$

這個式子告訴我們 $\sqrt{2}$ 的數值介於正整數 `1` 和 `2` 之間，不過我們還不能確定 $\sqrt{2}$ 的數值較靠近 `1` 抑或靠近 `2`。

這個情況下我們可以用一個方法更進一步確定 $\sqrt{2}$ 的數值。將 `1` 到 `2` 再細分為 10 等分，也就是 `1.1`, `1.2`, …, `1.9`, `2`，再計算每一個細分刻度的平方:

在數線上方我們列出 n 平方根的十等分值，在下方列出相對應的 n = (n平方根)^2^ 的值。我們可以看到 n = (n平方根)^2^ 有一個重要的變化點，也就是由 `1.96` 遞增到 `2.25` 之處，遞增的 n 值由小於 `2` 變成大於 `2`：

而我們要找的 $\sqrt{2}$ 答案應該就在這個區間。

藉由這樣細分刻度的方法，我們可以將 $\sqrt{2}$ 的答案逼近到小數點後第一位，知道了 $\sqrt{2}$ 的值介於 `1.4` 和 `1.5` 之間。這個方法就稱為十分逼近法。

再一次利用十分逼近法，我們將 `1.4` 和 `1.5` 之間再細分為 10 等分：

這一次跨越 `2` 的區間落在 `1.9881` 和 `2.0164` 之間：

在這個關鍵的區間我們可以再次觀察到：
$$
1.9881 < 2 < 2.0164 \\
\sqrt{1.9881} < \sqrt{2} < \sqrt{2.0164} \\
1.41 < \sqrt{2} < 1.42
$$

因此確定了 $\sqrt{2}$ 在小數點後第 2 位的數值，$\sqrt{2}$ 的值會大於 `1.41`，而小於 `1.42`。如此重複運用十分逼近法可以求得更精確的 $\sqrt{2}$ 數值。

### Google 內建的計算機

根據平方根的定義，一個數 n 的平方根會滿足 n = (n平方根)^2^。如果只考慮 n 平方根是正整數的部分，n 的平方根 = $\sqrt{n}$。

平方根的英文叫做 square root，在數學表示式裡常常會用縮寫 `sqrt` 表示平方根，$\sqrt{n}$ 也可以記成 n 平方根 = `sqrt(n)`。

Google 的網頁搜尋列支援計算機功能，我們可在 Google 搜尋列鍵入`sqrt(2)`:

這時 Google 會顯示出一個計算機的畫面，告訴我們 2 的平方根等於`1.41421356237`:

這個看似簡單的動作事實上包含了輸入輸出的步驟和一個內建的平方根演算法；我們一開始在計算機裡輸入 n 的值等於 2，計算機依據內建的平方根演算法幫我們算出 n 平方根 = `1.41421356237`，然後輸出在答案列。

Google 也提供函式繪圖的功能，我們試試 Google 搜尋列中鍵入`y=sqrt(x)`:

我們會看到 Google 返回一張 `sqrt(x)` 的圖:

這張圖的意義在於將函數 `y=sqrt(x)` 的幾何圖形畫在二維的卡氏座標系裡。可見縱軸座標表示的就是 n 平方根，橫軸座標表示的是 n。

`y=sqrt(x)` 的圖表畫出了一條抛物線，每一個 `x` 的點都可以透過平方根函式求得一個對應的 `y` 值，因此每一個橫軸上的 `x` 座標，對應到縱軸的y座標，就可以畫出一個點。例如圖中的 `x = 2.00000347`，平方根函式就會算出對應的 `y = 1.41421479`，然後在座標系描出這個點。

當程式掃描橫軸上很多的 `x` 座標，就可以繪製出這條 `y=sqrt(x)` 的抛物線。

### 二分逼近法求 $\sqrt{2}$

`y = sqrt(x)`，或說 $y = \sqrt{x}$ 是一條抛物線狀的連續函數，要在這種連續函數上找到 $\sqrt{2}$ 的值，[二分逼近法](https://en.wikipedia.org/wiki/Bisection_method)是種簡單而可靠的計算方法。

二分逼近法有 4 個基本步驟：
1. 先設定左右二個邊界值。先前已經知道 $\sqrt{2}$ 的值介於 1 和 2 之間，我們將第一個邊界值設在 `y = 1` 的地方，再將第二個邊界值設在 `y = 2` 的地方。
2. 再來以二個邊界點平均求中間點，在 1 到 2 這個範圍中間點是`y = 1.5`，我們將中間點也標註在線段上。
3. 接下來檢查`x` 值在落在左右那一個區間內。由於 $\sqrt{2}$^2^ = 2 介於 1 和 2.25，所以我們可以知道 $\sqrt{2}$ 的值會介於 1 到 `1.5`。
4. 重複步驟 1–3，直到 $\sqrt{2}$ 的值收斂到理想範圍。








用二分逼近法重複到第 8 次可以計算出 $\sqrt{2}$ 的值到小數點後第二位，實際的值會介於 `1.4140625` 和 `1.41796` 之間，不難發現二分計算法收斂的速度不快。

## 導入遞迴求解

遞迴是程式語言中的一個重要概念。用遞迴法解決問題首先需要找出問題的規律性，然後將問題分解為同類的子問題，重複求解直到確定滿足我們要的答案。

遞迴強大的地方在於允許使用者用有限的語句描述無限的物件。許多電腦程式語言支援遞迴法，允許一個函式呼叫它自己，來達到遞迴的目的。

在二分逼近法和十分逼近法的程式裡，我們都運用遞迴技巧。遞迴有三個重要的步驟：
1. 設定問題的初始值。在二分設定法中，我們製做了二個變數，上限值和下限值。上限值的初始設定為使用者輸入的 n 值，下限值的初始設定為 1。
2. 將一個問題分解為可重複解決的較小的子問題。計算出上、下限值的平均值當做一個 $\sqrt{n}$ 的檢查點，$\sqrt{n}$ 的檢查點會將數線分為左右兩半，檢查 n 是落在左段或右段。如果 n 落在左段，則重新設定上限值為 $\sqrt{n}$ 的檢查點。如果 n 落在右段，則重新設定下限值為 $\sqrt{n}$ 的檢查點。重複上面的二個步驟。如下圖所示，$\sqrt{n}$ 的檢查點的漸漸逼近真正的 $\sqrt{n}$。
3. 設定找到答案的條件。如果條件尚未滿足，進行下一次的遞迴呼叫，一旦滿足這個條件則離開遞迴呼叫。

如果 $\sqrt{n}$ 檢查點的平方和 n 的絶對值差距在 `0.0000001` 之內，$\sqrt{n}$ 已經有相當的精準度，則跳出遞迴呼叫。

展示影片: [二分逼近法和十分逼近法求平方根](https://youtu.be/PHPj0jLhcnM)

## 牛頓法求平方根
以下數學理論部分參考：
1. 臺灣師範大學數學系黃聰明教授的[數值分析](https://sites.google.com/gapps.ntnu.edu.tw/tmhuang/courses/numerical-analysis?authuser=0)
2. 〈[理論與計算的結合：迭代法](https://web.math.sinica.edu.tw/mathmedia/HTMLarticle18.jsp?mID=44203)〉，作者：林琦焜教授，發表於數學傳播第 44 卷第 2 期
3. 《[Numerical Analysis 10/e](https://www.amazon.com/Numerical-Analysis-Richard-L-Burden/dp/1305253663)》 第二章
4. [COMPUTING AND ESTIMATING THE RATE OF CONVERGENCE](https://www.math-cs.gordon.edu/courses/ma342/handouts/rate.pdf)

實作部分參考：
1. [Fundamentals of Numerical Computation: 4.2. Fixed-point iteration](https://tobydriscoll.net/fnc-julia/nonlineqn/fixed-point.html)
2. 《[Numerical Analysis 10/e](https://www.amazon.com/Numerical-Analysis-Richard-L-Burden/dp/1305253663)》 第二章

### 固定點定理
定義：如果一點 $p\in I$ 為單變數函數 $g:I\rightarrow I$ 的固定點，則 $g(p)=p$。

通常在數學中對函數、方程等等，我們會對其(1)解的存在性和(2)解的唯一性感興趣，例如線性方程組中從矩陣的零空間（null space）知道線性系統的解空間，解是否唯一（nullity is zero）或是多組解（nullity is nonzero, hence the solutions can be the linear combination of the basis in null space）。

所以對於一個單變數函數在一區間是否有固定點？如果有，那這固定點是否唯一？是我們關心的。

- [ ] 存在性

假設函數 $g$ 在一閉區間 $I=[a,b]$ 連續且 $g(x)\in[a,b]$ 對所有 $x\in[a,b]$，則函數 $g$ 最少有一固定點 $p\in[a,b]$。

存在性證明：
令函數 $g$ 在一閉區間 $I=[a,b]$ 連續且 $a\leq g(x)\leq b$ 對所有 $x\in[a,b]$ 分二個情況討論。
1. 固定點在區間邊界 $g(a)=a$ 或是 $g(b)=b$
2. 固定點在區間內部（interior），$g(a) > a$ 且 $g(b) < b$。

考慮一輔助函數 $h(x)=g(x)-x$ ，由於連續函數 $g$ 減連續函數 $x$ 還是連續函數，$h$ 在區間 $I$ 為連續函數，根據情況二的假設此連續函數
$$
h(a)=g(a)-a > 0 \quad \text{and}\quad h(b)=g(b)-b < 0 .
$$

看到此處有幾個關鍵點：
1. 函數 $h$ 在閉區間 $I$ 連續，
2. $h(a)h(b) < 0$

這就是在中學學過勘根定理的充分條件，勘根定理又稱 [Bolzano's theorem](https://en.wikipedia.org/wiki/Bolzano%E2%80%93Weierstrass_theorem) 為中間值定理的推廣。所以根據中間值定理，存在一點 $p\in[a,b]$ 使得連續函數 $h(p)=g(p)-p=0$，得 $p$ 為 $g$ 在閉區間內部的固定點。根據情況一與情況二，$g$ 在閉區間 $I$ 至少存在一固定點 $p$。

- [ ] 唯一性

承上述存在性，連續函數 $g$ 且 $g(x)\in[a,b]$ 在開區間 $(a,b)$ 一次可微，並且存在一正數 $M<1$ 使得 $$\vert g'(x) \vert\leq M, \quad \text{for all }x\in(a,b)$$ 則存在唯一固定點 $p\in[a,b]$。

唯一性證明：
令連續函數 $g$ 在開區間 $(a,b)$ 一次可微，與正數 $M<1$ 且 $\vert g'(x) \vert\leq M$ 對所有 $x\in(a,b)$ 。
**證明唯一性的常用手法之一：假設其不唯一，使用反證法得證其唯一。**

假設一次可微函數 $g$ 有兩固定點 $p$ 跟 $q$ 且 $p\neq q$，則 $g(p)=p$ 和 $g(q)=q$ 。從關鍵點(1)開區間函數可微，(2)閉區間函數連續。則對兩函數值，使用均值定理（mean value theorem）： 存在一點 $c\in(a,b)$ 使得
$$
g'(c)=\dfrac{g(p)-g(q)}{p-q}=\dfrac{p-q}{p-q}=1
$$

從 $\vert g'(c) \vert = 1$，得到跟一開始命題 $\vert g'(x) \vert\leq M <1$ 對所有 $x\in(a,b)$ 矛盾。故假設兩固定點 $p\neq q$ 不成立，證畢得 $p=q$ 固定點唯一。

結合上述存在性與唯一性便是**固定點定理（fixed point theorem）**。

- [ ] 固定點定理

令函數 $g$ 連續和 $g(x)\in [a,b]$。且 $g'$ 在開區間 $(a,b)$ 可微分及一常數 $0<M<1$ 使得
$$
\vert g'(x) \vert \leq M, \quad \text{ for all } x\in(a,b).
$$

藉由迭代法
\begin{align}
&p_0\in[a,b], \\
&p_n=g(p_{n-1}), \quad n\geq1.
\end{align}

得到的數列 $\lbrace p_i\rbrace$ 會收斂至唯一的固定點 $p\in[a,b]$。

固定點定理證明：
從存在性與唯一性得知，唯一固定點 $p$ 存在於 $[a,b]$ ，又 $g$ 的定義域與值域相同，藉由迭代法確實可以構造數列 $\lbrace p_i \rbrace^{\infty}_{i=0}$ 其中 $p_i\in[a,b]$ 對所有 $i$ 。

在證明收斂的部分，我們會想看 $\lbrace p_i \rbrace^{\infty}_{i=0}$ 數列中，當 $i$ 趨近無窮時 $p_i$ 跟 $p$ 的距離是否足夠靠近？如果是則說明數列最終收斂至 $p$ 。從 $p_i$ 跟 $p$ 距離得
$$
\lvert p_i-p\rvert=\lvert g(p_{i-1})-p\rvert,
$$

根據均值定理，存在 $\xi_{i}\in(a,b)$ 使得
$$
\lvert g(p_{i-1})-p\rvert = \lvert g'(\xi_{i}) \rvert \lvert p_{i-1}-p\rvert,
$$

從 $\lvert g'(\xi_{i}) \rvert\leq M$ 可以得到
$$
\lvert g'(\xi_{i}) \rvert \lvert p_{i-1}-p\rvert \leq M \lvert p_{i-1}-p \rvert.
$$

將上述推論的過程再做一次會得
$$
\lvert p_i-p\rvert \leq M \lvert p_{i-1}-p \rvert \leq M^2 \lvert p_{i-2}-p \rvert.
$$

以此類推最終可以得到一不等式
$$
\lvert p_i-p\rvert \leq M^{i} \lvert p_{0}-p \rvert.
$$

根據 $0<M<1$ ， 極限 $\displaystyle\lim_{i \to \infty} M^k=0$ 故
$$
\displaystyle\lim_{i \to \infty} \lvert p_i-p\rvert \leq M^k \lim_{i \to \infty}\lvert p_{0}-p \rvert =0.
$$

最後得證數列 $\lbrace p_i \rbrace^{\infty}_{i=0}$ 收斂至 $p$。

#### 例題 1
> 出自《[Numerical Analysis 10/e](https://www.amazon.com/Numerical-Analysis-Richard-L-Burden/dp/1305253663)》第 60 頁
> [程式碼](https://github.com/jouae/sqrt/blob/main/fixed_point.c)

考慮函數 $f(x)= x^3+4x^2-10$，使用 fixed-point iteration 在區間 $[1,2]$ 求根
$$
f(x) = x^3+4x^2-10 = 0
$$
數值解為何？

從固定點的定義，我們需要定義一個函數 $g(p)=p$ ，所以藉由原函數 $f(x)$ 可以構造出 $g(x)=x-f(x)$ ，當 $p$ 為固定點時 $g(p)=p-f(p)=p$ ，則 $f(p)=0$ ，固定點 $p$ 則為所求 $f$ 之根。書中給出的參考解為
$$
1.365230013
$$

則固定點迭代函數 $g_1$ 表示為
$$
g_1(x) = x - x^3-4x^2+10
$$

在單精度浮點數 `float` 運算下， `p_0` 為初始值，固定設定 $1.5$ ，其中相對誤差 (relative error) 為
$$
\text{relative error}=\dfrac{\lvert p_{exact} - p_i \rvert}{p_{exact}}
$$

使用這個方式來跑固定點迭代如下：

p_0: 1.500000
1 iteration:
    p_i: -0.875000
    relative error: 2.240230
2 iteration:
    p_i: 6.732422
    relative error: 5.367192
3 iteration:
    p_i: -469.720032
    relative error: 471.085266
4 iteration:
    p_i: 102754568.000000
    relative error: 102754568.000000
5 iteration:
    p_i: -1084934299940933513248768.000000
    relative error: 1084934299940933513248768.000000
6 iteration:
    p_i: -nan
    relative error: nan

> 從上圖看橘線部分為 $y=x$ 藍線部分為 $g(x)$，其實整個迭代過程圖像化就如圖中箭頭指向的過程。

何以程式輸出 `nan` 呢？固定點迭代不該會收斂嗎？

固定點的存在根據上述推導
1. 函數要連續
2. 函數要為 [contraction mapping](https://en.wikipedia.org/wiki/Contraction_mapping) ，白話講有一縮放係數 $0<M\leq1$ 使得 $\lvert\lvert f(x)-f(y) \rvert\rvert \leq M \lvert\lvert x-y \rvert\rvert$ 對 $x,y\in[1,2]$ 。

那來檢查下剛剛的 $g(x)=x-f(x)$ ，確實函數為連續的， $f(x)$ 和 $x$ 為多項式函數故連續函數減連續函數還是連續。那定義域與值域相同嗎？檢查 $g$ 一次微分
$$
g'(x) = -3x^2-8x-1<0, \quad \text{for }x\in[1,2].
$$
可以看出函數 $g$ 在區間 $[1,2]$ 為嚴格遞減，則由 $g(2)=-12$ 與 $g(1)=6$，得到值域為 $[-12,6]$ ，不滿足固定點的存在性。

所以這題無法使用固定點迭代求根了嗎？答案是可以的。但要改寫原函數。

從另一個角度來想，**原函數是因為非 contraction mapping 而導致固定點不存在，所以找一個有相同根並滿足存在性的函數即可**。

從原方程改寫
$$
x^3+4x^2-10=0
$$
將二次項以外移至等號右側後開根號
$$
4x^2 = 10-x^3
$$
由於要找的根在 $[1,2]$ 取正號部分即可
$$
x = \dfrac{1}{2}\sqrt{10-x^3}
$$
稍微改寫下形式，得到函數 $g_2$
$$
g(x)=x - \dfrac{1}{2}\sqrt{10-x^3} = x - g_2(x)
$$

再來確認函數 $g_2$ 的固定點存不存在：
$g_2$ 在區間 $[1,2]$ 連續，從一次微分判斷遞增遞減
$$
g_2'(x) = -\dfrac{3}{4}x^2(10-x^3)^{-\frac{1}{2}} < 0, \quad \text{for } x\in[1,2].
$$
得到 $g_2$ 在區間 $[1,2]$ 嚴格遞減，再來確認兩端點是否落在同區間內
$$
g_2(1) = \dfrac{3}{2}= 1.5,\quad g_2(2) = \dfrac{\sqrt{2}}{2}\approx 0.7071.
$$
得知值域 $[1.5, \frac{\sqrt{2}}{2}]$ ，固定點存在。

但此函數不滿足固定點定理的唯一性
$$
\lvert g_2'(x) \rvert\leq \frac{3\sqrt{2}}{2}\neq1, \quad \text{for }x\in[1,2].
$$

這種不滿足固定點唯一性，但固定點存在的函數，是不穩定 (unstable) 的。

同時還有兩個函數可以求解原函數之根：
$$
g_3(x)=\left(\dfrac{10}{4+x}\right)^{\frac{1}{2}},\quad g_4(x)=x-\dfrac{x^3+4x^2-10}{3x^2+8x}.
$$

其中一次微分為：
$$
\begin{align}
g'_3(x)&=-\dfrac{\sqrt{10}}{2}(4+x)^{-\frac{3}{2}}, \\
g'_4(x)&=\dfrac{(6x+8)(x^3+4x^2-10)}{(3x^2+8x)^2}.
\end{align}
$$

在數值方法的分析中關心兩件事情：
1. 收斂速率
2. 誤差分析

以下會介紹關於收斂速率的計算以及固定點誤差分析方式。

#### 收斂速率
使用 $g_1$ 到 $g_4$ 的固定點迭代，並計算其相對誤差如下：

relative errors:
iteration\function
                   g1              g2              g3              g4
 1-iter: 1.640918e+00    5.733551e-02    1.232774e-02    5.935535e-03
 2-iter: 3.931346e+00    2.732934e-02    1.572250e-03    2.348857e-05
 3-iter: 3.450593e+02    1.448223e-02    1.999584e-04    8.731810e-08
 4-iter: 7.526539e+07    7.281020e-03    2.549689e-05
 5-iter: 7.946898e+23    3.761838e-03    3.230770e-06
 6-iter:          nan    1.916894e-03
 7-iter:                 9.837258e-04
 8-iter:                 5.030396e-04
 9-iter:                 2.575884e-04
10-iter:                 1.319376e-04
11-iter:                 6.749689e-05
12-iter:                 3.457797e-05
13-iter:                 1.763826e-05

將上述數據取 log (避免尺度問題)繪製成圖如下：

從圖形上看 `g4` 取 log 後的相對誤差下降的速度遠小於 `g3` 和 `g2` ，從數據上看 `g2` 的相對誤差每一次迭代約小 2 倍； `g3` 的相對誤差每一次迭代約小 7 倍； `g4` 的相對誤差每一次迭代約小 252 倍。

從上述來看，這個迭代法在收斂至固定點的過程中，其誤差是有規律的乘上一個係數，同時在不同的函數下，由誤差構成的線斜率也不盡相同。那要怎麼分析且計算收斂方面的數據呢？請參考下數學推論過程。

- [ ] 收斂速率(Order of Convergence)定義：
令數列 $\lbrace x_i \rbrace^n_{i=1}$ 收斂到 $L$，假設有一實數 $\lambda>0$ 和 實數 $\alpha \geq 1$ 使得
$$\displaystyle
\lim_{n\rightarrow\infty}\frac{\lvert x_{n+1} - L \rvert}{\lvert x_n -L \rvert^{\alpha}} =\lambda.
$$ 其中 $\alpha$ 稱作為 **order of convegence**， $\lambda$ 稱作為 **asymptotic error constant**
> 某些書中討論收斂性的時候有些會使用 rate ，有些會使用 order ，而在中文裡頭這兩者有時候都會翻譯成「速率」。為避免誤會，本處採用《[Numerical Analysis 10/e](https://www.amazon.com/Numerical-Analysis-Richard-L-Burden/dp/1305253663)》的定義與英文名詞為基準，使用 $\alpha$ 稱作為 **order of convegence**， $\lambda$ 稱作為 **asymptotic error constant**。

- [ ] 假設二次微分連續函數 $g$ 在區間 $I$ 固定點 $p$ 存在且唯一，同時$g'(p)\neq 0$ ，則固定點為線性收斂，即 $\alpha=1$。
證明：
對固定點 $p$ 做泰勒展開
$$
g(x) = g(p) + g'(p)(x-p)+\dfrac{g''(\xi)}{2}(x-p)^2
$$ 其中 $\xi$ 介於 $x$ 和 $r$ 之間。藉由固定點的定義 $x_{n+1}=g(x_n)$ 和 $g(p)=p$ 得到
$$
x_{n+1} = p + g'(p) (x_n-p)+\dfrac{g''(\xi)}{2}(x_n-p)^2
$$ 兩邊減去 $p$ 後，同除 $(x_n-p)$ 得到
$$
\dfrac{x_{n+1}-p}{x_n-p} = g'(p) + \dfrac{g''(\xi)}{2}(x_n-p)
$$ 當 $n$ 趨近無窮時， $x_n$ 足夠靠近固定點 $p$ 就得到結果
$$
\displaystyle
\lim_{n\rightarrow\infty}\frac{\lvert x_{n+1}-p \rvert}{\lvert x_n-p \rvert} =\lvert g'(p)\rvert.
$$

藉由這個結果，我們可以很粗略的來計算在線性收斂的時候，誤差以 $\lvert g'(p)\rvert$ 照比例遞減。假設當 $x_n$ 足夠靠近固定點時，我們就可以藉由此方程求得遞減比例為：
$$
\dfrac{\lvert x_{n+1}-p \rvert}{\lvert x_n-p \rvert} \approx \lvert g'(p)\rvert.
$$

以上述的例子來看：
$$
\lvert g_2'(1.3625) \rvert \approx 0.511, \quad \lvert g_3'(1.3625) \rvert \approx 0.127.
$$

比較特別的是 $g'_4(p) =0$ 故無法使用泰勒展開的分析做結論。

### 牛頓法
從泰勒定理可以得知，假設函數 $f$ 為二次可微連續函數，和 $p$ 為 $f(p)=0$ 的一解，令 $p_0=p-h$ 為 $p$ 附近之一點，其中 $h$ 為足夠小的正數，則對 $p_0$ 的一階泰勒多項式為：
$$
0 = f(p) = f(p_0+h) =f(p_0) + hf'(p_0)+\dfrac{h^2}{2}f''(\xi).
$$

其中剩餘項 $\frac{h^2}{2}f''(\xi)$ 由均值定理得來， $\xi$ 位於 $p_0$ 和 $p$ 之間。將 $h$ 的一次項移至等號左側
$$
-hf'(p_0) = f(p_0) +\dfrac{h^2}{2}f''(\xi)
$$

只要 $f'(p_0)\neq 0$ 等式兩側同除 $-f'(p_0)$ 得
$$
h=-\dfrac{f(p_0)}{f'(p_0)}-\dfrac{h^2f''(\xi)}{2f'(p_0)}
$$

故 $p_0+h$ 可以寫為
$$
p = p_0+h = p_0 -\dfrac{f(p_0)}{f'(p_0)}-\dfrac{h^2f''(\xi)}{2f'(p_0)} = p_0 -\dfrac{f(p_0)}{f'(p_0)} + \mathcal{O}(h^2)
$$
從上述推導的過程中，得知 $p$ 點可以由附近點 $p_0$ 與該點斜率和函數值來逼近
$$
p \approx p_0 -\dfrac{f(p_0)}{f'(p_0)}
$$

#### 局部收斂

考慮 $f(x)=x^2-N = 0$，得解 $x= \pm \sqrt{N}$　
1. 給一定值 $a$
2. 在 $( a,f(a) )$ 做切線，與 $x$ 軸交於 $b$
3. 在 $( b,f(b) )$ 做切線，與 $x$ 軸交於 $c$
4. 重複上面步驟，可得近似 $\sqrt{N}$ 的值

公式:
$$
b = \frac{a^2+N}{2a} = (a+\frac{N}{a}) \div 2
$$

算術-幾何平均值不等式
$$
\frac{a+b}{2} >= \sqrt{ab}, \space a, b>0 \\
\Rightarrow \frac{a+1}{2} >= \sqrt{a}, \space a>0 \\
$$

練習 LeetCode [69. Sqrt(x)](https://leetcode.com/problems/sqrtx)，對於求開根號的近似值，可用上述[牛頓法](https://en.wikipedia.org/wiki/Newton%27s_method)來求解。

int mySqrt(int n)
{
    if (n == 0)
        return 0;
    if (n < 4)
        return 1;

    unsigned int ans = n / 2;
    if (ans > 65535)  // 65535 = 2^16 - 1
        ans = 65535;
    
    while (ans * ans > n  || (ans + 1) * (ans + 1) <= n)
        ans = (ans + n / ans) / 2;
    return ans;
}

注意到以下:
* `ans` 可初始為 $\dfrac{n}{2}$，是因假設 $\sqrt n = x$，當 $x >= 2$ 時，$n = x^2 >= 2x$ 恆成立，也就是 $\dfrac{n}{2} >= x$，符合的初始狀態
* `if (ans > 65535) ans = 65535;` 是避免 `ans * ans` 可能會溢位從而導致錯誤，因此需要對初始狀態做額外的調整

為了討論的便利，我們將上述實作稱為 `newton`。

## 二進位系統的平方根

參考〈[How to calculate a square root without a calculator](https://www.homeschoolmath.net/teaching/square-root-algorithm.php)〉，其方法的證明可見〈[Why the square root algorithm works](https://www.homeschoolmath.net/teaching/sqr-algorithm-why-works.php)〉，應用到二進位〈[Methods of computing square roots](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_(base_2))〉也相似，我們則可利用該方法把平方根換成較快的逐位元位移的運算。

以下是 $x$ 的平方根 $\sqrt{x}$ 逐位元運算的數學推導。

令 $q_i=\sqrt{x},\ q_i$ 代表到小數點後面 $i$ 位精確值 $i \ge 0$
假設 $1 \le x < 4$，於是 $\sqrt{x}$ 的整數位必定是 1，所以 $q_0=1$。

考慮 $q_{i+1}$，亦即逐位元逼近:
1. 若 $(q_{i} + 2^{-(i+1)})^2 \le x$，則 $q_{i+1}=q_i+2^{-(i+1)}$
- 即 $q_i$ 下一位補 1
2. 若 $(q_{i} + 2^{-(i+1)})^2 > x$，則 $q_{i+1}=q_i$
- 即 $q_i$ 下一位補 0

整理上面第一種情況 $(q_{i} + 2^{-(i+1)})^2 \le x$：
$$
\begin{aligned}
{(q_i + 2^{-(i+1)})}^2 \leq x &\Rightarrow q_i^2 + 2\times q_i \times (2^{-(i+1)}) + (2^{-(i+1)})^2 \leq x \\
&\Rightarrow q_i \times (2^{-i}) + 2^{-2(i+1)} \leq x - q_i^2 \\
&\Rightarrow q_i \times 2 \space + 2^{-(i+1)} \leq (x - q_i^2) \times 2^{(i+1)} \\
&\Rightarrow s_i + r_i \leq x_i, \\
&\begin{split} where \space &s_i = 2 \times q_i, \\
&r_i = 2^{-(i+1)}, \\
&x_i = (x - q_i^2) \times 2^{(i+1)} = (x - q_i^2) \times r_i^{-1}
\end{split}
\end{aligned}
$$

因此
- 若 $s_i + r_i \leq x_i$，則
* $q_{i+1}= q_i + r_i$
* $s_{i+1}=(s_i + r_i) + r_i$
* $x_{i+1}=2x_i - 2(s_i + r_i) = 2[x_i - (s_i + r_i)]$
$$
\begin{split}
s_{i+1} &= 2 \times q_{i+1} \\
&= 2 \times [q_i + 2^{-(i+1)}] \\
&= (2 \times q_i) + 2^{-i} \\
&= s_i + 2^{-i} \\
&= s_i + 2^{-(i+1)} + 2^{-(i+1)} \\
&= (s_i + r_i) + r_i \\
x_{i+1} &=[x - q_{i+1}^2] \times 2^{((i+1)+1)} \\
&= [x - (q_i + r_i)^2] \times 2^{(i+2)} \\
&= [x - (q_i^2 + 2q_ir_i + r_i^2)] \times 2r_i^{-1} \\
&= [x - q_i^2] \times 2r_i^{-1} - 2q_ir_i \times 2r_i^{-1} - r_i^2 \times 2r_i^{-1} \\
&= 2x_i - 2s_i - 2r_i \\
&= 2x_i - 2(s_i + r_i) \\
\end{split}
$$
- 若 $s_i + r_i > x_i$，則
* $q_{i+1} = q_i$
* $s_{i+1}=s_i$
* $x_{i+1}=2x_i$
$$
\begin{split}
s_{i+1} &= 2 \times q_{i+1} \\
&= 2 \times q_i \\
&= s_i \\
x_{i+1} &=[x - q_{i+1}^2] \times 2^{((i+1)+1)} \\
&= [x - q_i^2] \times 2^{(i+2)} \\
&= [x - q_i^2] \times 2r_i^{-1} \\
&= 2x_i
\end{split}
$$

LeetCode [69. Sqrt(x)](https://leetcode.com/problems/sqrtx) 可運用上述想法，搭配二分搜尋法來實作:

typedef union {
    float value;
    uint32_t v_int;
} ieee_float_shape_type;

/* Get a 32 bit int from a float. */
#define EXTRACT_WORDS(ix0, d)        \
    do {                             \
        ieee_float_shape_type ew_u;  \
        ew_u.value = (d);            \
        (ix0) = ew_u.v_int;          \
    } while (0)

int mySqrt(int n)
{ 
    int32_t sign = 0x80000000;
    int32_t ix0, m, i;

    float x = n;
    EXTRACT_WORDS(ix0, x);
    
    /* take care of zero */
    if (ix0 <= 0) {
        if ((ix0 & (~sign)) == 0)
            return x; /* sqrt(+-0) = +-0 */
        if (ix0 < 0)
            return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
    }

    /* normalize x */
    m = (ix0 >> 23);
    if (m == 0) { /* subnormal x */
        for (i = 0; (ix0 & 0x00800000) == 0; i++)
            ix0 <<= 1;
        m -= i - 1;
    }
    m -= 127; /* unbias exponent */
    ix0 = (ix0 & 0x007fffff) | 0x00800000;
    if (m & 1) /* odd m, double x to make it even */
        ix0 += ix0;
    m >>= 1; /* m = [m/2] */
     
    /* binary search 'm' */
    unsigned int head = 1 << m;       // 2^m
    unsigned int tail = 1 << (m + 1);  // 2^(m+1)
    unsigned int mid = (head + tail) >> 1; // same as (2^m + 2^(m+1)) / 2
    while (1) {
        if (n > (mid + 1) * (mid + 1)) {
            head = mid; 
            mid = (head + tail) >> 1;
        } else if (n < mid * mid) {
            tail = mid; 
            mid = (head + tail) >> 1;
        } else
            break;
    }

    return mid;
}

先求出整數 `n` 開根號後結果的 $1.FRACTION \times 2^{EXPONENT}$ 的 $EXPONENT$ 部份，則我們知道 n 開根號後的結果 $k$ 應滿足 $2^{EXPONENT} \leq k < 2^{EXPONENT+1}$，因此後續可以用二分搜尋法找出結果。此實作不需要支援浮點數運算，甚至是整數的除法就可得到結果。

為了討論的便利，我們將上述實作稱為 `ieee754 + binary search`。

巨集 `EXTRACT_WORDS` 能將 32 位元整數轉成 `float`，可用 bitmask 對應各部分的位置，參考底下:
- sign: 1 bit `0x80000000`
- exponent: 8 bits `0x7f800000`
- fraction: 23 bits `0x007fffff`

float 表示式: $S\times1.Frac\times2^{127+E}$
> 延伸閱讀: 〈[你所不知道的 C 語言: 浮點數運算](https://hackmd.io/@sysprog/c-floating-point)〉

將 `float x` 表示為 $x=y\times2^{m}$，當要算 $\sqrt{x}$ 時，就計算$\sqrt{y}\times2^{m/2}$。其中 $2^{m/2}$ 用位移即可達成；而要處理的部分就剩下 $\sqrt{y}$ 的估算。

我們來比較前述二個整數平方根的效率。

首先考慮 n 較小的情況下，即 0 到 1000 的範圍，可見 `ieee754 + binary search` 的執行時間低於 `newton`。

考慮到在數值範圍較大時，二分搜尋可能會導致時間變長，觀察在 `x = 1,000,000` 到 `x = 10,000,000` 間的執行時間差異。在 n 達 1,000,000 時，`ieee754 + binary search` 的執行時間低於 `newton` 實作。

## Fast Inverse Square Root （平方根倒數）

3D 繪圖經常要計算平方根倒數，原理是牛頓法，避開開根號和除法運算，節省很多計算時間。

《雷神之鎚 III》(Quake III) 的[原始程式碼](https://github.com/id-Software/Quake-III-Arena)直到 QuakeCon 2005 才正式釋出，但早在 2002 年，反平方根快速演算法的程式碼就出現在論壇中。最初人們猜測是 John D. Carmack 寫下這段程式碼，但他否定這個觀點，並猜測可能是先前曾幫 id Software 最佳化雷神之鎚的資深程式設計師 Terje Mathisen 寫下這段程式碼，而後者表示，他在 1990 年代初，曾有過類似的實作，確切來說這段程式碼亦非他所作。現在所知的最早實作是由 Gary Tarolli 在 SGI Indigo 中實作的，但他亦坦承他僅對常數 R 的取值做一定的改進，實際上他不是原作者。向發明 MATLAB 而聞名的 Cleve Moler 查證後，Rys Sommefeldt 則認為原始的演算法是 Ardent Computer 公司的 Greg Walsh 所發明，但他也不能確認。


* Wikipedia: [Fast inverse square root](http://en.wikipedia.org/wiki/Fast_inverse_square_root)
* [Chris Lomont 的分析](http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf)
* [A Brief History of InvSqrt](https://mrober.io/papers/rsqrt.pdf)
* [揭祕·變態的平方根倒數演算法](https://i.hsfzxjy.site/uncover-the-secret-of-fast-inverse-square-root-algorithm/)
* [解說錄影](https://youtu.be/g1r3iLejTw0)
* [Fast Inverse Square Root — A Quake III Algorithm](https://youtu.be/p8u_k2LIZyo?si=Km0V1RVe0UT6VKzD)

float InvSqrt(float x)
{  
    float xhalf = 0.5f * x;
    int i = *(int *) &x;
    i = 0x5f3759df - (i >> 1);
    x = *(float *) &i;
    x = x * (1.5f - xhalf * x * x);
    return x;
}