LoginSignup
1
1

More than 1 year has passed since last update.

【円周率計算】ガウス=ルジャンドルのアルゴリズムの証明

Posted at

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

ガウス=ルジャンドルのアルゴリズムは、円周率の計算に用いられるアルゴリズムです。
円周率計算のアルゴリズムの中でもトップクラスの収束速度なので、精度がよいです。

そのアルゴリズムがこちら。

初期値を

\begin{align*}
a_0 = 1,\ \ b_0 = \frac{1}{\sqrt{2}}
\end{align*}

として,下の漸化式で数列を生成する。

\begin{align*}
a_n &= \frac{a_{n-1}+b_{n-1}}{2}\\
b_n &= \sqrt{a_{n-1} \cdot b_{n-1}}\\
c_n &= a_n^2 - b_n^2
\end{align*}

このとき,数列

p_n := \frac{(a_{n} + b_{n})^2}{1 - 2\sum_{j=1}^{n}2^jc_j^2}

は円周率$\pi$に収束する。

これが正しいことの証明方法はいくつか与えられているそうですが,2019年にLorenz Millaによってなされた初等的(高校生でもとける,一部に大学範囲あり)な証明を簡単に紹介します(詳細は以下の論文を見てください)。

参考にさせていただいた論文はこちら。
Milla, L. (2019). Easy Proof of Three Recursive π-Algorithms--Einfacher Beweis dreier rekursiver π-Algorithmen. arXiv preprint arXiv:1907.04110.

使用する命題たち

本当はこれらも証明する必要があるのですが,ここでは詳しくは書きません。

命題1

$a_n,b_n$はともに同じ値($:=m$)に収束する。

$a_n,b_n$の単調性から示す。
$c_n \to 0$から収束先が同じことを示す。

命題2

I(a,b) := \int_{0}^{\pi/2} \frac{d\phi}{\sqrt{a^2\cos^2\phi + b^2\sin^2\phi}}

とすると,任意の自然数$n$で $I(a_n,b_n) = I(a_0,b_0) = \displaystyle\frac{\pi}{2m}$ が成り立つ。

$t=b\tan\phi$とおいて整理したあと,$x=\frac{1}{2}\left(t-\frac{ab}{t}\right)$とおいて計算。
途中で,$f(x)=\frac{(t^2+a^2)(t^2+b^2)}{t^2}(=(2x)^2+(a+b)^2)$とおくとよい。
$I(a,b) = \frac{\pi}{2m}$の証明では,被積分関数列が一様収束するので積分と極限の順序が交換できることを用いる。

命題3

L(a,b) := \int_{0}^{\pi/2} \frac{\cos^2\phi}{\sqrt{a^2\cos^2\phi + b^2\sin^2\phi}} d\phi

とすると,

\left\{
\begin{array}{ll}
L(a,b) + L(b,a) = I(a,b) & \cdots (1)\\
L(a,b) - L(b,a) = \displaystyle\frac{a-b}{a+b}L(b_1,a_1) & \cdots (2)
\end{array}\right.

命題2と同じような置き方をして計算する。

命題4

$S:=\sum_{j=1}^{\infty}2^jc_j^2$ とおくと,$2c_0^2L(a,b) = (c_0^2 - S)I(a.b)$

命題2,3を使う。

命題5

ガンマ関数$\left(\Gamma(x) := \int_{0}^{\infty}t^{x-1}e^{-t}dt\right)$について,次が成り立つ。

\Gamma(x+1) = x\Gamma(x)\ \ \ (x>0),\ \ \ \Gamma\left(\frac{1}{2}\right) = \sqrt{\pi},\ \ \ 

部分積分の公式。
2つ目は,二乗したものを計算する。

命題6

ベータ関数$\left(B(u,v) := \int_{0}^{1}t^{u-1}(1-t)^{v-1}dt\right)$について,次が成り立つ。

B(u,v) = \frac{\Gamma(u)\Gamma(v)}{\Gamma(u+v)}\ \ \ (u,v>0)

ヤコビアンをかけるのを忘れずに。

命題7

L(\sqrt{2},1)\cdot I(\sqrt{2},1) = \frac{\pi}{4}

命題5,6を使う。

証明

$\lambda$を定数として

\begin{align*}
I(\lambda a, \lambda b) &= \int_{0}^{\pi/2} \frac{d\phi}{\sqrt{\lambda^2a^2\cos^2\phi + \lambda^2b^2\sin^2\phi}} \\
&= \frac{1}{|\lambda|}\int_{0}^{\pi/2} \frac{d\phi}{\sqrt{a^2\cos^2\phi + b^2\sin^2\phi}}\\
&= \frac{1}{|\lambda|} I(a,b)\\

L(\lambda a, \lambda b) &= \int_{0}^{\pi/2} \frac{\cos^2\phi}{\sqrt{\lambda^2a^2\cos^2\phi + \lambda^2b^2\sin^2\phi}}d\phi \\
&= \frac{1}{|\lambda|}\int_{0}^{\pi/2} \frac{\cos^2\phi}{\sqrt{a^2\cos^2\phi + b^2\sin^2\phi}} d\phi\\
&= \frac{1}{|\lambda|} L(a,b)\\
\end{align*}

だから,命題7より

\begin{align*}
L\left(1,\frac{1}{\sqrt{2}}\right) \cdot I\left(1,\frac{1}{\sqrt{2}}\right)
&= L\left(\frac{1}{\sqrt{2}}\cdot\sqrt{2}, \frac{1}{\sqrt{2}}\cdot 1\right) \cdot I\left(\frac{1}{\sqrt{2}}\cdot\sqrt{2}, \frac{1}{\sqrt{2}}\cdot 1\right)\\
&= (\sqrt{2})^2 L(\sqrt{2},1) \cdot I(\sqrt{2},1)\\
&= 2\cdot\frac{\pi}{4} \\
&= \frac{\pi}{2}
\end{align*}

これを用いると,命題4より

\begin{align*}
2c_0^2L(a,b)I(a,b) &= (c_0^2 - S)I(a.b)^2\\
2\cdot \frac{1}{2} \cdot L\left(1,\frac{1}{\sqrt{2}}\right) I\left(1,\frac{1}{\sqrt{2}}\right) &= \left(\frac{1}{2} - S\right) I\left(1,\frac{1}{\sqrt{2}}\right)^2\\
\frac{\pi}{2} &= \left(\frac{1}{2} - S\right) I\left(1,\frac{1}{\sqrt{2}}\right)^2\\
\frac{\pi}{2} &= \left(\frac{1}{2} - S\right)\cdot\left(\frac{\pi}{2m}\right)^2\\
1 &= \left(\frac{1}{2} - S\right)\cdot\frac{\pi}{2m^2}\\
\pi &= \frac{4m^2}{1-2S} = \frac{4m^2}{1-2\sum_{j=1}^{\infty}2^jc_j^2}
\end{align*}

となる。$m$ は $a_n = \frac{a_{n-1} + b_{n-1}}{2}$ の極限なので,

4m^2 = 4\left(\frac{\lim a_n + \lim b_n}{2}\right)^2 = (\lim a_n + \lim b_n)^2 = \lim(a_n+b_n)^2

したがって

\pi = \lim_{n\to\infty} \frac{(a_{n} + b_{n})^2}{1 - 2\sum_{j=1}^{n}2^jc_j^2}

(おまけ)Pythonコード

gauss_legendre.py
# ガウス=ルジャンドルのアルゴリズム
a0 = 1.0
b0 = 0.5**0.5
c0 = 0.5
p0 = 1.0

def A(n):
    if n==0:
        return a0
    else:
        return (A(n-1) + B(n-1)) / 2

def B(n):
    if n==0:
        return b0
    else:
        return (A(n-1) * B(n-1))**0.5

def C2(n):
    if n==0:
        return c0
    else:
        return A(n)**2 - B(n)**2
    
def S(n):
    res = 0
    for j in range(1,n+1):
        res += 2**j * C2(j)
    return res

N = 3
pi = (A(N) + B(N))**2 / (1 - 2*S(N))
print(pi)
1
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
1
1