back to table · edit · history · where entries came from · files · download
2516 bytes, as of the version from 2026-11-02 11:01 (current). Recorded here, not run.
"""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()