back to table · edit · history · where entries came from · files · download
9018 bytes, as of the version from 2026-11-09 16:19 (current). Recorded here, not run.
"""Independent checks for T693, runnable as: sage -python checks.py.
With NUMBERDB_STORED=1, read T693 and test the written intervals as well.
The separate formal-series computation uses reciprocal/exponential convolutions,
not the generator's Bernoulli routine or logarithmic-derivative recurrence.
The Kummer differential equation and direct Arb U computations are additional
checks external to the defining finite sum. See derivation.md for limitations.
"""
import json
import os
import sys
from math import factorial
from functools import reduce
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.real_mpfi import RealIntervalField
from sage.rings.complex_arb import ComplexBallField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.polynomial.laurent_polynomial_ring import LaurentPolynomialRing
from sage.env import SAGE_VERSION
from numberdb._write import to_text
from check import exactness, measure, agrees_with, names_its_rings
import generate
N = 60
S = PolynomialRing(QQ, 'z')
z = S.gen()
def independent_polynomials(limit):
# Invert (exp(u)-1)/u using exact convolution, avoiding bernoulli().
inv = [QQ.one()]
for k in range(1, limit+2):
inv.append(-sum(inv[k-r]/QQ(factorial(r+1)) for r in range(1,k+1)))
lam = [QQ.zero()]+inv[2:limit+2]
# Sum (-z lambda(u))^m/m! by truncated Cauchy products.
power = [QQ.one()]+[QQ.zero()]*limit
p = [S.one()]+[S.zero()]*limit
for m in range(1,limit+1):
power = [sum(power[k-r]*lam[r] for r in range(1,k+1)) for k in range(limit+1)]
for k in range(m,limit+1):
p[k] += power[k]*(-z)**m/QQ(factorial(m))
c=[]
R=PolynomialRing(QQ,'x'); x=R.gen()
for ell in range(limit+1):
v=S.zero()
for k in range(ell+1):
j=ell-k
# H(nu,j)=(-1)^j (1/2-nu)_j (1/2+nu)_j/j!.
poch=QQ.one()
for r in range(j): poch *= (QQ(1)/2-(k+1)+r)*(QQ(1)/2+(k+1)+r)
v += QQ((-1)**j)*poch/QQ(factorial(j)*4**j)*z**k*p[k]
assert all(v[r]==0 for r in range(1,v.degree()+1,2))
c.append(R([v[2*m]*(-4)**m for m in range(ell+1)]))
return c,inv
def written_bounds(text):
from decimal import Decimal
if '/' in text: return (QQ(text),QQ(text))
d=Decimal(text); n,dn=d.as_integer_ratio(); center=QQ(n)/QQ(dn)
if '.' not in text and 'e' not in text.lower(): return center,center
e=d.as_tuple().exponent
ulp=QQ(10)**e
return center-ulp,center+ulp
def encloses(text,v):
lo,hi=written_bounds(text)
if isinstance(v,type(QQ.one())):return lo<=v<=hi
return bool(v.is_finite() and lo < QQ(v.lower()) and QQ(v.upper()) < hi)
def kummer_residuals(polys):
L=LaurentPolynomialRing(QQ,'t');t=L.gen()
a=[sum(L(c)*QQ(-1)**m/QQ(4)**m*t**(4*m-ell) for m,c in enumerate(p.list())) for ell,p in enumerate(polys)]
return [4*a[l+1].derivative()-a[l].derivative(2)+(t*t+QQ(3)/4/t**2)*a[l] for l in range(len(a)-1)]
def run():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
numberdb.configure(api_key=sys.stdin.read().strip())
out={'table':'T693','sage':SAGE_VERSION,'numberdb':numberdb.__version__,
'skill_manifest':'agents/runs/20261109T154709Z-build-skills/manifest.json',
'skill_sha256':'886390324dd1971e381a43cf3b46b822ab0bc67d62f31d8481e29ec79556746b'}
ps=generate.coefficient_polynomials(N)
ips,inv=independent_polynomials(N)
assert inv[:4]==[QQ(1),-QQ(1)/2,QQ(1)/12,QQ(0)]
assert encloses('1.25',QQ(5)/4) and not encloses('1.25',QQ(3)/2)
assert not exactness(ps) and not exactness(ips)
assert not names_its_rings(generate.__file__)
assert not agrees_with(dict(enumerate(ps)),lambda l:ips[l])
bad=list(ps[:2]);bad[1]=bad[1]+1
assert any(kummer_residuals(bad))
assert all(r==0 for r in kummer_residuals(ps))
expected=[ [1], [QQ(3)/16,QQ(1)/3],[-QQ(15)/512,QQ(5)/16,QQ(1)/18],
[QQ(105)/8192,QQ(35)/512,QQ(23)/160,QQ(1)/162],
[-QQ(4725)/524288,-QQ(105)/8192,QQ(973)/5120,QQ(137)/4320,QQ(1)/1944],
[QQ(72765)/8388608,QQ(3465)/524288,QQ(2849)/16384,QQ(17365)/193536,QQ(229)/51840,QQ(1)/29160] ]
assert all(ps[i]==generate.R(p) for i,p in enumerate(expected))
const=QQ.one()
for ell,p in enumerate(ps):
if ell:const*=QQ(4-(2*ell-1)**2)/QQ(16*ell)
assert p.degree()==ell and p[ell]==QQ(1)/QQ(factorial(ell)*3**ell) and p[0]==const
R=RealBallField(700); pi2=R.pi()**2
main={i:generate.evaluate(p) for i,p in enumerate(ps)}
other={i:ips[i](pi2) for i in range(N+1)}
texts={i:to_text(v,digits=100) for i,v in main.items()}
for i in range(N+1):
assert other[i].is_finite() and (i==0 or other[i].accuracy()>650)
assert encloses(texts[i],main[i]) and (i==0 or encloses(texts[i],other[i]))
neg=[i for i in range(N+1) if main[i]<0]
assert neg[0]==50 and all(main[i]>0 for i in range(50))
cancellation=[]
for p,v in zip(ps,main.values()):
if p.degree()>0:
abs_sum=sum(abs(c)*pi2**k for k,c in enumerate(p.list()))
cancellation.append(float((abs_sum/abs(R(v))).log()/R(10).log()))
out.update(controls='reciprocal series, exact interval inclusion, corrupted coefficient and nonzero U identity passed',
entries=N+1,exact_polynomial_agreements=N+1,printed_formula_agreements=6,kummer_ode_residuals=60,
minimum_accuracy_bits=min(int(v.accuracy()) for i,v in main.items() if i),
working_bits=numberdb.bits(100,losing=generate.WORKING_GUARD),
minimum_independent_accuracy_bits=min(int(other[i].accuracy()) for i in range(1,N+1)),
max_cancellation_decimal_digits=max(cancellation),negative_orders=neg,
raw_ball_measurement=measure(main),written_measurement=measure(texts),
values={str(k):v for k,v in texts.items()},polynomials={str(i):str(p) for i,p in enumerate(ps)})
# Corpus ingredients read through numberdb, never altered.
bern=numberdb.table('T13')['Numbers']
hank=numberdb.table('T569')['Numbers']
def entries(block):
if isinstance(block,list):
return [(next(iter(r['params'].values())),r['number']) for r in block]
return [(k,v['number'] if isinstance(v,dict) else v) for k,v in block.items()]
bc=0
for k,v in entries(bern):
k=int(k)
if 0<=k<=N+1:
assert QQ(str(v))==inv[k]*QQ(factorial(k));bc+=1
hc=0
Hring=PolynomialRing(QQ,'v');v=Hring.gen()
for j,expr in entries(hank):
j=int(j)
if j>N:continue
poly=Hring(str(expr))
for k in range(N-j+1):
h=QQ.one()
for r in range(1,j+1):h*=QQ(4*(k+1)**2-(2*r-1)**2)/QQ(4*r)
assert 2**j*poly(k+1)==h;hc+=1
out['bernoulli_ingredient_checks']=bc;out['hankel_ingredient_checks']=hc
# Direct complex Tricomi values validate amplitude, phase, and truncation.
C=ComplexBallField(1200);rb=RealBallField(1200);zz=2*C.pi()*C.gen(0)
control=C(2).hypergeometric_U(1,2)
assert control.real().is_finite() and control.imag().is_finite() and (QQ(1)/2 in control.real()) and (0 in control.imag())
assert control.real().rad() < QQ(2)**-1000 and control.imag().rad() < QQ(2)**-1000
out['tricomi_control']={'real':str(control.real()),'imaginary':str(control.imag())}
used=main
if os.environ.get('NUMBERDB_STORED')=='1':
pairs=entries(numberdb.table('T693')['Numbers'])
assert [int(k) for k,v in pairs]==list(range(N+1))
used={}
for k,text in pairs:
k=int(k); assert encloses(str(text),main[k])
if k: assert encloses(str(text),other[k])
lo,hi=written_bounds(str(text))
used[k]=rb(RealIntervalField(1200)(lo,hi))
assert used[k].is_finite() and (k==0 or used[k].accuracy()>320)
assert [i for i in range(N+1) if used[i]<0]==neg
assert str(pairs[0][1])=='1'
out['stored_interval_checks']=len(pairs)
out['stored_sign_checks']=len(pairs)
tests=[]
for n in [25,100,400,1600]:
val=(-4*C(factorial(n-1))*zz.hypergeometric_U(n,0)).real()
assert val.is_finite() and val.accuracy()>500
r=rb(n);pi=rb.pi();phase=2*(pi*r).sqrt()
amplitude=rb(2)**(QQ(9)/4)*pi**(QQ(3)/4)*r**(-QQ(3)/4)*(-phase).exp()
target=val/amplitude
for L in [0,1,3,5,10,20,40,60]:
approx=sum(rb(used[l])*(2*pi*r)**(-QQ(l)/2)*(phase+pi*l/4+3*pi/8).sin() for l in range(L+1))
error=abs(target-approx)
tests.append({'n':n,'L':L,'normalized_error':str(error),'error_mid':float(error.mid()),'U_accuracy_bits':int(val.accuracy())})
out['tricomi_tests']=tests
for n in [25,100,400,1600]:
errors={row['L']:row['error_mid'] for row in tests if row['n']==n}
assert errors[5]<errors[0] and errors[20]<errors[5]
return out
if __name__=='__main__':
print(json.dumps(run(),indent=2))