back to table · edit · history · where entries came from · files · download
4077 bytes, as of the version from 2026-10-28 22:35 (current). Recorded here, not run.
"""Read-only checks for T457; run: sage -python check_T457.py T457-before.json.
Independent closed forms are enclosed using 600-bit Arb arithmetic. Numerical
CDF bracket checks use mpmath at 150 decimal digits; these cross-checks do not
certify all of the generator's remainder bounds. No values are generated or
submitted by this script. In the campaign use agents/sage.sh with both files.
"""
import json, sys
from pathlib import Path
from fractions import Fraction
import numberdb.sage as numberdb
from sage.rings.real_arb import RealBallField
from sage.rings.rational_field import QQ
import mpmath as mp
path = Path(sys.argv[1]) if len(sys.argv) > 1 else Path(__file__).with_name('T457-before.json')
doc = json.loads(path.read_text())
B = RealBallField(600)
mp.mp.dps = 150
nus = list(range(1,11)) + [20,30,60,120]
ps = ['1/200','1/100','1/40','1/20','1/10','1/4','1/2','3/4','9/10','19/20','39/40','99/100','199/200','999/1000']
assert list(doc['Parameters']) == ['nu','k','p']
rows = {}
for nu, by_k in doc['Numbers'].items():
for k, by_p in by_k.items():
for p, entry in by_p.items():
spelling = entry['number'] if isinstance(entry, dict) else entry
prob = Fraction(p)
assert str(prob) == p and 0 < prob < 1
center = Fraction(spelling)
fraction_digits = len(spelling.split('.')[1])
ulp = Fraction(1, 10**fraction_digits)
rows[int(nu),int(k),p] = (center,ulp)
assert set(rows) == {(nu,k,p) for nu in nus for k in [2,3] for p in ps}
assert len(rows) == 392
assert Fraction(1)-Fraction(19,20) == Fraction(5,100)
assert Fraction(1)-Fraction(99,100) == Fraction(1,100)
def arb(r):
return B(QQ(r.numerator)/QQ(r.denominator))
def mpq(r):
return mp.mpf(r.numerator)/r.denominator
closed_count = 0
max_closed_error = mp.mpf(0)
for (nu,k,p), (c,u) in rows.items():
pball = arb(Fraction(p))
if nu == 1 and k == 2:
v = B(2).sqrt()*(B.pi()*pball/2).tan()
elif nu == 1 and k == 3:
v = (1/(2*(B.pi()*(1-pball)/6).sin()**2)-2).sqrt()
elif nu == 2 and k == 2:
v = 2*pball/(1-pball**2).sqrt()
else:
continue
assert v.is_finite()
assert v > arb(c-u) and v < arb(c+u), (nu,k,p)
closed_count += 1
assert closed_count == 42
print('600-bit Arb closed-form enclosures: 42/42 inside stored intervals', flush=True)
def cdf2(x, nu):
return mp.betainc(mp.mpf(1)/2,mp.mpf(nu)/2,0,x*x/(2*nu+x*x),regularized=True)
def cdf3(x, nu):
return 1-6/mp.pi*mp.quad(lambda theta:(1+x*x/(2*nu*mp.cos(theta)**2))**(-mp.mpf(nu)/2),[0,mp.pi/6])
controls = [abs(cdf2(mp.sqrt(2),1)-mp.mpf(1)/2),
abs(cdf2(2/mp.sqrt(3),2)-mp.mpf(1)/2),
abs(cdf3(mp.sqrt(2),1)-(1-6/mp.pi*mp.asin(1/mp.sqrt(8)))),
abs(cdf3(mp.sqrt(2),2)-6/(mp.pi*mp.sqrt(3))*mp.atan(mp.mpf(1)/3))]
assert max(controls) < mp.mpf('1e-140')
checked = []
for (nu,k,p),(c,u) in rows.items():
if k == 3 and not (nu in [1,2,10,120] and p in ['1/200','1/4','1/2','19/20','999/1000']):
continue
f = cdf2 if k == 2 else cdf3
prob = mpq(Fraction(p))
low=f(mpq(c-u),nu)-prob
high=f(mpq(c+u),nu)-prob
assert mp.isfinite(low) and mp.isfinite(high)
assert low < 0 < high, (nu,k,p,mp.nstr(low),mp.nstr(high))
checked.append({'nu':nu,'k':k,'p':p,'lower_cdf_minus_p':mp.nstr(low,12),'upper_cdf_minus_p':mp.nstr(high,12)})
print('150-digit numerical endpoint checks: %d/%d passed' % (len(checked),len(checked)), flush=True)
result = {'table':'T457','entries':len(rows),'selection_matches':True,'exact_rational_probabilities':True,
'upper_tail_examples_exact':True,'arb_precision_bits':600,'closed_form_enclosures':closed_count,
'mpmath_decimal_digits':150,'known_answer_control_count':4,'maximum_control_residual':mp.nstr(max(controls),15),
'numerical_brackets':checked,'limitations':'The 216 mpmath endpoint checks are numerical cross-checks, not a fresh proof of all Arb remainder bounds. No values were modified.'}
print('RESULT_JSON='+json.dumps(result, separators=(',',':')))