# buffon.py 
# Fig. 22.12 

import math, random
import matplotlib.pyplot as plt

NUM_NEEDLES = 1000000

def truncL(x):  
  # truncate to the left for [-1, 2)
  return math.trunc(x+2) - 2

ns = range(1, NUM_NEEDLES+1)
probs = []
count = 0
for n in ns: 
  x = random.random() 
  y = random.random() 
  angle = math.pi*random.random() 
  b = math.cos(angle)/2
  c = math.sin(angle)/2
  if truncL(x-c) != truncL(x+c): 
    # vertical intersects
    count += 1
  if truncL(y-b) != truncL(y+b): 
    # horizonatal intersects
    count += 1
  prob = count/n
  probs.append(prob)
print("(no. intersects)/(total throws) approaches:", prob)
print(f"PI estimate: {4/prob:.5f}")

plt.semilogx(ns[100:], probs[100:], 'b')
plt.xlabel('No. of throws')
plt.ylabel('Prob. of intersection')
plt.xlim(100, NUM_NEEDLES)
plt.title(f'Final prob. = {prob:5}')
plt.grid('on')
plt.show()