import numpy as np 
import scipy.stats as st 
import matplotlib.pyplot as plt 
import scipy.optimize as so

np.random.seed(42)

p = np.pi/4
sigma = ( p*(1-p) )**(1/2)
gamma = p*(1-p)*((1-p)**2 + p**2)**(1/3)
K_BE = 0.4748

alpha = 1e-4
eps = 1e-2

# Central Limit theorem bound
N_CLT = int(np.ceil( ( sigma*st.norm.ppf(1-alpha/2)/eps )**2 ))
print("N_CLT \t= ", N_CLT)

# Berry Essen bound
# Note - this interval needs to contain the root of f, otherwise the bisection method
# iterations will fail as coded here. 
Nmin = 100000
Nmax = 300000
maxIter = 1000
Ns = [Nmin, Nmax]
f = lambda N : 1-alpha/2-st.norm.cdf(np.sqrt(N)*eps/sigma) \
        + K_BE*gamma**3/sigma**3/np.sqrt(N)
i=1
while True and i < maxIter:
    Nmid = int(np.ceil((Nmin+Nmax)/2))
    if Nmid==Nmax-1 or Nmid==Nmin+1 or Nmid==Nmax or Nmid==Nmin:
        break
    if f(Nmin)*f(Nmid) < 0:
        Nmax = Nmid-1
    elif f(Nmax)*f(Nmid) < 0:
        Nmin = Nmid+1
        
    else:
        import sys
        sys.exit('Error during bisection method')
    i = i+1

if f(Nmid)<0:
    Nmid=Nmid+1

N_BE = Nmid
print("N_BE \t= ", N_BE)

# Chebyshev bound
N_CHEB = int(np.ceil(  sigma**2/(alpha*eps**2) ) )
print("N_CHEB \t= ", N_CHEB)