これの続き。ガウス・ルジャンドル法を使うとより早く円周率計算ができる。
初期値を
\begin{align}
a_0=1 \hspace{10mm} b_0=\frac{1}{\sqrt{2}} \hspace{10mm} t_0=\frac{1}{4} \hspace{10mm} p_0 = 1
\end{align}
を与えて次を反復計算すればいい。
\begin{align}
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} & = 2 p_n
\end{align}
その後得られたパラメータを使えば
\begin{align}
\pi \approx \frac{(a+b)^2}{4t}
\end{align}
で近似できる。
import numpy as np
n = 3
a = np.zeros(n + 1)
b = np.zeros(n + 1)
t = np.zeros(n + 1)
p = np.zeros(n + 1)
Pi = np.zeros(n + 1)
a[0] = 1
b[0] = 1 / np.sqrt(2)
t[0] = 1 / 4
p[0] = 1
Pi[0] = (a[0] + b[0]) / (4 * t[0])
for i in range(n):
a[i + 1] = (a[i] + b[i]) / 2
b[i + 1] = np.sqrt(a[i] * b[i])
for i in range(n):
t[i + 1] = t[i] - p[i] * (a[i] - a[i + 1]) ** 2
p[i + 1] = 2 * p[i]
Pi[i + 1] = (a[i + 1] + b[i + 1]) ** 2 / (4 * t[i + 1])
print(Pi)
コメント