back to table · edit · history · where entries came from · files · download
3613 bytes, as of the version from 2026-10-30 06:50 (current). Recorded here, not run.
"""Independent radial-coordinate check of T613, without importing its solver.
Run with SageMath: sage -python independent_checks.py
Uses theta and -xi² theta' as functions of xi, and stops while theta>0.
At theta<1e-20, continue the homogeneous exterior solution to its zero.
The omitted density changes this zero by O(theta**(n+2)); the approximation
is heuristic. Positive controls use exact n=0,1 solutions before other rows.
Input evidence.json is optional and is used only for comparisons, never IVP
initial conditions or step control. This is build validation, not a critique.
"""
import numberdb.sage as numberdb
from sage.rings.real_mpfr import RealField
from sage.rings.rational_field import QQ
from fractions import Fraction
from pathlib import Path
import time,json
def evalpoly(coeff,x):
ans=x.parent()(0)
for c in reversed(coeff):ans=ans*x+c
return ans
def initial(index,F):
# Exact rational central coefficients, computed independently.
a=[Fraction(1)];power=[Fraction(1)];n=Fraction(index)
for k in range(200):
a.append(-power[k]/((2*k+2)*(2*k+3)))
j=k+1
power.append(sum((n*i-(j-i))*a[i]*power[j-i] for i in range(1,j+1))/j)
r=F(1)/4;u=r*r
af=[F(QQ(c.numerator)/QQ(c.denominator)) for c in a]
y=evalpoly(af,u)
q=-r*evalpoly([2*k*v for k,v in enumerate(af)],u)
return r,y,q
def radial(index,dps=70,degree=80):
F=RealField(numberdb.bits(dps));n=F(QQ(index));r,y,q=initial(index,F)
tol=F(10)**(-dps+20);stop=F(10)**-20;steps=0
while y>stop:
ys=[y];qs=[q];power=[y**n]
recip=[F((-1)**k*(k+1))/r**(k+2) for k in range(degree)]
for k in range(degree):
ys.append(-sum(qs[i]*recip[k-i] for i in range(k+1))/(k+1))
rhs=r*r*power[k]
if k>=1:rhs+=2*r*power[k-1]
if k>=2:rhs+=power[k-2]
qs.append(rhs/(k+1))
j=k+1
power.append(sum(((n+1)*i-j)*ys[i]*power[j-i] for i in range(1,j+1))/(j*y))
h=min(r/3, y*r*r/(4*q),F(1)/8)
for coeff,scale in [(ys,abs(y)),(qs,abs(q))]:
for k in range(degree-5,degree+1):
if coeff[k]:h=min(h,(tol*scale/abs(coeff[k]))**(F(1)/k)/2)
new_y=evalpoly(ys,h);new_q=evalpoly(qs,h)
assert 0<new_y<y and new_q>0
r+=h;y=new_y;q=new_q;steps+=1
if steps>6000:raise RuntimeError('too many radial steps')
# Exact zero of theta_ext(rho)=y+q(1/rho-1/r).
root=r/(1-r*y/q)
return root,{'steps':steps,'terminal_theta':str(y),'dps':dps,'degree':degree,'exterior_correction':str(root-r)}
def main():
evidence_path=Path(__file__).with_name('numerical-evidence.json')
expected=json.loads(evidence_path.read_text())['values'] if evidence_path.exists() else {}
output={}
for n in ['0','1','1/2','3/2','2','5/2','3','7/2','4','9/2']:
started=time.time();value,meta=radial(n);F=value.parent()
if n in ['0','1']:
control=F(6).sqrt() if n=='0' else F.pi()
err=abs(value-control)/control
assert err<F(10)**-35,(n,err)
meta['control_relative_error']=str(err)
assert abs(value*(1+F(10)**-25)-control)/control>F(10)**-35
if n in expected:
error=abs(value-F(expected[n]['high']))/value
assert error<F(10)**-35,(n,error)
meta['relative_difference_from_surface_solver']=str(error)
meta['value']=str(value);meta['seconds']=time.time()-started;output[n]=meta
print('CHECK',n,json.dumps(meta),flush=True)
print('DATA '+json.dumps(output),flush=True)
if __name__=='__main__':main()