back to table · edit · history · where entries came from · files · download
7215 bytes, as of the version from 2026-10-30 19:31 (current). Recorded here, not run.
"""Independent integral and energy checks for numberdb.org/T628.
Run with SageMath: sage -python checks.py
Set NUMBERDB_CHECK_STORED=1 to check the live draft (key required) or call
main(document) with its exported JSON document. The polynomial/integral
checking method is adapted from the public T627 checks.py attachment.
The root and beta enclosures alone come from generate.py; polynomial
coefficients use factorials and oblate functions use validated quadrature.
The energy ratio is independently assembled from the surface potential,
equilibrium condition, inertia, and the virial energy relation in AKM.
"""
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 (MaclaurinEnergyRatios, enclose, root_enclosure,
polynomials, client_from_stdin)
PRECISION = 600
C = ComplexBallField(PRECISION)
F = RealBallField(PRECISION)
COARSE_BETA = ('0.35890','0.40345','0.42664','0.44084','0.45044',
'0.45736','0.46259','0.46667','0.46995','0.47265')
COARSE_E = ('0.98523','0.99375','0.99657','0.99784','0.99851',
'0.99891','0.99917','0.99935','0.99947','0.99957')
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,uncertain=False):
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)
tol = F(2)**(-400 if uncertain else -510)
value = explicit_p(n,x)*C.integral(integrand,0,1,abs_tol=tol,rel_tol=tol)
assert value.real().is_finite() and value.imag().is_finite()
assert F(value.real().rad()) < F(10)**(-115 if uncertain else -140)
assert value.imag().contains_zero()
return value.real()
def original_equation(k,x):
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 independent sign.')
def physical_energy(x,uncertain=False):
# Focal length a0=1, G=rho=1. AKM surface potential in Section 2.
q0=integral_q(0,C(x),uncertain)
q2=integral_q(2,C(x),uncertain)
a2=1+x*x
mass=4*F.pi()*x*a2/3
phi_pole=-mass*(q0+q2)
phi_equator=-mass*(q0-q2/2)
omega_squared=2*(phi_equator-phi_pole)/a2
inertia=2*mass*a2/5
kinetic=inertia*omega_squared/2
# AKM Section 4: W=(3/5)M Phi0+T; Phi0=Phi_pole.
gravitational=3*mass*phi_pole/5+kinetic
assert kinetic > 0 and gravitational < 0
beta=kinetic/(-gravitational)
assert beta.is_finite() and F(beta.rad()) < F(10)**(-110 if uncertain else -135)
return beta,omega_squared,kinetic,gravitational
def main(document=None):
for n,known in [(0,F.pi()/4),(2,F.pi()/2-F(3)/2)]:
got=integral_q(n,C(1))
assert got.overlaps(known) and not got.contains_zero()
control=physical_energy(F(1))[0]
assert control.overlaps(2-6/F.pi()) and not control.contains_zero()
print('CONTROLS: q0(1), q2(1), and beta(xi=1)=2-6/pi certified.',flush=True)
stored=document is not None or os.environ.get('NUMBERDB_CHECK_STORED')=='1'
if document is None and stored:
document=numberdb.table('T628',client=client_from_stdin())
if stored:
records=document['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':'T628','stored_readback':stored,
'skill_manifest':'agents/runs/20261030T190951Z-build-skills/manifest.json',
'skill_sha256':'dedecbdca15a728b6e594bfa79f0d149896fca336080ac27b53f513e0b16cd2e',
'sage_version':SAGE_VERSION,'python_version':platform.python_version(),
'root_working_bits':[500,700],'quadrature_bits':PRECISION,
'controls_passed':True,'entries':[]}
previous=F(0)
for k in range(1,11):
beta,lo,hi=enclose(k)
other,lo2,hi2=enclose(k,working_bits=700)
assert (lo,hi)==(lo2,hi2) and beta.overlaps(other)
assert previous < beta
previous=beta
x=F((lo+hi)/2).add_error(F((hi-lo)/2))
p,_,_,_,derivative=polynomials(k)
assert list(explicit_coefficients(2*k+2))==list(p.list()[::2])
f_lo=original_equation(k,F(lo));f_hi=original_equation(k,F(hi))
assert sign(f_lo)==1 and sign(f_hi)==-1
independent,omega2,kinetic,gravitational=physical_energy(x,uncertain=True)
assert beta.overlaps(independent)
center=decimal_rational(texts[k]);ulp=QQ(1)/10**100
assert len(texts[k].split('.')[1])==100
assert beta>F(center-ulp) and beta<F(center+ulp)
assert independent>F(center-ulp) and independent<F(center+ulp)
e=1/(1+x*x).sqrt();s=(1-e*e).sqrt();theta=e.arcsin()
beta_e=F(3)/(2*e*e)-1-3*s/(2*e*theta)
omega2_e=2*F.pi()*s*((3-2*e*e)*theta-3*e*s)/(e**3)
assert beta_e.overlaps(beta) and omega2_e.overlaps(omega2)
mass=4*F.pi()*x*(1+x*x)/3
w_e=-3*mass*mass*theta/(5*(1+x*x).sqrt()*e)
assert w_e.overlaps(gravitational)
for got,coarse in [(beta,COARSE_BETA[k-1]),(e,COARSE_E[k-1])]:
assert abs(got-F(decimal_rational(coarse)))<F(QQ(1)/200000)
evidence['entries'].append({'k':k,'beta':texts[k],
'xi_lower':str(lo),'xi_upper':str(hi),'beta_radius':str(beta.rad()),
'independent_beta':str(independent),'independent_beta_radius':str(independent.rad()),
'independent_root_endpoint_residuals':[str(f_lo),str(f_hi)],
'root_derivative_certificate_degree':int(derivative.degree()),
'exact_certificate_passed':True,'same_bracket_at_700_bits':True,
'all_displayed_identities_checked':True,'coarse_beta_and_e_match':True,
'independent_enclosure_inside_written_interval':True})
print('PASS',k,texts[k],flush=True)
print('EVIDENCE_JSON='+json.dumps(evidence,sort_keys=True),flush=True)
if __name__=='__main__':main()