fork download
  1. # Python 3.9
  2.  
  3. import random
  4. import math
  5. import time
  6.  
  7. twopi = 2.0*math.pi
  8.  
  9. def trial():
  10. """
  11. place a point at (1,0) and four other points randomly on the unit circle
  12. return boolean a**2<b*c; a = area of petal touching point A, b = etc.
  13. """
  14. theta = [0.] + sorted([twopi*random.random() for _ in range(4)])
  15. A,B,C,D,E = [(math.cos(theta[i]), math.sin(theta[i])) for i in range(5)]
  16. return area(A,B,C,D,E)**2 < area(B,C,D,E,A) * area(C,D,E,A,B) # a**2<bc
  17. #return area(A,B,C,D,E)**2 < area(C,D,E,A,B) * area(D,E,A,B,C) # a**2<cd
  18.  
  19. def area(A,B,C,D,E):
  20. """
  21. return the area of triangle APQ (the petal touching A) where A,B,C,D,E
  22. are points on the perimeter of the circle, in *counterclockwise* order
  23. """
  24. P = intersection(A,C,B,E)
  25. Q = intersection(A,D,B,E)
  26. return 0.5 * ( A[0]*(P[1]-Q[1]) + P[0]*(Q[1]-A[1]) + Q[0]*(A[1]-P[1]) )
  27.  
  28. def intersection(P,Q,R,S):
  29. """
  30. return the point of intersection of line(P,Q) and line(R,S)
  31. each line is specified by a pair of points, e.g. P=(x,y) in rect. coords.
  32. """
  33. t = ( (S[0]-R[0])*(R[1]-P[1]) - (S[1]-R[1])*(R[0]-P[0]) ) / \
  34. ( (S[0]-R[0])*(Q[1]-P[1]) - (S[1]-R[1])*(Q[0]-P[0]) )
  35. return (P[0] + t*(Q[0]-P[0]), P[1] + t*(Q[1]-P[1]))
  36.  
  37. nsamp = 10**5 #<--------- sample-size
  38. print(f"{nsamp=}")
  39.  
  40. start = time.time()
  41.  
  42. nsucc = sum(trial() for _ in range(nsamp))
  43. print(f"{nsucc=}")
  44.  
  45. phat = nsucc / nsamp
  46. print(f"{phat=:.6f}")
  47.  
  48. sdev = math.sqrt(phat*(1-phat)/nsamp)
  49. print(f"99.9% CI for p: {phat:.6f} +- {3.29053*sdev:.6f} \
  50. = ({phat-3.29053*sdev:.6f}, {phat+3.29053*sdev:.6f})")
  51.  
  52. tsecs = time.time() - start
  53. print(f"time elapsed \
  54. = ({tsecs:9.2f} s) = ({tsecs/60:7.2f} m) = ({tsecs/3600:5.2f} h)")
Success #stdin #stdout 1.25s 10252KB
stdin
Standard input is empty
stdout
nsamp=100000
nsucc=49986
phat=0.499860
99.9% CI for p: 0.499860 +- 0.005203 =  (0.494657, 0.505063)
time elapsed = (     1.23 s) = (   0.02 m) = ( 0.00 h)