back to table · edit · history · where entries came from · files · download
7376 bytes, as of the version from 2026-11-02 01:23 (current). Recorded here, not run.
"""Independent checks for T663. Run: sage -python checks.py [--stored].
Requires generate.py next to this file, numberdb and mpmath. --stored reads
all live rows with NUMBERDB_API_KEY when the table is still a private draft.
No write or review action is performed. JSON on stdout is the numerical record.
"""
import sys, json, math, time, platform
import numberdb.sage as numberdb
import mpmath as mp
from generate import ARGUMENTS, compute, BertinP
def q_parameter(tau):
q = mp.exp(2j * mp.pi * tau)
w = mp.exp(mp.pi*1j*tau) * (mp.qp(q)*mp.qp(q**6)/(mp.qp(q**2)*mp.qp(q**3)))**6
return w
def modular(k, working=140):
"""Bertin Theorem 1.1, with fixed CM tau or an independently solved tau."""
with mp.workdps(working):
known = {0: (-3+mp.sqrt(-3))/6, 2: (-2+mp.sqrt(-2))/6,
3: (-3+mp.sqrt(-15))/12, 6: 1j/mp.sqrt(6),
10: 1j/mp.sqrt(2), 18: 1j*mp.sqrt(mp.mpf(5)/6)}
if k in known:
tau = known[k]
else:
target = (k-mp.sqrt(mp.mpf(k)**2-4))/2
guesses = {1: mp.mpc('-0.4','0.28'), 4: mp.mpc('-0.2','0.36'),
5: mp.mpc('-0.13','0.387')}
start = guesses.get(k, 1j*mp.log(k)/mp.pi)
tau = mp.findroot(lambda t: q_parameter(t)-target,
(start, start+mp.mpc('0.001','0.002')),
tol=mp.power(10,-working+10), maxsteps=80)
w=q_parameter(tau)
residual=abs(w+1/w-k)
assert residual < mp.power(10,-working+10)
q=mp.exp(2j*mp.pi*tau); r=abs(q)
assert 0<r<mp.mpf('0.25')
nmax=int(mp.ceil((working+20)*mp.log(10)/(-mp.log(r))))+30
sigma=[0]*(nmax+1)
for d in range(1,nmax+1):
for n in range(d,nmax+1,d):sigma[n]+=d**3
value=-mp.pi*1j*tau
qn=q
for n in range(1,nmax+1):
value+=mp.mpf(sigma[n])/n*(4*qn-8*qn**2+12*qn**3-24*qn**6)
qn*=q
tail=96*r**(nmax+1)*((nmax+1)**2/(1-r)+2*(nmax+1)*r/(1-r)**2+r*(1+r)/(1-r)**3)
assert tail < mp.power(10,-working-10)
return mp.nstr(value.real,working),mp.nstr(residual,6),mp.nstr(tail,6),mp.nstr(tau,25)
def relative_agreement(a,b):
a,b=mp.mpf(a),mp.mpf(b)
if a==b:return mp.inf
return -mp.log10(abs(a-b)/max(1,abs(a),abs(b)))
def controls():
with mp.workdps(140):
# Independent elliptic-RF form of the square-lattice density.
rho=lambda s: mp.elliprf(0,s*s/16,1)/(2*mp.pi**2)
mass=2*mp.quad(rho,[0,4])
moment=2*mp.quad(lambda s:s*s*rho(s),[0,4])
assert abs(mass-1)<mp.mpf('1e-125')
assert abs(moment-4)<mp.mpf('1e-125')
# Smyth's independent Dirichlet L-value, via Hurwitz zeta.
L=(mp.zeta(2,mp.mpf(1)/3)-mp.zeta(2,mp.mpf(2)/3))/9
smyth=3*mp.sqrt(3)*L/(4*mp.pi)
zero=compute(0,150)
score=relative_agreement(zero,smyth)
assert score>125
# Negative control: a perturbation of a size that our comparison rejects.
assert relative_agreement(zero,smyth+mp.mpf('1e-90'))<95
return {'density_mass_error':mp.nstr(abs(mass-1),6),
'density_second_moment_error':mp.nstr(abs(moment-4),6),
'smyth_agreement':mp.nstr(score,8),'perturbation_rejected':True,
'smyth':mp.nstr(smyth,135)}
def exterior_series(k, working=140):
"""Independent closed-walk expansion with exact integer coefficients."""
with mp.workdps(working):
r=mp.mpf(36)/(k*k)
nmax=int(mp.ceil((working+10)*mp.log(10)/(-mp.log(r))))
total=mp.log(k)
power=mp.mpf(1)
for n in range(1,nmax+1):
power/=k*k
c=math.comb(2*n,n)*sum(math.comb(n,j)**2*math.comb(2*j,j) for j in range(n+1))
total-=c*power/(2*n)
# Since every walk has at most 6 choices at each step, c_n <= 36^n.
tail=r**(nmax+1)/(2*(nmax+1)*(1-r))
assert tail<mp.mpf('1e-145')
return mp.nstr(total,working),mp.nstr(tail,6)
def derivative_control():
with mp.workdps(140):
# d/dk Jensen integral at k=6. t=4-s cancels its endpoint square root.
def f(u):
if u==0:return 1/(2*mp.pi*mp.sqrt(4))
s=4-u*u
if s==0:return mp.mpf(0)
rho=mp.elliprf(0,s*s/16,1)/(2*mp.pi**2)
return 2*rho/mp.sqrt(4+u*u)
derivative=mp.quad(f,[0,2,mp.sqrt(8)])
# Watson's closed gamma evaluation, normalized as T393.
W=mp.sqrt(6)/(32*mp.pi**3)*mp.gamma(mp.mpf(1)/24)*mp.gamma(mp.mpf(5)/24)*mp.gamma(mp.mpf(7)/24)*mp.gamma(mp.mpf(11)/24)
score=relative_agreement(6*derivative,W)
assert score>125
return {'six_times_derivative':mp.nstr(6*derivative,130),'watson_agreement':mp.nstr(score,8)}
def main():
report={'table':'T663','manifest':'agents/runs/20261101T235407Z-build-skills/manifest.json',
'skill_sha256':'a2b28d9953264d7b0ca62663785c2952190c5be6d6d3597fd695d9a54d606278',
'versions':{'python':platform.python_version(),'mpmath':mp.__version__},'controls':controls(),'rows':[]}
stored=None
if '--stored' in sys.argv:
stored=numberdb.table('T663')['Numbers']
assert list(stored)==[str(k) for k in ARGUMENTS]
for k in ARGUMENTS:
start=time.time()
a,b=compute(k,150),compute(k,210)
c,res,tail,tau=modular(k)
with mp.workdps(230):
agreement=relative_agreement(a,b); independent=relative_agreement(b,c)
assert agreement>125 and independent>125
row={'k':k,'value':b,'precision_agreement':mp.nstr(agreement,10),
'modular_agreement':mp.nstr(independent,10),'eta_residual':res,
'q_tail_bound':tail,'tau':tau}
negative=compute(-k,150)
assert relative_agreement(a,negative)>125
if stored is not None:
v=stored[str(k)]
text=v['number'] if isinstance(v,dict) else v
row['stored']=text
# Decimal convention is +/- one unit of the last printed place.
from decimal import Decimal
ulp=mp.mpf(10)**Decimal(text).as_tuple().exponent
assert abs(mp.mpf(text)-mp.mpf(c))<ulp
row['seconds']=round(time.time()-start,3)
report['rows'].append(row)
print('checked k=%d'%k,file=sys.stderr,flush=True)
report['derivative']=derivative_control()
# Check the exterior expansion at every selected exterior argument.
report['exterior_checks']=[]
for k in ARGUMENTS:
if k<=6:continue
c,tail=exterior_series(k)
with mp.workdps(160):
score=relative_agreement(c,report['rows'][ARGUMENTS.index(k)]['value'])
assert score>125
report['exterior_checks'].append({'k':k,'agreement':mp.nstr(score,8),'tail':tail})
print('series k=%d'%k,file=sys.stderr,flush=True)
# m(P_10)=3m(P_2)+2m(P_0), from the two L(g8,3) evaluations in BF section 1.
with mp.workdps(200):
byk={r['k']:mp.mpf(r.get('stored',r['value'])) for r in report['rows']}
residual=abs(byk[10]-3*byk[2]-2*byk[0])
assert residual<mp.mpf('1e-98' if stored else '1e-190')
report['arithmetic_relation_residual']=mp.nstr(residual,8)
print(json.dumps(report,indent=2))
if __name__=='__main__':main()