reset()
def f(m,n,k) :
ell = m/n
if (ell / (k+1) >= 1) : return 1
return (2*n*exp(-ell)*(ell^k)*(ell+1)) / (factorial(k)*(k+1-ell))
def bin_search_f_k_bdd(m,n,t,k) :
while(true) :
if (f(m,n,k) <= 2^(-t)) : return k
else : k = (k<<1)
def bin_search_f_k_range(m,n,t,k_min,k_max) :
if (k_max == 1) : return k_max
while(true):
if (k_max == k_min+1) : return k_max
k_mid = floor((k_max + k_min)>>1)
if (f(m,n,k_mid) <= 2^(-t)) : k_max = k_mid
else: k_min = k_mid
def bin_search_f_n_bdd(m,n,t,k) :
while(true) :
if (f(m,n,k) <= 2^(-t)) : return n
else : n = (n<<1)
def bin_search_f_n_range(m,k,t,n_min,n_max) :
if (n_max == 2) : return n_max
while(true):
if (n_max == n_min+1) : return n_max
n_mid = floor((n_max + n_min)>>1)
if (f(m,n_mid,k) <= 2^(-t)) : n_max = n_mid
else: n_min = n_mid
for (m,k,t) in [ (10^3,3,8), (10^3,5,8), (10^3,9,8), (10^3,17,8)
, (10^4,3,8), (10^4,5,8), (10^4,9,8), (10^4,17,8)
, (10^6,3,8), (10^6,5,8), (10^6,9,8), (10^6,17,8)
] :
n_max = bin_search_f_n_bdd(m,1,t,k)
n_min = max(2,n_max>>1)
n_just = bin_search_f_n_range(m,k,t,n_min,n_max)
print "m=", m
print "\t t=", 2^(-t)
print "\t k=", k
print "\t n=", n_just