# Python 3.9

import random
import math 
import time

twopi = 2.0*math.pi

def trial():
    """
    place a point at (1,0) and four other points randomly on the unit circle
    return boolean a**2<b*c; a = area of petal touching point A, b = etc.
    """
    theta = [0.] + sorted([twopi*random.random() for _ in range(4)])
    A,B,C,D,E = [(math.cos(theta[i]), math.sin(theta[i])) for i in range(5)]
    return area(A,B,C,D,E)**2 < area(B,C,D,E,A) * area(C,D,E,A,B) # a**2<bc
    #return area(A,B,C,D,E)**2 < area(C,D,E,A,B) * area(D,E,A,B,C) # a**2<cd

def area(A,B,C,D,E):
    """
    return the area of triangle APQ (the petal touching A) where A,B,C,D,E  
    are points on the perimeter of the circle, in *counterclockwise* order
    """
    P = intersection(A,C,B,E)
    Q = intersection(A,D,B,E)
    return 0.5 * ( A[0]*(P[1]-Q[1]) + P[0]*(Q[1]-A[1]) + Q[0]*(A[1]-P[1]) )

def intersection(P,Q,R,S):
    """
    return the point of intersection of line(P,Q) and line(R,S)
    each line is specified by a pair of points, e.g. P=(x,y) in rect. coords.
    """
    t = ( (S[0]-R[0])*(R[1]-P[1]) - (S[1]-R[1])*(R[0]-P[0]) ) / \
        ( (S[0]-R[0])*(Q[1]-P[1]) - (S[1]-R[1])*(Q[0]-P[0]) )
    return (P[0] + t*(Q[0]-P[0]), P[1] + t*(Q[1]-P[1]))

nsamp = 10**5    #<--------- sample-size
print(f"{nsamp=}")

start = time.time() 

nsucc = sum(trial() for _ in range(nsamp)) 
print(f"{nsucc=}")

phat = nsucc / nsamp 
print(f"{phat=:.6f}")

sdev = math.sqrt(phat*(1-phat)/nsamp) 
print(f"99.9% CI for p: {phat:.6f} +- {3.29053*sdev:.6f} \
=  ({phat-3.29053*sdev:.6f}, {phat+3.29053*sdev:.6f})")

tsecs = time.time() - start
print(f"time elapsed \
= ({tsecs:9.2f} s) = ({tsecs/60:7.2f} m) = ({tsecs/3600:5.2f} h)")