高斯-勒讓德演算法 是一種用於計算圓周率 (π)的演算法 。它以迅速收斂 著稱,只需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 也使用此演算法。
演算法
設定初始值:
a
0
=
1
b
0
=
1
2
t
0
=
1
4
p
0
=
1.
{\displaystyle a_{0}=1\qquad b_{0}={\frac {1}{\sqrt {2}}}\qquad t_{0}={\frac {1}{4}}\qquad p_{0}=1.\!}
反覆執行以下步驟直到
a
n
{\displaystyle a_{n}\!}
與
b
n
{\displaystyle b_{n}\!}
之間的誤差到達所需精度:
a
n
+
1
=
a
n
+
b
n
2
,
b
n
+
1
=
a
n
b
n
,
t
n
+
1
=
t
n
−
p
n
(
a
n
−
a
n
+
1
)
2
,
p
n
+
1
=
2
p
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}}}
則π的近似值為:
π
≈
(
a
n
+
1
+
b
n
+
1
)
2
4
t
n
+
1
.
{\displaystyle \pi \approx {\frac {(a_{n+1}+b_{n+1})^{2}}{4t_{n+1}}}.\!}
下面給出前三個迭代結果(近似值精確到第一個錯誤的位數):
3.140
…
{\displaystyle 3.140\dots \!}
3.14159264
…
{\displaystyle 3.14159264\dots \!}
3.1415926535897932382
…
{\displaystyle 3.1415926535897932382\dots \!}
該演算法具有二階收斂性,本質上說就是演算法每執行一步正確位數就會加倍。
數學背景
算術-幾何平均數的極限
a0 和b0 兩個數的算術-幾何平均數 ,是通過計算它們的序列極限得到的:
a
n
+
1
=
a
n
+
b
n
2
,
b
n
+
1
=
a
n
b
n
,
{\displaystyle {\begin{aligned}a_{n+1}&={\frac {a_{n}+b_{n}}{2}},\\b_{n+1}&={\sqrt {a_{n}b_{n}}},\end{aligned}}}
兩者匯聚於同一極限。
若
a
0
=
1
{\displaystyle a_{0}=1\!}
且
b
0
=
cos
φ
{\displaystyle b_{0}=\cos \varphi \!}
,則極限為
π
2
K
(
sin
φ
)
{\displaystyle {\pi \over 2K(\sin \varphi )}\!}
,其中
K
(
k
)
{\displaystyle K(k)\!}
為第一類完全橢圓積分 :
K
(
k
)
=
∫
0
π
2
d
θ
1
−
k
2
sin
2
θ
.
{\displaystyle K(k)=\int _{0}^{\pi \over 2}{\frac {d\theta }{\sqrt {1-k^{2}\sin ^{2}\theta }}}.\!}
若
c
0
=
sin
φ
{\displaystyle c_{0}=\sin \varphi \!}
,
c
i
+
1
=
a
i
−
a
i
+
1
{\displaystyle c_{i+1}=a_{i}-a_{i+1}\!}
,則
∑
i
=
0
∞
2
i
−
1
c
i
2
=
1
−
E
(
sin
φ
)
K
(
sin
φ
)
{\displaystyle \sum _{i=0}^{\infty }2^{i-1}c_{i}^{2}=1-{E(\sin \varphi ) \over K(\sin \varphi )}\!}
其中
E
(
k
)
{\displaystyle E(k)\!}
為第二類完全橢圓積分 :
E
(
k
)
=
∫
0
π
2
1
−
k
2
sin
2
θ
d
θ
.
{\displaystyle E(k)=\int _{0}^{\pi \over 2}{\sqrt {1-k^{2}\sin ^{2}\theta }}\,d\theta .\!}
高斯知道以上這兩個結果。[ 3]
[ 4]
[ 5]
勒讓德恆等式
對於滿足
φ
+
θ
=
1
2
π
{\displaystyle \varphi +\theta ={1 \over 2}\pi \!}
的
φ
{\displaystyle \varphi \!}
與
θ
{\displaystyle \theta \!}
,勒讓德證明了以下恆等式:
K
(
sin
φ
)
E
(
sin
θ
)
+
K
(
sin
θ
)
E
(
sin
φ
)
−
K
(
sin
φ
)
K
(
sin
θ
)
=
1
2
π
.
{\displaystyle K(\sin \varphi )E(\sin \theta )+K(\sin \theta )E(\sin \varphi )-K(\sin \varphi )K(\sin \theta )={1 \over 2}\pi \!.}
[ 3]
高斯-歐拉法
φ
=
θ
=
π
4
{\displaystyle \varphi =\theta ={\pi \over 4}\!}
的值可以代入勒讓德恆等式,且K、E的近似值可通過
a
0
=
1
{\displaystyle a_{0}=1\!}
與
b
0
=
sin
π
4
=
1
2
{\displaystyle b_{0}=\sin {\pi \over 4}={\frac {1}{\sqrt {2}}}\!}
的算術-幾何平均數的序列項得到。[ 6]
參考文獻
^ Brent, Richard , Old and New Algorithms for pi , Letters to the Editor, Notices of the AMS 60(1), p. 7
^ [円周率計算桁数世界記録樹立について (PDF) . [2009-11-29 ] . (原始內容存檔 (PDF) 於2020-09-25). 円周率計算桁数世界記録樹立について ]
^ 3.0 3.1 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)
^ Salamin, Eugene , Computation of pi , Charles Stark Draper Laboratory ISS memo 74–19, 30 January, 1974, Cambridge, Massachusetts
^ Salamin, Eugene , Computation of pi Using Arithmetic-Geometric Mean, Mathematics of Computation 30 (135), 1976, 30 (135): 565–570, ISSN 0025-5718
^ Adlaj, Semjon, An eloquent formula for the perimeter of an ellipse , Notices of the AMS 59(8), p. 1096