高斯-勒讓德演算法

高斯-勒讓德演算法是一種用於計算圓周率（π）的演算法。它以迅速收斂著稱，只需25次迭代即可產生π的4500萬位正確數字。不過，它的缺點是主記憶體密集，因此有時它不如梅欽類公式使用廣泛。

該方法基於德國數學家卡爾·弗里德里希·高斯（Johann Carl Friedrich Gauß，1777–1855）和法國數學家阿德里安-馬里·勒讓德（Adrien-Marie Legendre，1752–1833）的個人成果與乘法和平方根運算之現代演算法的結合。該演算法反覆替換兩個數值的算術平均數和幾何平均數，以接近它們的算術-幾何平均數。

下文的版本也被稱為高斯-歐拉，布倫特-薩拉明（或薩拉明-布倫特）演算法^[1]；它於1975年被理查德·布倫特和尤金·薩拉明獨立發現。日本筑波大學於2009年8月17日宣佈利用此演算法計算出π小數點後2,576,980,370,000位數字，計算結果用波溫演算法檢驗。^[2]

知名的電腦效能測試程式Super PI也使用此演算法。

演算法

1. 設置初始值：
   ${\displaystyle a_{0}=1\qquad b_{0}={\frac {1}{\sqrt {2}}}\qquad t_{0}={\frac {1}{4}}\qquad p_{0}=1.\!}$
2. 反覆執行以下步驟直到${\displaystyle a_{n}\!}$與${\displaystyle b_{n}\!}$之間的誤差到達所需精度：
   ${\displaystyle {\begin{aligned}a_{n+1}&={\frac {a_{n}+b_{n}}{2}},\\b_{n+1}&={\sqrt {a_{n}b_{n}}},\\t_{n+1}&=t_{n}-p_{n}(a_{n}-a_{n+1})^{2},\\p_{n+1}&=2p_{n}.\end{aligned}}}$
3. 則π的近似值為：
   ${\displaystyle \pi \approx {\frac {(a_{n+1}+b_{n+1})^{2}}{4t_{n+1}}}.\!}$

下面給出前三個迭代結果（近似值精確到第一個錯誤的位數）：

${\displaystyle 3.140\dots \!}$

${\displaystyle 3.14159264\dots \!}$

${\displaystyle 3.1415926535897932382\dots \!}$

該演算法具有二階收斂性，本質上說就是演算法每執行一步正確位數就會加倍。

數學背景

算術-幾何平均數的極限

a₀和b₀兩個數的算術-幾何平均數，是通過計算它們的序列極限得到的：

${\displaystyle {\begin{aligned}a_{n+1}&={\frac {a_{n}+b_{n}}{2}},\\b_{n+1}&={\sqrt {a_{n}b_{n}}},\end{aligned}}}$

兩者匯聚於同一極限。
若${\displaystyle a_{0}=1\!}$且${\displaystyle b_{0}=\cos \varphi \!}$，則極限為${\displaystyle {\pi \over 2K(\sin \varphi )}\!}$，其中${\displaystyle K(k)\!}$為第一類完全橢圓積分：

${\displaystyle K(k)=\int _{0}^{\pi \over 2}{\frac {d\theta }{\sqrt {1-k^{2}\sin ^{2}\theta }}}.\!}$

若${\displaystyle c_{0}=\sin \varphi \!}$，${\displaystyle c_{i+1}=a_{i}-a_{i+1}\!}$，則

${\displaystyle \sum _{i=0}^{\infty }2^{i-1}c_{i}^{2}=1-{E(\sin \varphi ) \over K(\sin \varphi )}\!}$

其中${\displaystyle E(k)\!}$為第二類完全橢圓積分：

${\displaystyle E(k)=\int _{0}^{\pi \over 2}{\sqrt {1-k^{2}\sin ^{2}\theta }}\,d\theta .\!}$

高斯知道以上這兩個結果。^{[3] [4] [5]}

勒讓德恆等式

對於滿足${\displaystyle \varphi +\theta ={1 \over 2}\pi \!}$的${\displaystyle \varphi \!}$與${\displaystyle \theta \!}$，勒讓德證明了以下恆等式：

${\displaystyle K(\sin \varphi )E(\sin \theta )+K(\sin \theta )E(\sin \varphi )-K(\sin \varphi )K(\sin \theta )={1 \over 2}\pi \!.}$^[3]

高斯-歐拉法

${\displaystyle \varphi =\theta ={\pi \over 4}\!}$的值可以代入勒讓德恆等式，且K、E的近似值可通過${\displaystyle a_{0}=1\!}$與${\displaystyle b_{0}=\sin {\pi \over 4}={\frac {1}{\sqrt {2}}}\!}$的算術-幾何平均數的序列項得到。^[6]

參考文獻

1. Brent, Richard, Old and New Algorithms for pi, Letters to the Editor, Notices of the AMS 60(1), p. 7
2. [円周率計算桁数世界記録樹立について (PDF). [2009-11-29]. （原始內容存檔 (PDF)於2020-09-25）. 円周率計算桁数世界記録樹立について]
3. Brent, Richard, Traub, J F , 編, Multiple-precision zero-finding methods and the complexity of elementary function evaluation, Analytic Computational Complexity (New York: Academic Press), 1975: 151–176 [8 September 2007], （原始內容存檔於2008-07-23）
4. Salamin, Eugene, Computation of pi, Charles Stark Draper Laboratory ISS memo 74–19, 30 January, 1974, Cambridge, Massachusetts
5. Salamin, Eugene, Computation of pi Using Arithmetic-Geometric Mean, Mathematics of Computation 30 (135), 1976, 30 (135): 565–570, ISSN 0025-5718 引文格式1維護：日期與年 (link)
6. Adlaj, Semjon, An eloquent formula for the perimeter of an ellipse, Notices of the AMS 59(8), p. 1096