back to table · edit · history · where entries came from · files · download
3613 bytes, as of the version from 2026-11-02 11:02 (current). Recorded here, not run.
"""Independent quadratic-root integrals for T651 (Sage or Python + mpmath).
Run: sage -python integral_checks.py, with computed.jsonl beside this file.
This never calls the analytic period evaluator. It tests the nonzero N(0)
control before evaluating all thirty scalar inputs at 125 working digits,
and compares both the signed periods and the fifteen final determinants.
Quadrature convergence estimates are numerical diagnostics, not proofs.
"""
import json,re,time
from pathlib import Path
import mpmath as mp
def argument(key):
# Parse the exact machine identity independently of analytic.py.
if 'sqrt' not in key:return mp.mpf(key)
u,v,d=re.fullmatch(r'(-?\d+(?:/\d+)?)([+-]\d+(?:/\d+)?)\*sqrt\((\d+)\)',key).groups()
def rat(s):
parts=s.split('/')
return mp.mpf(parts[0])/(mp.mpf(parts[1]) if len(parts)==2 else 1)
return rat(u)+rat(v)*mp.sqrt(int(d))
def conjugate(key):
u,v,d=re.fullmatch(r'(-?\d+(?:/\d+)?)([+-]\d+(?:/\d+)?)\*sqrt\((\d+)\)',key).groups()
return u+('+' if v[0]=='-' else '-')+v[1:]+'*sqrt('+d+')'
# Independently transcribed from Table 1, rows 13--27.
COEFFICIENTS={13:(2,1),14:(2,1),15:(2,1),16:(4,-1),17:(5,1),18:(5,1),
19:(25,-1),20:(2,1),21:(2,1),22:(17,1),23:(49,-1),
24:(11,1),25:(11,1),26:(41,1),27:(89,1)}
def integral(key,working=125,ordinary=False):
mp.mp._tanh_sinh.clear()
with mp.workdps(working):
t=argument(key);k=mp.sign(t)*mp.root(abs(t),3)
def f(theta):
x=mp.exp(mp.j*theta);B=x*x-k*x;root=mp.sqrt(B*B-4*x)
return mp.log(max(abs((-B+root)/2),abs((-B-root)/2)))
tolerance=mp.power(10,-working+15)
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>=24:raise ArithmeticError('Integration did not converge: '+key)
c=(a+b)/2
return integrate(a,c,depth+1)+integrate(c,b,depth+1)
def interval(a,b):
c=(a+b)/2
return (integrate(a,c)+integrate(c,b))/mp.pi
if -1<t<27:
c=mp.acos((k-1)/2);I=interval(0,c);J=interval(c,mp.pi)
result=I+J if ordinary else I-2*J
else:result=interval(mp.mpf(0),mp.pi)
assert mp.isfinite(result)
return mp.nstr(result,working-5)
def agreement(x,y):
with mp.workdps(210):
delta=abs((mp.mpf(x)-mp.mpf(y))/mp.mpf(y))
return float(-mp.log10(delta)) if delta else 999.
def main():
control=integral('0',90,ordinary=True)
with mp.workdps(100):
known=3*mp.sqrt(3)/(4*mp.pi)*(mp.zeta(2,mp.mpf(1)/3)-mp.zeta(2,mp.mpf(2)/3))/9
assert abs(mp.mpf(control)-known)<mp.mpf('1e-80')
print(json.dumps({'control':'N(0)=3sqrt(3)L(2,chi_-3)/(4pi)','computed':control,'known':mp.nstr(known,90),'agreement_digits':agreement(control,known)}),flush=True)
for line in Path(__file__).with_name('computed.jsonl').read_text().splitlines():
row=json.loads(line);key=row['t'];start=time.time()
v=integral(key);w=integral(conjugate(key));a,b=COEFFICIENTS[row['row']]
with mp.workdps(150):value=mp.nstr(a*mp.mpf(v)**2+b*mp.mpf(w)**2,115)
dv=agreement(v,row['periods'][1][0]);dw=agreement(w,row['periods'][1][1]);dd=agreement(value,row['high'])
assert min(dv,dw,dd)>110,(key,dv,dw,dd)
print(json.dumps({'t':key,'row':row['row'],'periods':[v,w],'determinant':value,'agreement_digits':[dv,dw,dd],'working_digits':125,'seconds':time.time()-start}),flush=True)
if __name__=='__main__':main()