"""Independent quadratic-root quadrature for T650 (mpmath, Python 3).

Run: python3 integral_checks.py > integral-checks.jsonl
This checks the defining period directly, sharing no hypergeometric code with
our generator. Numerical quadrature estimates are not rigorous error bounds.
"""
import json,re,time
import mpmath as mp
KEYS=('24','18-6*sqrt(3)','-27/2+27/2*sqrt(5)',
 '-4320+1944*sqrt(5)','-864+216*sqrt(17)',
 '-163296+35640*sqrt(21)','4+4*sqrt(33)',
 '-16740+2916*sqrt(33)','-55296+8640*sqrt(41)',
 '-13500000+1431000*sqrt(89)')

def rational(s):
 a,sep,b=s.partition('/')
 return mp.mpf(a)/mp.mpf(b) if sep else mp.mpf(a)

def argument(key):
 if 'sqrt' not in key:return rational(key)
 a,b,d=re.fullmatch(r'(-?\d+(?:/\d+)?)([+-]\d+(?:/\d+)?)\*sqrt\((\d+)\)',key).groups()
 return rational(a)+rational(b)*mp.sqrt(int(d))

def periods(key,working):
 with mp.workdps(working):
  t=argument(key);k=mp.sign(t)*mp.root(abs(t),3)
  c=mp.acos((k-1)/2)
  calls=0
  def f(theta):
   nonlocal calls
   calls+=1
   x=mp.exp(mp.j*theta);B=x*x-k*x;D=mp.sqrt(B*B-4*x)
   return mp.log(max(abs((-B+D)/2),abs((-B-D)/2)))
  tolerance=mp.power(10,-working+20)
  def integrate(a,b,depth=0):
   value,error=mp.quad(f,[a,b],error=True,maxdegree=8)
   if mp.isfinite(value) and error<=tolerance*(b-a)/mp.pi:return value
   if depth>=22:raise ArithmeticError('Quadrature did not converge')
   m=(a+b)/2
   return integrate(a,m,depth+1)+integrate(m,b,depth+1)
  def interval(a,b):
   m=(a+b)/2
   return integrate(a,m)+integrate(m,b)
  I=interval(0,c)/mp.pi;J=interval(c,mp.pi)/mp.pi
  return I,J,calls

def run():
 # Known nonzero control before believing the numerical instrument.
 I,J,calls=periods('0',70)
 with mp.workdps(80):
  known=3*mp.sqrt(3)/(4*mp.pi)*(mp.zeta(2,mp.mpf(1)/3)-mp.zeta(2,mp.mpf(2)/3))/9
  err=abs(I+J-known)
  assert err<mp.mpf('1e-60')
  print(json.dumps({'control':'N(0)=3sqrt(3)L(2,chi_-3)/(4pi)','value':mp.nstr(I+J,65),'known':mp.nstr(known,65),'absolute_error':mp.nstr(err,6)}),flush=True)
 for key in KEYS:
  start=time.time();low=periods(key,135);high=periods(key,185)
  with mp.workdps(200):
   a=low[0]-2*low[1];b=high[0]-2*high[1]
   agreement=-mp.log10(abs(a-b)/abs(b))
   assert agreement>115
   print(json.dumps({'t':key,'modified':mp.nstr(b,175),'I':mp.nstr(high[0],175),'J':mp.nstr(high[1],175),'ordinary':mp.nstr(high[0]+high[1],175),'agreement_digits':mp.nstr(agreement,8),'quadrature_evaluations':[low[2],high[2]],'seconds':time.time()-start}),flush=True)
if __name__=='__main__':run()