Power-of-2-choices
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)

For n= 1000 t= 20 [500, 324, 160, 56, 16, 7, 4, 3, 2] End correction: 3 Chain length: 13 For n= 10000 t= 20 [5000, 2730, 887, 130, 12, 3, 2] End correction: 2 Chain length: 10 For n= 1000000 t= 20 [500000, 252282, 64935, 4563, 49, 2, 1] End correction: 2 Chain length: 10 For n= 1000000000000 t= 20 [500000000000, 250002280047, 62502414628, 3906880306, 15284290, 319, 1, 0] End correction: 0 Chain length: 8 For n= 10000 t= 40 [5000, 2827, 1008, 184, 24, 7, 4] End correction: 4 Chain length: 12 For n= 1000000 t= 40 [500000, 253229, 65957, 4849, 68, 4, 2] End correction: 3 Chain length: 11 For n= 1000000000 t= 40 [500000000, 250101971, 62608027, 3934489, 16416, 10, 1] End correction: 2 Chain length: 10 For n= 1000000000000 t= 40 [500000000000, 250003224474, 62503414812, 3907141400, 15294858, 357, 1] End correction: 2 Chain length: 10