π を明示的には使わない，ビュフォンの針のシミュレーション(Python)

ビュフォンの針のシミュレーションプログラムはたくさん書かれているが，プログラム内で定数 π を使っているものがほとんどである。

π を使わないシミュレーションプログラムを書いたので記録に残しておく。

方針は，原点を左下の頂点とする一辺  1 の正方形内のランダムな点の座標（x2，y2） のうち，原点からの距離が 1 以下の点について sin(theta) は三角関数の定義により y2/sqrt(x2**2 + y2**2) であることを使う。

以下のプログラムで DEBUG = True として実行すると，theta の分布は一様分布になっていることが確認できる。

なお，プログラム中の REMARK と書かれた if 文を外すと，半径が 1 より大きい点も採用されてしまうので，theta は一様分布にはならない。

import random
import math

DEBUG = False
if DEBUG:
    import numpy as np
    import matplotlib.pyplot as plt

n = 100000
d = 10
l = d / 2
i = 0
cross = 0
if DEBUG: theta = np.zeros(n)
while i < n:
  y = random.uniform(0, d / 2) # 針の中心から最も近い平行線までの距離
  # 0 <= theta <= π/2 の一様分布を求める
  x2 = random.uniform(0, 1) # 一辺 1 の正方形内の座標 (x2, y2)
  y2 = random.uniform(0, 1)
  r2 = x2**2 + y2**2 # 原点から (x2, y2) までの平方距離
  if r2 < 1: ## REMARK　円の内部の点に限定する
    r = math.sqrt(r2) # 原点から (x2, y2) までの距離
    sine_theta = y2 / r # sine の定義により
    if DEBUG: theta[i] = math.asin(sine_theta)
    i += 1
    if y <= l / 2 * sine_theta: # 平行線と交わる
      cross += 1

if DEBUG: plt.hist(theta, bins=np.arange(0, np.pi/2, np.pi/60))
print(2 * l * n / (cross * d))