reset()

# Assuming that 1 >= q >= p >0
# Using the D_KL version of the Chernoff Bound
def deviation_bound(p,q,n) :
    q_bar = 1 - q
    p_bar = 1 - p
    return ( (q/p)^q * (q_bar/p_bar)^q_bar )^(-n)

# Assuming that 1 >= q >= p > 0
# Using the D_KL version of the Chernoff Bound
def log_deviation_bound(p,q,n) :
    q_bar = 1 - q
    p_bar = 1 - p
    S1 = 0 if (q==0) else (-n)*q*log(q/p)
    S2 = 0 if (q_bar==0) else (-n)*q_bar*log(q_bar/p_bar)
    return  S1 + S2

# Given p, n, and eps, find smallest q such that log_deviation_bound(p,q,n) <= log(eps), up to granularity
def binary_search_q(p,n,gran,eps) :
    if (eps<=0 or eps >=1) : return -1  ### Sanity check of the input
    if (log_deviation_bound(p,1,n) > log(eps) ) : return -1 ### Sanity check at the right boundary case 
    q_min = p
    q_max = 1.0
    ### Assumption: At q_min deviation probabiliy is > eps, and at q_max deviation probability is <= eps
    while(true) :
        if (q_max <= q_min + gran) : return q_max
        q_mid = (q_min + q_max)/2
        #print q_mid
        val = log_deviation_bound(p,q_mid,n)
        if ( val > log(eps) ) : q_min = q_mid
        else : q_max = q_mid



# Given n, eps, find the ``recursion of weights''
def recursive_probability_list(n,eps,init_balls) :
    gran = 1/n
    ball_list = [ floor(n/init_balls) ]
    p_list = [ ball_list[-1]/n ]
    old_ball_num = ball_list[-1]
    #print old_ball_num
    while (true) :
        p = p_list[-1]^2
        q = binary_search_q(p,n,gran,eps)
        new_ball_num = floor(q*n)
        to_abort = (old_ball_num == new_ball_num)
        if (to_abort) : return (ball_list,p_list)
        #p_list.append(q)
        p_list.append(new_ball_num/n)
        ball_list.append(new_ball_num)
        old_ball_num = new_ball_num
        #print old_ball_num
        if (new_ball_num==0) : return (ball_list,p_list)

def edge_case(n,a,eps) : 
    if (a==0) : return 0 
    return ceil(log((a^2)/eps)/log(n))
    #p = (a/n)^2
    #q = (4/n^2)
    #return 1+min(ceil(log((eps*q)/(n*p^2)) / log(((n^2)*(q^2))/2)), floor(n/a))


def chain_calculate(n,t,init_balls) : 
    print "For n=" , n, "t=", t 
    gran = 1/n
    eps = 2^(-t)
    
    (ball_list,p_list) = recursive_probability_list(n,eps,init_balls)
    print "\t" , ball_list
    at_end = edge_case(n,ball_list[-1],eps)
    print "\t" , "End correction:" , at_end
    print "\t", "Chain length:", init_balls + len(ball_list) - 1 + at_end + (-1 if (ball_list[-1]==0) else 0)

#n = 10^(4)
#t = 20
#init_balls=2


#gran = 1/n
#eps = 2^(-t)


#(ball_list,p_list) = recursive_probability_list(n,eps,init_balls)
#print ball_list
#print p_list
#print "Chain Length:", init_balls + (len(ball_list)-1) + ball_list[-1]


for (n,t,init_balls) in [(10^3,20,2), (10^4,20,2), (10^6,20,2), (10^(12),20,2), (10^4,40,2), (10^6,40,2), (10^9,40,2), (10^(12),40,2)] :
    chain_calculate(n,t,init_balls)
    #gran = 1/n
    #eps = 2^(-t)

    #(ball_list,p_list) = recursive_probability_list(n,eps,init_balls)
    #print "For n=" , n, "t=", t,
    #print "Chain Length:", init_balls + (len(ball_list)-1) + (ball_list[-1]-1),
    #print ball_list
    #print edge_case(n,ball_list[-1],eps)