code

import math

from decimal import *

context = Context(prec=1000)

setcontext(context)

def sum1(D, Delta): #this is for the first term in equation (16): sum 1/y * log(y+1) when y goes through D+1 to Floor[Delta]

temp = Decimal(0)

for i in range(int(D+1), int(Delta+1)):

temp = temp + Decimal(1.0/i*math.log(i+1, 2))

return temp

def bigsum(c, D, Delta): # this evaluates the right hand side of equation (16), assuming we know D

c = Decimal(c)

D = Decimal(D)

Delta = Decimal(Delta)

val = Decimal(Decimal(1.0/2)/(c-1)*sum1(D, Delta) + (Decimal(1.0/2)/(c-1)-1/D)*Decimal(math.log(D+1,2)) + 1 -Decimal(1.0-math.log(1+2**(-1/D),2))*(c-1)/D)

return val

def optBigsum(c): # This evaluates the maximum possible value for the right hand side of equation (16), among all feasible values of D (when fixed C)

Dlow = math.ceil(2*(c-1)) # This is the lower bound for D. Note that D is an integer, thus D >= ceiling(2*(C-1))

Delta = math.floor(2*c*(c-1)) # This is the upper bound for D. Note that D is an integer, thus D <= floor(2*C*(C-1))

if Delta < Dlow:

return (None, None)

else:

value = 0

tempD = Dlow

for D in range(int(Dlow), int(Delta+1)):

bigs = bigsum(c,D,Delta)

if bigs> value:

value = bigs

tempD = D

return (round(value,3), int(tempD))

for counter in range(1, 3): # This function prints out for each C, the upper bounf for the left hnad side of equation (16), log(f(Ct, (C-1)t, t), i.e., the exponent we want to bound.

c = Decimal(4.85+ 0.01*counter)

print('value C is {}; the value is: {}; the optimized D is: {}'.format(str(round(c,3)),str(optBigsum(c)[0]),str(optBigsum(c)[1])) )