ガウス=ルジャンドルのアルゴリズム

ガウス=ルジャンドルのアルゴリズム英語: Gauss–Legendre algorithm)は、円周率を計算する際に用いられる数学の反復計算アルゴリズムである。円周率を計算するものの中では非常に収束が速く、2009年にこの式を用いて2,576,980,370,000桁(約26000億桁)の計算がなされた。

このアルゴリズムはカール・フリードリヒ・ガウスアドリアン=マリ・ルジャンドルがそれぞれ別個に研究したものである。これは2つの数値の算術幾何平均を求めるために、それぞれの数値を算術平均(相加平均)と幾何平均(相乗平均)で置き換えていくものである。

アルゴリズム

これによる円周率の計算方法は以下の通りである。

初期値の設定

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}

反復式

ab が希望する精度(桁数)になるまで以下の計算を繰り返す。小数第n位まで求めるとき log2n 回程度の反復でよい。

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

π の算出

円周率 π は、abt を用いて以下のように近似される。

π ( a + b ) 2 4 t {\displaystyle \pi \approx {\frac {(a+b)^{2}}{4t}}}

最初の3回の反復で得られる数値(最後の桁は真値とは異なる)は以下の通りである。

3.140 {\displaystyle 3.140\dots }  (小数点以下2桁目までが正しい)
3.14159264 {\displaystyle 3.14159264\dots }  (小数点以下7桁目までが正しい)
3.1415926535897932382 {\displaystyle 3.1415926535897932382\dots } (小数点以下18桁目までが正しい)

この計算過程は二次収束する。つまり反復のたびに正しい桁数が直前のもののほぼ2倍になるのである。ガウス自身もこの式を用いて反復を4回まで行って12桁まで正しいことを確認したことが知られている。

Python3による実装例

#!/usr/bin/python3
import sys
from decimal import *

def picalc():
  a = Decimal(1)
  b = Decimal(1) / Decimal(2).sqrt()
  t = Decimal(1) / Decimal(4)
  p = Decimal(1)
  r = Decimal(0)
  rn = Decimal(3)
  while r != rn:
    r = rn
    an = (a + b) / 2
    bn = (a * b).sqrt()
    tn = t - p * (a - an) * (a - an)
    pn = 2 * p
    rn = ((a + b) * (a + b)) / (4 * t)
    a = an
    b = bn
    t = tn
    p = pn
  return rn

if __name__ == "__main__":
  if len(sys.argv) < 2:
    getcontext().prec = 10000
  else:
    try:
      getcontext().prec = int(sys.argv[1])
    except:
      print("引数が不正です。桁数を数値で指定して下さい。")
      sys.exit(1)
  print(picalc())

数学的背景

算術幾何平均

二つの数 a 0 , b 0 {\displaystyle a_{0},b_{0}} 算術幾何平均とは、数列

a n + 1 = a n + b n 2 ,   b n + 1 = a n b n {\displaystyle a_{n+1}={\frac {a_{n}+b_{n}}{2}},\ b_{n+1}={\sqrt {a_{n}b_{n}}}}

の極限のことである。両数列は同一の極限値を持つ。

a 0 = 1 ,   b 0 = cos φ {\displaystyle a_{0}=1,\ b_{0}=\cos \varphi } のとき、算術幾何平均は π 2 K ( sin φ ) {\displaystyle {\frac {\pi }{2K(\sin \varphi )}}} となる。ここで、 K ( k ) {\displaystyle K(k)} 第一種完全楕円積分である。また、 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-{\frac {E(\sin \varphi )}{K(\sin \varphi )}}}

が成立する。ここで、 E ( k ) {\displaystyle E(k)} 第二種完全楕円積分である。

ガウスはこれらの結果を知っていたとされる[1][2][3]

ルジャンドル恒等式

φ + θ = π 2 {\displaystyle \varphi +\theta ={\frac {\pi }{2}}} なる φ , θ {\displaystyle \varphi ,\theta } について、ルジャンドルは次の恒等式を得た。

K ( sin φ ) E ( sin θ ) + K ( sin θ ) E ( sin φ ) K ( sin φ ) K ( sin θ ) = π 2 . {\displaystyle K(\sin \varphi )E(\sin \theta )+K(\sin \theta )E(\sin \varphi )-K(\sin \varphi )K(\sin \theta )={\frac {\pi }{2}}.}

この式は、任意の φ {\displaystyle \varphi } について、以下のようにも書き表すことができる。

K ( sin φ ) { E ( cos φ ) K ( cos φ ) } + K ( cos φ ) E ( sin φ ) = π 2 . {\displaystyle K(\sin \varphi )\left\{E(\cos \varphi )-K(\cos \varphi )\right\}+K(\cos \varphi )E(\sin \varphi )={\frac {\pi }{2}}.}

初等的な証明

ガウス=ルジャンドルのアルゴリズムは楕円モジュラー関数を用いずに示すこともできる。初等的な積分を用いた証明が Lord (1992)[4]、Milla (2019)[5] によりなされている。

脚注

[脚注の使い方]
  1. ^ Brent, Richard P. (1976-01-01). Traub, J. F.. ed (英語). Analytic Computational Complexity. Academic Press. pp. 151–176. doi:10.1016/b978-0-12-697560-4.50014-9. ISBN 978-0-12-697560-4. http://www.sciencedirect.com/science/article/pii/B9780126975604500149 
  2. ^ Salamin, Eugene (1976-07). “Computation of π Using Arithmetic-Geometric Mean”. Mathematics of Computation 30 (135): 565. doi:10.2307/2005327. https://www.jstor.org/stable/2005327?origin=crossref. 
  3. ^ Salamin, Eugene, Computation of pi, Charles Stark Draper Laboratory ISS memo 74–19, 30 January 1974, Cambridge, Massachusetts
  4. ^ Lord, Nick (1992-07). “Recent Calculations of p: The Gauss-Salamin Algorithm”. The Mathematical Gazette 76 (476): 231. doi:10.2307/3619132. https://www.jstor.org/stable/3619132?origin=crossref. 
  5. ^ Milla, Lorenz (2019), Easy Proof of Three Recursive π-Algorithms, arXiv:1907.04110

関連項目

  • スーパーπ
  • 表示
  • 編集