back to table · edit · history · where entries came from · files · download
5248 bytes, as of the version from 2026-11-02 11:04 (current). Recorded here, not run.
"""Independent checks for numberdb.org/T658.
Run: sage -python checks.py (beside numerics.jsonl)
No generator code is imported. The full-range comparison integrates a direct
quadratic-root Jensen reduction of P_k as a polynomial in y, using x=e^(it).
The generator instead integrates a real trigonometric expression over a
half-period. Algebraic identities and PARI elliptic L-values supply checks
outside this shared Jensen principle. Exterior hypergeometric evaluations
use Rogers--Zudilin (7) only on its stated domain k outside [-4,2].
"""
import sys,json,time
from pathlib import Path
import numberdb.sage as numberdb
import mpmath as mp
from sage.libs.pari import pari
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.rational_field import QQ
mp.mp.dps=130
HERE=Path(__file__).resolve().parent
rows=[json.loads(x) for x in (HERE/'numerics.jsonl').read_text().splitlines()]
V={r['k']:mp.mpf(r['at190']) for r in rows if 'k' in r}
def argument(key):
if 'sqrt' not in key: return mp.mpf(key)
a={'4+2*sqrt(5)':(4,2),'4-2*sqrt(5)':(4,-2),
'8+4*sqrt(5)':(8,4),'8-4*sqrt(5)':(8,-4)}[key]
return a[0]+a[1]*mp.sqrt(5)
def compare(label,actual,expected,digits=110):
assert mp.isfinite(actual) and mp.isfinite(expected)
error=abs(actual-expected)/max(1,abs(expected))
assert error<mp.power(10,-digits),(label,mp.nstr(error,8))
print(json.dumps({'check':label,'relative_error':mp.nstr(error,8),'required_digits':digits}),flush=True)
# A known, nonzero control before the elliptic computations.
pari.set_real_precision(140)
lz=mp.mpf(str(pari('lfun(1,2)')))
compare('PARI lfun control zeta(2)',lz,mp.pi**2/6,120)
L3=(mp.zeta(2,mp.mpf(1)/3)-mp.zeta(2,mp.mpf(2)/3))/9
mu=3*mp.sqrt(3)*L3/(4*mp.pi)
compare('factorization k=-1',V['-1'],2*mu,120)
compare('spanning tree k=8',V['8'],5*mu,120)
assert V['0']==0
R=PolynomialRing(QQ,names=('x','y'));x,y=R.gens()
P=(1+x)*(1+y)*(x+y)
assert P+x*y==(x+y+1)*(x*y+x+y)
assert P==(1+x)*(1+y)*(x+y)
print(json.dumps({'check':'exact factorizations k=0,-1','passed':True}),flush=True)
# Three independent elliptic L-values, with conductors checked rather than assumed.
curves={14:[1,0,1,4,-6],20:[0,0,0,-108,297],36:[0,0,0,0,1]}
Ls={}
for N,ainvs in curves.items():
E=pari(ainvs).ellinit()
assert int(E.ellglobalred()[0])==N
Ls[N]=mp.mpf(str(pari.lfun(E,2,precision=470)))
print(json.dumps({'check':'elliptic L value','conductor':N,'ainvs':ainvs,'L2':mp.nstr(Ls[N],125)}),flush=True)
for key,c in [('1',mp.mpf(7)/2),('7',21),('-8',35)]:
compare('conductor 14 k='+key,V[key],c*Ls[14]/mp.pi**2,120)
compare('conductor 20 k=4',V['4'],10*Ls[20]/mp.pi**2,120)
compare('conductor 20 k=-2',V['-2'],15*Ls[20]/mp.pi**2,120)
compare('conductor 36 k=2',V['2'],9*Ls[36]/(2*mp.pi**2),120)
compare('positive quadratic identity',2*V['4+2*sqrt(5)']-V['8+4*sqrt(5)'],V['4'],120)
assert abs(V['1']-V['-1'])>mp.mpf('0.1')
# This control must fail: a perturbation larger than a stored unit is detected.
assert abs(V['4']+mp.mpf('1e-90')-10*Ls[20]/mp.pi**2)>mp.mpf('1e-100')
print(json.dumps({'check':'deliberate 1e-90 perturbation rejected','passed':True}),flush=True)
# Direct quadratic root integrator; no scalar expression from generate.py.
def roots_integral(key):
k=argument(key)
ends=[mp.mpf(0),mp.pi]
if -1<=k<0:
ends.extend(2*mp.acos((1+s*mp.sqrt(1+k))/2) for s in [-1,1])
if 0<k<8:
ends.append(2*mp.acos((mp.sqrt(1+k)-1)/2))
ends=sorted(set(ends))
def integrand(t):
x=mp.exp(mp.j*t)
a=1+x; b=(1+x)**2-k*x; c=x*(1+x)
disc=mp.sqrt(b*b-4*a*c)
r1=(-b+disc)/(2*a);r2=(-b-disc)/(2*a)
return mp.log(abs(a))+mp.log(max(1,abs(r1)))+mp.log(max(1,abs(r2)))
return mp.quad(integrand,ends,method='tanh-sinh',maxdegree=9)/mp.pi
with mp.workdps(125):
compare('quadratic integrator nonzero control k=-1',roots_integral('-1'),2*mu,115)
for key in V:
t=time.monotonic()
result=roots_integral(key)
compare('direct quadratic Jensen k='+key,result,V[key],110)
# An independently implemented hypergeometric continuation; real part is used
# only on the source's permitted domain. Boundary k=8 uses its factorization
# control instead of a slowly convergent hypergeometric series at z=1.
with mp.workdps(125):
def H(z):
if abs(z) < 1:
term=mp.mpf(1); total=term; n=0
while True:
term *= z*(n+mp.mpf(4)/3)*(n+mp.mpf(5)/3)*(n+1)/(n+2)**3
total += term
n += 1
if abs(term)/(1-abs(z)) < mp.mpf('1e-123'):return total
if n > 2000000:raise ArithmeticError('Hypergeometric series too slow')
return mp.hyper([mp.mpf(4)/3,mp.mpf(5)/3,1,1],[2,2,2],z)
for key in V:
k=argument(key)
if not (k < -4 or k > 2) or k==8:continue
result=mp.re(mp.log((k+4)*(k-2)**4/k**2)-2*k*k/(k+4)**3*H(27*k*k/(k+4)**3)-8*k/(k-2)**3*H(27*k/(k-2)**3))/3
compare('RZ equation 7 k='+key,result,V[key],110)
n32=mp.log(32)/3-H(mp.mpf(27)/32)/16
compare('Hessian normalization g(4)=3*N(32)/4',V['4'],3*n32/4,110)
print(json.dumps({'check':'all independent checks passed','entries':len(V)}),flush=True)