旧式SEの計算機室 平行線の間隔を\(d\)、針の長さを\(l (l\leq d)\)とし、針の始端Aのｘ座標、ｙ座標は一様乱数で与えます。
針の終端Bを次のようにして求めます。
始点Aを中心とした単位円内にランダムに点を取り、Aからの距離をdx、dyとします。
三平方の定理から、r を求められるので、\(\sin\theta\)、\(\cos\theta\)を求めることができます。（\(\theta\)は、X軸と針の角度）
\(\sin\theta = \large{\frac{dy}{r}}\)、\(\cos\theta = \large{\frac{dx}{r}}\)
Aの位置を\((x,y)\)とすると、終端Bの位置は、\((x+l\large{\frac{dx}{r}} \normalsize{,y+l} \large{\frac{dy}{r}}\))
針が線と交差するのは
・AかBのy座標が、いずれかの平行線のy座標と等しい
・AとBのy座標をd（間隔）で割った商が異なる
のいずれかのときです。
平行線の間隔\(d\)を10、針の長さ\(l\)を5とし、200×200の平面上でＮ回（Ｎ本）試行するプログラムを掲載します。
'''
ビュホンの針で円周率計算(三角関数は使わない）
平行線の間隔 10cm
針の長さ 5cm
200cm×200cm内で（N回 or N本）試行
'''
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
# 試行回数 or 本数の指定
N=1000
# 線と交差する針 ゼロクリア
count_red=0
# 針の始端セット （一様乱数）
x=np.random.uniform(5.0,195.0,N)
y=np.random.uniform(5.0,195.0,N)
# グラフのサイズ指定
plt.figure(figsize=(10,10))
# 全部の針の終端を計算 → 描画
for i in range(N):
r=2 # 単位円内にランダムに点を取る (0,0)は除く
while r>1:
dx=np.random.uniform(-1.0,1.0)
dy=np.random.uniform(-1.0,1.0)
if dx==0 and dy==0:
r=2
else:
r=math.sqrt(dx*dx+dy*dy)
x_end=x[i]+5*dx/r #針の終端の座標
y_end=y[i]+5*dy/r
X=[x[i],x_end] #針の座標
Y=[y[i],y_end]
if y[i]%10==0 or y_end%10==0: #端が線上 (赤）
plt.plot(X,Y,color="r")
count_red+=1
else:
if y[i]//10!=y_end//10: #線と交差 (赤）
plt.plot(X,Y,color="r")
count_red+=1
else:
plt.plot(X,Y,color="c") #線に触れない （シアン）
print("N = ",N)
print("線と交わる針：",count_red)
print("円周率の近似値:",N/count_red)
#平行線描画(黒）
X=[0,210]
for i in range(21):
Y=[i*10,i*10]
plt.plot(X,Y,color="k")
ax=plt.gca()
ax.xaxis.set_major_locator(ticker.MultipleLocator(10)) #X軸 目盛り指定
ax.yaxis.set_major_locator(ticker.MultipleLocator(10)) #Y軸 目盛り指定
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
