pythonで円周率を求める

pythonで円周率を求めるアルゴリズム
　
○ライプニッツの方法
収束が遅い。
100万回計算して、6桁くらい。

pi4 = 0
for i in range(1000000):
    pi4 += (1 / (i * 4 + 1) - 1 / (i * 4 + 3))
print(pi4 * 4)    # 3.141592153589902

　
○区分求積を使う。
収束が速い。
25回の計算で15桁くらいまで出ます。

pi4 = 0
for i in range(1, 50, 2):
    pi4 += (-1)**(i//2) * (1/2**i + 1/3**i) / i
print(pi4 * 4)        # 3.1415926535897922

　
○ガウス・ルジャンドルの方法
収束が非常に速く、4回も計算すれば15桁まで出ます。

a, b, t, p = 1, 1/2**0.5, 1/4, 1
for i in range(4):
    a, b, t, p = (a+b)/2, (a*b)**0.5, t-p*((b-a)/2)**2, p*2 
    print((a + b)**2 / (4 * t))

# 3.1405792505221686
# 3.141592646213543
# 3.141592653589794
# 3.141592653589794

　
○モンテカルロ法
乱数を使って円周率を求める方法。
原点中心、半径１の円とその円に外接する正方形の面積比がπ/4であることを利用し、1万個の座標を乱数で生成して円周率を求めてみた。
サンプルコード(1)　randomモジュールで。

import random

inner = 0
for i in range(10000):
    x, y = random.random(), random.random()
    if x**2 + y**2 < 1:
        inner += 1
print(inner *4 / 10000)   # 3.1446

　
サンプルコード(2)　NumPyのrandomで。

np.random.seed(1)               # 乱数の初期化。好みの数で。
arr = np.random.rand(2, 10000)  # 2行10000列の乱数ndarray
distance = (arr**2).sum(axis=0) # 全要素を2乗し縦に合計すると原点からの距離
count = (distance < 1).sum()    # 距離が1未満か判断しTrueの個数を調べる
print(4 * count / 10000)        # 3.1348

　
1万個の散布図を描くとキャンバスが埋まってしまうので、1000個で描画

　

import numpy as np
import matplotlib.pyplot as plt

c1 = plt.Circle((0, 0), radius=1, fc="None", ec="r", linewidth=3)
ax = plt.gca()
ax.add_patch(c1)

plt.axis("scaled")
plt.xlim(0, 1)
plt.ylim(0, 1)

x = np.random.rand(1000)
y = np.random.rand(1000)
plt.scatter(x, y, marker="x")
plt.show()

　
○逆三角関数を使う方法
cos(π)＝－1なので
π＝arccos(-1)

import math
print(math.acos(-1)) # 3.141592653589793

　
参考：
円周率を計算してみよう - 金沢工業大学
Wikipedia　ガウス・ルジャンドルのアルゴリズム