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

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

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

>初期値を
>```math
>\begin{align*}
>a_0 = 1,\ \ b_0 = \frac{1}{\sqrt{2}}
>\end{align*}
>```
>として，下の漸化式で数列を生成する。
>```math
>\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*}
>```
>このとき，数列
>```math
>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.](https://arxiv.org/abs/1907.04110)

## 使用する命題たち
本当はこれらも証明する必要があるのですが，ここでは詳しくは書きません。

### 命題１
$a_n,b_n$はともに同じ値（$:=m$）に収束する。

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

### 命題２
```math
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}$ が成り立つ。

:::note
$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}$の証明では，被積分関数列が一様収束するので積分と極限の順序が交換できることを用いる。
:::

### 命題３
```math
L(a,b) := \int_{0}^{\pi/2} \frac{\cos^2\phi}{\sqrt{a^2\cos^2\phi + b^2\sin^2\phi}} d\phi
```
とすると，
```math
\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.
```

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

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

:::note
命題２，３を使う。
:::

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

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

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

### 命題６
ベータ関数$\left(B(u,v) := \int_{0}^{1}t^{u-1}(1-t)^{v-1}dt\right)$について，次が成り立つ。
```math
B(u,v) = \frac{\Gamma(u)\Gamma(v)}{\Gamma(u+v)}\ \ \ (u,v>0)
```

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

### 命題７
```math
L(\sqrt{2},1)\cdot I(\sqrt{2},1) = \frac{\pi}{4}
```
:::note
命題５，６を使う。
:::

## 証明
$\lambda$を定数として
```math
\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*}
```
だから，命題７より

```math
\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*}
```
これを用いると，命題４より
```math
\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}$ の極限なので，
```math
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
```
したがって
```math
\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)
```