back to table · edit · history · where entries came from · files · download
5984 bytes, as of the version from 2026-10-30 14:29. Recorded here, not run.
"""Independent checks for numberdb.org/T627.
Run: sage -python checks.py
Set NUMBERDB_CHECK_STORED=1 to read the draft/public table instead of local
values (an unpublished draft requires NUMBERDB_API_KEY). This script shares
only the generated values with generate.py: p_n is computed from factorial
coefficients, q_n from its defining integral using certified Arb quadrature.
The original AKM equation is evaluated at both endpoints of every written
100-digit eccentricity. The two endpoint signs must differ.
"""
import json
import os
from math import factorial
import platform
from functools import lru_cache
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.rational_field import QQ
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_arb import RealBallField
from sage.env import SAGE_VERSION
from generate import MaclaurinEccentricities, enclose, polynomials, client_from_stdin
PRECISION = 500
C = ComplexBallField(PRECISION)
F = RealBallField(PRECISION)
COARSE_E = ('0.98523','0.99375','0.99657','0.99784','0.99851',
'0.99891','0.99917','0.99935','0.99947','0.99957')
COARSE_AXIS = ('0.17126','0.11160','0.082750','0.065744','0.054534',
'0.046589','0.040664','0.036075','0.032417','0.029432')
def decimal_rational(text):
if 'e' in text.lower():
mant, power = text.lower().split('e')
return decimal_rational(mant)*QQ(10)**int(power)
a,b = text.split('.')
return QQ(int(a+b))/QQ(10)**len(b)
@lru_cache(None)
def explicit_coefficients(n):
assert n % 2 == 0
return tuple(QQ(factorial(n+2*j))/(2**n*factorial(n//2-j)*
factorial(n//2+j)*factorial(2*j)) for j in range(n//2+1))
def explicit_p(n,x):
value = x.parent()(0)
for c in explicit_coefficients(n)[::-1]:
value = value*x*x+c
return value
def integral_q(n,x):
# Substitute t=x/u in the integral over [x,infinity].
coefficients = [a*x**(2*j) for j,a in enumerate(explicit_coefficients(n))]
def integrand(u,analytic):
p_scaled = C(0)
for a in coefficients:
p_scaled = p_scaled*u*u+a
return x*u**(2*n)/((u*u+x*x)*p_scaled*p_scaled)
value = explicit_p(n,x)*C.integral(integrand,0,1)
assert value.real().is_finite() and value.imag().is_finite()
assert F(value.real().rad()) < F(10)**(-130)
assert value.imag().contains_zero()
return value.real()
def original_equation(k,e):
x = ((1-e*e)/(e*e)).sqrt()
xc = C(x)
return (2*explicit_p(2*k+2,x)*integral_q(2*k+2,xc)
-x/explicit_p(2,x)*(1-2*x*integral_q(2,xc)))
def sign(v):
assert v.is_finite()
if v > 0: return 1
if v < 0: return -1
raise ArithmeticError('Unresolved sign in independent integral check.')
def main(document=None):
# Controls must succeed before results are believed.
controls = [(0,F.pi()/4),(2,F.pi()/2-F(3)/2)]
for n,known in controls:
got = integral_q(n,C(1))
assert got.overlaps(known)
assert not got.contains_zero()
print('CONTROL: q0(1)=pi/4 and q2(1)=pi/2-3/2 certified.', flush=True)
stored = os.environ.get('NUMBERDB_CHECK_STORED') == '1'
table = document if document is not None else (
numberdb.table('T627',client=client_from_stdin()) if stored else None)
if document is not None:
assert stored, 'An exported document is a stored-value check.'
if stored:
records=table['Numbers']
if isinstance(records, dict):
assert [str(k) for k in records]==[str(k) for k in range(1,11)]
texts={int(k):(v['number'] if isinstance(v,dict) else v)
for k,v in records.items()}
else:
assert [str(r['params']['k']) for r in records]==[str(k) for k in range(1,11)]
texts={int(r['params']['k']):r['number'] for r in records}
else:
texts={k:to_text(enclose(k)[0],100) for k in range(1,11)}
evidence={'table':'T627','stored_readback':stored,
'skill_manifest':'agents/runs/20261030T141911Z-build-skills/manifest.json',
'skill_sha256':'28ecbeeccb0c5e5d6eb541b5a5b0f88aec47b154374d909085116d34f5629199',
'sage_version':SAGE_VERSION,'python_version':platform.python_version(),
'working_bits':500,'second_working_bits':700,'entries':[]}
previous = F(0)
for k in range(1,11):
e,lo,hi=enclose(k)
other,lo2,hi2=enclose(k,working_bits=700)
assert (lo,hi)==(lo2,hi2) and e.overlaps(other)
assert previous < e
previous=e
p,_,_,_,derivative=polynomials(k)
assert list(explicit_coefficients(2*k+2)) == list(p.list()[::2])
center=decimal_rational(texts[k]); ulp=QQ(1)/10**100
# Each eccentricity lies in (0.1,1), so 100 digits means 100 decimals.
assert len(texts[k].split('.')[1])==100
assert e > F(center-ulp) and e < F(center+ulp)
minus=original_equation(k,F(center-ulp))
plus=original_equation(k,F(center+ulp))
assert sign(minus)==-1 and sign(plus)==1
coarse=decimal_rational(COARSE_E[k-1]); coarse_axis=decimal_rational(COARSE_AXIS[k-1])
assert abs(e-F(coarse)) < F(QQ(1)/200000)
axis=(1-e*e).sqrt()
halfunit=QQ(1)/(2*10**len(COARSE_AXIS[k-1].split('.')[1]))
assert abs(axis-F(coarse_axis)) < F(halfunit)
row={'k':k,'eccentricity':texts[k], 'xi_lower':str(lo),'xi_upper':str(hi),
'eccentricity_radius':str(e.rad()),'independent_endpoint_residuals':[str(minus),str(plus)],
'global_uniqueness_certificate_degree':int(derivative.degree()),
'positive_coefficient_pattern_verified':True,
'coarse_eccentricity_and_axis_ratio_match':True,'same_bracket_at_700_bits':True}
evidence['entries'].append(row)
print('PASS',k,texts[k],flush=True)
print('EVIDENCE_JSON='+json.dumps(evidence,sort_keys=True),flush=True)
if __name__=='__main__':
main()