back to table · edit · history · where entries came from · files · download
5456 bytes, as of the version from 2026-10-16 06:19 (current). Recorded here, not run.
"""Independent checks for Student t quantiles -- numberdb.org/T454.
Run with SageMath (after installing numberdb):
sage -python repair-20261016T061049Z-checks.py
Keep repair-20261016T061049Z-before.json and repair-20261016T061049Z-nist.txt beside this file.
The checks write JSON to stdout and never write to NumberDB.
The CDF calculation substitutes x = sqrt(nu)*tan(theta) in the density.
For m=nu-1, integrate cos(theta)^m by the recurrence
I_m = sin(theta)*cos(theta)^(m-1)/m + (m-1)*I_(m-2)/m,
with I_0=theta, I_1=sin(theta), and normalize by 2*I_m(pi/2).
This uses no incomplete-beta, gamma or hypergeometric routine. All CDF
comparisons use Arb balls at 768 bits. Decimal endpoints are exact rationals.
"""
import json
import hashlib
import re
from pathlib import Path
from fractions import Fraction
from decimal import Decimal, localcontext
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
HERE = Path(__file__).resolve().parent
raw = (HERE / 'repair-20261016T061049Z-before.json').read_bytes()
doc = json.loads(raw)
rows = doc['Numbers']
B = RealBallField(768)
def rational(text):
f = Fraction(str(text))
return QQ(f.numerator) / QQ(f.denominator)
def number(entry):
return str(entry['number'] if isinstance(entry, dict) else entry)
def cdf(nu, x):
x = B(x)
root_nu = B(nu).sqrt()
theta = (x/root_nu).arctan()
den = (B(nu) + x*x).sqrt()
sine, cosine = x/den, root_nu/den
m = nu - 1
if m % 2 == 0:
integral, whole, start = theta, B.pi()/2, 2
else:
integral, whole, start = sine, B(1), 3
for k in range(start, m+1, 2):
integral = sine*cosine**(k-1)/k + B(k-1)/k*integral
whole = B(k-1)/k*whole
result = B(1)/2 + integral/(2*whole)
assert result.is_finite()
assert B(result.rad()) < B(2)**(-700)
return result
# Controls: symmetry center, Cauchy quartiles, and the nu=2 quartiles.
controls = [(1, QQ(1), QQ(3)/4), (1, QQ(-1), QQ(1)/4),
(2, (B(2)/3).sqrt(), QQ(3)/4),
(2, -(B(2)/3).sqrt(), QQ(1)/4), (100, QQ(0), QQ(1)/2)]
for nu, x, p in controls:
assert B(p) in cdf(nu, x), (nu, x, p)
nus = list(range(1,31)) + [40,50,60,100]
nist_probabilities = ['9/10','19/20','39/40','99/100','199/200','999/1000']
negative_probabilities = ['1/10','1/20','1/40','1/100','1/200']
probabilities = ['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(rows) == [str(nu) for nu in nus]
assert {1-rational(p) for p in nist_probabilities[:-1]} == {rational(p) for p in negative_probabilities}
assert rational('1/1000') not in {rational(p) for p in probabilities}
exact_count = decimal_count = symmetry_pairs = negative_count = 0
brackets = []
for nu in nus:
block = rows[str(nu)]
assert list(block) == probabilities
for p, entry in block.items():
q = rational(p)
text = number(entry)
value = rational(text)
if p in negative_probabilities:
assert value < 0
negative_count += 1
if q < QQ(1)/2:
other = rational(number(block[str(1-q)]))
assert value == -other
symmetry_pairs += 1
if '.' not in text and 'e' not in text.lower():
assert (q == QQ(1)/2 and value == 0) or (nu == 1 and (q,value) in [(QQ(1)/4,-1),(QQ(3)/4,1)])
exact_count += 1
else:
places = Decimal(text).as_tuple().exponent
ulp = QQ(10)**places
lo = cdf(nu, value-ulp)
hi = cdf(nu, value+ulp)
assert lo < B(q), (nu,p,'lower')
assert hi > B(q), (nu,p,'upper')
assert len(Decimal(text).as_tuple().digits) == 100
decimal_count += 1
brackets.append({'nu':nu,'p':p,'lower_cdf_less':True,'upper_cdf_greater':True})
nist_text = (HERE/'repair-20261016T061049Z-nist.txt').read_text()
header = re.search(r'Probability less.*?(0\.90\s+0\.95\s+0\.975\s+0\.99\s+0\.995\s+0\.999)',nist_text,re.S)
assert header
assert [rational(s) for s in header.group(1).split()] == [rational(p) for p in nist_probabilities]
parsed = re.findall(r'(?<![\d.])(\d+)\.\s+((?:\d+\.\d{3}\s+){5}\d+\.\d{3})',nist_text)
# The HTML repeats a '100.' for normal values; retain the first row of each degree.
nist = {}
for key, vals in parsed:
nist.setdefault(int(key), vals.split())
nist_count = 0
with localcontext() as ctx:
ctx.prec = 150
for nu in nus:
for p, expected in zip(nist_probabilities,nist[nu]):
assert Decimal(number(rows[str(nu)][p])).quantize(Decimal('0.001')) == Decimal(expected)
nist_count += 1
print(json.dumps({'tid':'T454','document_sha256':hashlib.sha256(raw).hexdigest(),
'degrees':len(nus),'probabilities':len(probabilities),'entries':exact_count+decimal_count,
'exact_integers':exact_count,'certified_decimal_brackets':decimal_count,
'working_bits':768,'cdf_radius_bound_bits':700,'controls':len(controls),
'symmetry_pairs':symmetry_pairs,'five_complements_negative_rows':negative_count,
'nist_shared_entries_rounded_to_3dp':nist_count,'nist_header':header.group(1).split(),
'unpaired_probability':'999/1000','omitted_complement':'1/1000',
'all_checks_passed':True,'brackets':brackets},indent=2))