Pythonでとりあえず何かしたいと思い立ち、Geminiに相談したところ出してきたアイデアの一つが円周率の計算だった。AIは統計的にそれっぽい案を提示しただけではあると思うが、アイデアを出す部分は人間の強みとして維持していたいところ。機械に操られている感もあるが、円周率の計算をやってみた。
 |
| Gemini生成 |
円周率の計算方法は思いのほか色々な手法が提案されており、様々なサイトで紹介されている。このBlogはあくまで自分宛のメモ程度なので、もっと役立ちそうな内容は他のサイトを参照：
色々試すと面白そうであるが、特に面白そうだったモンテカルロシミュレーションと、高速に計算できるというチュドノフスキー（Chudnovsky）の公式を試してみた。それぞれの特徴は以下の通り：
- モンテカルロシミュレーション： 「確率的・幾何学的」アプローチ。概念が直感的にわかりやすいが、精度を1桁上げるのに試行回数を100倍にする必要があり収束が極めて遅い。
- Chudnovskyの公式： 「解析的・数理的」アプローチ。数式自体は複雑だが、1項ごとに14桁伸びるため収束が異次元に速い。
モンテカルロシミュレーション
モンテカルロシミュレーションは、”ランダム・サンプリングを繰り返し実行することによって、ある範囲の結果が発生する可能性を算出する計算アルゴリズム”（IBM）、”不確実な事象の起こり得る結果を予測する数学的手法”（AWS）、などと定義されている。乱数を使って繰り返し実行した結果を統計的に分析して、求めたい事象の近似値を求める手法と理解した。
円周率を求めるには、正方形の中にランダムに点を打ち、内接する円に入った点の数と天の総数の比から円周率の近似値を求めることができる。
\[{\pi }≒ 4*\frac{円の中の点}{点の総数} \]
import numpy as np
import matplotlib.pyplot as plt
points = 1000
x = np.random.uniform(-1.0, 1.0, points)
y = np.random.uniform(-1.0, 1.0, points)
inside = (x**2 + y**2 <= 1.0)
pi_est = 4 * np.sum(inside) / points
# --- グラフの描画 ---
plt.figure(figsize=(6, 6))
# 円の内側の点を青（blue）、外側の点を赤（red）でプロット
plt.scatter(x[inside], y[inside], color='blue', s=1, label='Inside Circle')
plt.scatter(x[~inside], y[~inside], color='red', s=1, label='Outside Circle')
# 形状を正円にするためにアスペクト比を1:1に固定
plt.gca().set_aspect('equal', adjustable='box')
# グラフのタイトルとラベル
plt.title(f'Monte Carlo Pi Simulation (Total: {points:,} points)\nEstimated Pi = {pi_est:.5f}', fontsize=12)
plt.xlabel('X')
plt.ylabel('Y')
plt.legend(loc='upper right')
plt.grid(True, linestyle='--', alpha=0.5)
# 描画結果を表示
plt.show() |
| 1,000点の描画。結果は3.204 |
 |
| 10,000点の描画。結果は3.1324 |
import numpy as np
import matplotlib.pyplot as plt
num_points = 1000
num_trials = 100000
pi_estimates = []
for _ in range(num_trials):
x = np.random.uniform(-1.0, 1.0, num_points)
y = np.random.uniform(-1.0, 1.0, num_points)
inside = (x**2 + y**2 <= 1.0)
pi_est = 4 * np.sum(inside) / num_points
pi_estimates.append(pi_est)
mean_pi = np.mean(pi_estimates)
std_pi = np.std(pi_estimates)
# --- グラフの描画 ---
plt.figure(figsize=(10, 6))
# 近似値の散らばりをヒストグラムで描画
plt.hist(pi_estimates, bins=40, color='skyblue', edgecolor='black', alpha=0.7)
# 真の円周率(π)のライン
plt.axvline(np.pi, color='red', linestyle='--', linewidth=2, label=f'True Pi ({np.pi:.5f})')
# 複数回試行した平均値のライン
plt.axvline(mean_pi, color='green', linestyle=':', linewidth=2, label=f'Mean Estimate ({mean_pi:.5f})')
plt.title(f'Distribution of Pi Estimates\n({num_points} points per trial, Repeated {num_trials} times)', fontsize=12)
plt.xlabel('Estimated Pi Value')
plt.ylabel('Frequency')
plt.legend()
plt.grid(True, alpha=0.3)
# 描画結果を表示
plt.show()
print(f"真の円周率: {np.pi:.6f}")
print(f"{num_trials}回の平均値: {mean_pi:.6f}")
print(f"標準偏差（ばらつきの大きさ）: {std_pi:.6f}") |
| 100点を1,000回実行。平均：3.151520 |
 |
| 1,000点を10,000回実行。平均： 3.141897 |
 |
| 10,000点を1,000,000回実行。平均： 3.141598 |
チュドノフスキーの公式
公式を調べると様々な形態で紹介されていたが、とりあえず以下のようにまとめることができる：
\[\frac{1}{\pi }=12\sum _{k=0}^{\infty }\frac{(-1)^{k}(6k)!(545140134k+13591409)}{(3k)!(k!)^{3}(640320)^{3k+3/2}}\]
この公式の最大の特徴は収束の圧倒的な速さ。計算を1項進めるごとに、円周率の正しい桁数が約14.18桁増える。数回の反復計算だけで数百〜数千桁の精度が得られる、とのこと。
import math
import decimal
from decimal import Decimal, getcontext
def pi_chudnovsky(digits):
getcontext().prec = digits + 10
pi = Decimal(0)
k = 0
while True:
numerator = (Decimal(-1) ** k
* decimal.Decimal(math.factorial(6 * k))
* (13591409 + 545140134 * k))
denominator = (decimal.Decimal(math.factorial(3 * k))
* decimal.Decimal(math.factorial(k)) ** 3
* Decimal(640320) ** (3 * k + Decimal(3) / 2))
term = numerator / denominator
pi += term
if abs(term) < Decimal(1) / (10 ** digits):
break
k += 1
pi = 1 / (12 * pi)
return +pi
print(pi_chudnovsky(100))1,000桁の計算も5秒弱で完了：
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989380952574
今回は円周率をPythonで計算してみたが、モンテカルロシミュレーションを活かせるようなネタがあったら試してみたい。
