back to table · edit · history · where entries came from · files · download
7879 bytes, as of the version from 2026-11-10 02:09 (current). Recorded here, not run.
"""Independent checks for numberdb.org/T703.
Run: sage -python checks.py
Download generate.py and check.py alongside this file. This checks all 494
entries before submission. For a database check call run(stored_values), where
stored_values maps (Delta,N,r,s) to an interval from the database.
The independent mpmath computation sums the logarithm index k after performing
the n sum analytically, using the geometric moments 1/(1-t), t/(1-t)^2,
t(1+t)/(1-t)^3. It does not call the generator's summand or tail bound.
Its two precisions provide agreement evidence, not rigorous error bounds.
The independent square-lattice beta value and direct Eisenstein-Kronecker
sums check conventions outside the defining logarithmic series.
"""
import math
import json
import sys
import platform
import mpmath as mp
import numberdb.sage as numberdb
from sage.rings.real_arb import RealBallField
from sage.rings.rational_field import QQ
from generate import Table,compute,DISCRIMINANTS
from check import exactness,measure,names_its_rings,agrees_with
def midpoint_mp(v):
q=v.mid().exact_rational()
return mp.mpf(str(q.numerator()))/mp.mpf(str(q.denominator()))
def independent(D,N,r,s,dps):
with mp.workdps(dps):
a=mp.pi*mp.sqrt(-D);xi=mp.mpf(r)/N
q=(-1 if D%4 else 1)*mp.exp(-a)
x=mp.exp(-a*xi+mp.j*mp.pi*(2*s-int(D%4!=0)*r)/N)
rate=min(xi,1-xi) if r else mp.mpf(1)
K=int(mp.ceil((dps+35)*mp.log(10)/(a*rate)))+4
total=-a**3*xi**2*(1-xi)**2/4
xp=mp.mpc(1);qp=mp.mpf(1)
for k in range(1,K+1):
xp*=x;qp*=q
den=1-qp
S1=qp/den**2;S2=qp*(1+qp)/den**3
# The xi=0 branch avoids subtracting 1 from S0 at tiny q.
if r:
term=xp*(xi**2/den+2*xi*S1+S2)+(S2-2*xi*S1+xi**2*qp/den)/xp
else:
term=(xp+1/xp)*S2
total-=a*a*mp.re(term)/k
return +total
def canonical(N,r,s):
d=math.gcd(N,math.gcd(r,s));N//=d;r//=d;s//=d
r,s=min((r,s),((-r)%N,(-s)%N))
return N,r,s
def assert_narrow(v,bits=315):
assert v.is_finite() and not v.contains_zero(),str(v)
assert v.accuracy()>bits,(v.accuracy(),v)
def identities(values,R):
count=0;largest=R(0)
for D in DISCRIMINANTS:
c=-(R.pi()*R(-D).sqrt())**3/120
origin=values[D,1,0,0]
for M in range(1,7):
total=R(0)
for r in range(M):
for s in range(M):
total+=values[(D,)+canonical(M,r,s)]
target=origin/(M*M)+(M*M-R(1)/(M*M))*c
residual=total-target
assert residual.contains_zero(),(D,M,residual)
assert residual.rad_as_ball()<R(10)**(-90),(D,M,residual)
largest=max(largest,R(residual.rad_as_ball()));count+=1
return {'full_torsion_distribution_checks':count,'max_residual_radius':str(largest)}
def lattice_checks():
# Fixed double-precision lattice sums are only a coarse convention test.
# The omitted square-shell tail is <=48*y^3/(pi*M^2): the least Gram
# eigenvalue is >=1/2 in this CM slice, and there are 8k shell points.
import numpy as np
M=160;axis=np.arange(-M,M+1,dtype=float)
m,n=np.meshgrid(axis,axis,indexing='ij');nz=(m!=0)|(n!=0);m=m[nz];n=n[nz]
radius2=m*m+n*n
control=float(np.sum(1/radius2**2))
target=float(4*mp.zeta(2)*mp.catalan)
control_bound=4/M**2+1e-10
assert abs(control-target)<control_bound
results=[]
for D,N,r,s in [(-4,1,0,0),(-4,2,1,0),(-4,3,1,1),(-3,2,1,0),(-3,3,1,2),(-7,4,1,1)]:
y=math.sqrt(-D)/2;re=n-int(D%4!=0)*m/2;im=m*y;den=re*re+im*im
phase=np.cos(2*math.pi*(n*r-m*s)/N)
value=float(y**3/math.pi*np.sum(phase*(1/den**2+2*(re*re-im*im)/den**3))-(2*math.pi*y)**3/120)
exact=float(independent(D,N,r,s,70));bound=48*y**3/(math.pi*M*M)+1e-9
error=abs(value-exact);assert error<bound,(D,N,r,s,error,bound)
results.append({'params':[D,N,r,s],'absolute_error':error,'tail_bound_plus_roundoff_allowance':bound})
return {'control_value':control,'control_expected':target,'control_error_bound':control_bound,'checks':results,'limitation':'Numerical roundoff has a conservative 1e-9 allowance, not a ball proof; this is only a coarse normalization check.'}
def run(stored_values=None):
mp.mp.dps=200;R=RealBallField(768)
# The nonzero control tests the independent geometric-moment identity.
t=mp.mpf(1)/3
geom=sum(mp.mpf(k*k)*t**k for k in range(1,600))
assert abs(geom-mp.mpf(3)/2)<mp.mpf('1e-190')
# mpmath evaluates Hurwitz zeta with s first.
beta_mp=(mp.zeta(2,mp.mpf(1)/4)-mp.zeta(2,mp.mpf(3)/4))/16
assert abs(beta_mp-mp.catalan)<mp.mpf('1e-190')
square_target=4*mp.zeta(2)*beta_mp/mp.pi-mp.pi**3/15
square=compute(-4,1,0,0,digits=150)
assert abs(midpoint_mp(square)-square_target)<mp.mpf('1e-170')
enum=list(Table().enumerate());assert len(enum)==494
counts={N:sum(p['Delta']==-3 and p['N']==N for p in enum) for N in range(1,7)}
assert counts=={1:1,2:3,3:4,4:6,5:12,6:12}
values={};details=[];max_mp=mp.mpf(0);max_two=mp.mpf(0)
by_key={}
for index,p in enumerate(enum):
key=tuple(int(p[k]) for k in Table.parameters)
v,diag=compute(*key,diagnostics=True);assert_narrow(v)
hi=compute(*key,digits=130);assert_narrow(hi,415)
assert v.overlaps(hi)
lo_mp=independent(*key,150);hi_mp=independent(*key,190)
mp_error=abs(lo_mp-hi_mp)/abs(hi_mp)
ball_error=abs(midpoint_mp(hi)-hi_mp)/abs(hi_mp)
assert mp_error<mp.mpf('1e-120'),(key,mp_error)
assert ball_error<mp.mpf('1e-155'),(key,ball_error)
max_mp=max(max_mp,ball_error);max_two=max(max_two,mp_error)
# Enclose the independent 190-digit midpoint conservatively at 160
# relative digits; its convergence is checked, not certified.
q=QQ(str(hi_mp._mpf_[1]))*QQ(2)**int(hi_mp._mpf_[2])*(-1 if hi_mp._mpf_[0] else 1)
by_key[key]=R(q).add_error(abs(R(q))*R(10)**(-160))
D,N,r,s=key
inv=compute(D,N,(-r)%N,(-s)%N,digits=110)
conj=compute(D,N,r,(int(D%4!=0)*r-s)%N,digits=110)
assert_narrow(inv);assert_narrow(conj);assert v.overlaps(inv) and v.overlaps(conj)
values[key]=v;details.append({'params':list(key),**diag,'value':str(v)})
if (index+1)%100==0:print('checked',index+1,flush=True)
assert not exactness(values)
assert not agrees_with(values,lambda k:by_key[k])
id_result=identities(stored_values or values,R)
if stored_values is not None:
assert set(stored_values)==set(values)
for k,v in stored_values.items():
assert_narrow(v);assert v.overlaps(values[k]) and v.overlaps(by_key[k])
assert abs(midpoint_mp(stored_values[-4,1,0,0])-square_target)<mp.mpf('1e-98')
measurements=measure(values)
report={'table':'T703','skill_manifest':'../20261110T015857Z-build-skills/manifest.json','skill_sha256':'886390324dd1971e381a43cf3b46b822ab0bc67d62f31d8481e29ec79556746b','entries':len(values),'counts_per_CM_order':counts,'minimum_relative_bits':min(x['relative_bits'] for x in details),'terms_each_side_range':[min(x['terms_each_side'] for x in details),max(x['terms_each_side'] for x in details)],'max_independent_relative_discrepancy':mp.nstr(max_mp,12),'max_two_precision_relative_discrepancy':mp.nstr(max_two,12),'geometric_moment_control':'sum k^2/3^k=3/2','square_lattice_origin':mp.nstr(square_target,110),'identities':id_result,'symmetry_checks':2*len(values),'size':measurements,'stored_values_checked':stored_values is not None,'versions':{'python':platform.python_version(),'numberdb':numberdb.__version__,'mpmath':mp.__version__},'lattice_checks':lattice_checks(),'entries_evidence':details}
return report
if __name__=='__main__':
report=run()
print('EVIDENCE_JSON '+json.dumps(report))