back to table · edit · history · where entries came from · files · download
5017 bytes, as of the version from 2026-10-30 02:35 (current). Recorded here, not run.
"""Independent radial-coordinate checks for T615.
Run: sage -python independent_checks.py, with numerical-evidence.json beside it.
Adapted from the live T614 independent_checks.py; imports no production solver.
Integrates theta, enclosed mass, and the pressure integral with respect to xi,
stopping at positive theta<1e-30. A leading mass-tail correction and exterior
root estimate give independent concentration values. The virial relation tests
the integrated density profile. Endpoint error remains heuristic.
This is build validation, not the separate campaign 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)**-30;steps=0
# Integrate the central series to initialize J independently.
# J is initialized by high-precision quadrature of the analytic local IVP.
import mpmath as mp
with mp.workdps(dps):
nn=mp.mpf(str(n));rr=mp.mpf(str(r))
aa=[mp.mpf(1)];bb=[mp.mpf(1)]
for k in range(200):
aa.append(-bb[k]/((2*k+2)*(2*k+3)))
j=k+1
bb.append(sum(((nn+1)*i-j)*aa[i]*bb[j-i] for i in range(1,j+1))/j)
def central(z):
val=mp.polyval(list(reversed(aa)),z*z)
return z*z*val**(nn+1)
J=F(str(mp.quad(central,[0,rr])))
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))
jcoeff=[]
yp=[sum(ys[i]*power[k-i] for i in range(k+1)) for k in range(degree)]
for k in range(degree):
jcoeff.append((r*r*yp[k]+(2*r*yp[k-1] if k else 0)+(yp[k-2] if k>=2 else 0))/(k+1))
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
J+=h*evalpoly(jcoeff,h)
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)
tail=r**4*y**(n+1)/(q*(n+1))
mass=q+tail
J+=r**4*y**(n+2)/(q*(n+2))
virial=(n+1)*mass**2/((5-n)*root)
err=abs(J-virial)/abs(virial)
assert err<F(10)**-40, ('virial',index,err)
return mass,{'radius':str(root),'steps':steps,'terminal_theta':str(y),'dps':dps,'degree':degree,'mass_tail':str(tail),'pressure_integral':str(J),'virial_relative_error':str(err)}
def main():
expected=json.loads(Path(__file__).with_name('numerical-evidence.json').read_text())['values']
output={'table':'T615','skill_manifest':'agents/runs/20261030T020254Z-build-skills/manifest.json','values':{}}
for n in ['0','1','1/2','3/2','2','5/2','3','7/2','4','9/2']:
started=time.time();mass,meta=radial(n);F=mass.parent();x=F(meta['radius'])
ratio=x**3/(3*mass)
if n in ['0','1']:
control=F(1) if n=='0' else F.pi()**2/3
err=abs(ratio-control)/control
assert err<F(10)**-40,(n,err)
meta['control_relative_error']=str(err)
assert abs(ratio*(1+F(10)**-25)-control)/control>F(10)**-40
error=abs(ratio-F(expected[n]['210']['ratio']))/ratio
assert error<F(10)**-40,(n,error)
meta['relative_difference_from_surface_solver']=str(error)
meta['mass']=str(mass);meta['ratio']=str(ratio);meta['seconds']=time.time()-started
if n in ['3/2','3']:
reference=F('5.99') if n=='3/2' else F('54.18')
assert abs(ratio-reference)<F('0.005')
meta['princeton_page_3_rounding_check']=True
output['values'][n]=meta
print('CHECK',n,json.dumps(meta),flush=True)
print('DATA '+json.dumps(output),flush=True)
if __name__=='__main__':main()