back to table · edit · history · where entries came from · files · download
4947 bytes, as of the version from 2026-10-16 06:58 (current). Recorded here, not run.
"""Independent repair checks for Quantiles of the F-distribution, numberdb.org/T455.
Run: sage -python repair_checks.py [document.json] [nist-overlap.json]
Default fixture names are T455-document-final.json and nist-overlap.json.
No API credential is needed. Install numberdb and mpmath first.
The NIST comparison uses exact rational arithmetic, including the entire
last-place interval represented by each stored decimal. mpmath at 140
decimal places independently checks strict CDF signs at both ends of all
stored intervals; this fixed-precision check does not recertify ball bounds.
Exact rows are checked by rational identities, without inferring exactness
from decimal digits. This file imports no code from the table generator.
"""
import json
import sys
from decimal import Decimal
from fractions import Fraction
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.real_mpfr import RealField
from sage.rings.rational_field import QQ
import mpmath as mp
root = Path(__file__).resolve().parent
doc_path = Path(sys.argv[1]) if len(sys.argv) > 1 else root / 'T455-document-final.json'
nist_path = Path(sys.argv[2]) if len(sys.argv) > 2 else root / 'nist-overlap.json'
doc = json.loads(doc_path.read_text())
nist = json.loads(nist_path.read_text())
mp.mp.dps = 140
def interval(text):
q = Fraction(text)
if '.' not in text and 'e' not in text.lower():
return q, q
unit = Fraction(10) ** Decimal(text).as_tuple().exponent
return q - unit, q + unit
def m(q):
q = Fraction(q)
return mp.mpf(q.numerator) / q.denominator
def cdf(x, d1, d2):
x = m(x)
return mp.betainc(mp.mpf(d1)/2, mp.mpf(d2)/2, 0,
d1*x/(d1*x+d2), regularized=True)
# Known-answer controls, including signs that must be nonzero.
assert abs(cdf(Fraction(19), 2, 2) - m(Fraction(19, 20))) < mp.mpf('1e-135')
assert abs(cdf(Fraction(1, 3), 1, 1) - m(Fraction(1, 3))) < mp.mpf('1e-135')
assert cdf(Fraction(18), 2, 2) < m(Fraction(19, 20)) < cdf(Fraction(20), 2, 2)
p_sage = QQ('99/100')
exact_example = ((QQ(1)-p_sage)**(-2)-QQ(1))/QQ(2)
assert exact_example == QQ('9999/2')
print('Controls passed; Sage exact F(2,1,99/100) =', exact_example, flush=True)
rows = {}
for d1, by_d2 in doc['Numbers'].items():
for d2, by_p in by_d2.items():
for p, entry in by_p.items():
text = entry['number'] if isinstance(entry, dict) else entry
rows[(int(d1), int(d2), p)] = (text, *interval(text))
probabilities = ('1/100', '1/40', '1/20', '1/2', '19/20', '39/40', '99/100')
assert set(rows) == {(a,b,p) for a in range(1,11) for b in range(1,11) for p in probabilities}
exact_count = 0
decimal_count = 0
for (d1, d2, p), (text, lo, hi) in rows.items():
prob = Fraction(p)
if lo == hi:
assert lo > 0
assert ((d1 == d2 and prob == Fraction(1,2) and lo == 1)
or (d1 == 2 and (1 + 2*lo/d2)**d2 == (1-prob)**-2)
or (d2 == 2 and (d1*lo/(d1*lo+2))**d1 == prob**2))
exact_count += 1
else:
left = cdf(lo, d1, d2)-m(prob)
right = cdf(hi, d1, d2)-m(prob)
assert mp.isfinite(left) and mp.isfinite(right) and left < 0 < right, (d1,d2,p)
decimal_count += 1
assert (exact_count, decimal_count) == (32, 668)
print('All 700 entries checked: 32 exact identities, 668 last-place brackets.', flush=True)
report = {'table': 'T455', 'entry_count': len(rows), 'exact_checks': exact_count,
'decimal_brackets_checked': decimal_count, 'cdf_precision_decimal_digits': mp.mp.dps,
'cdf_method': 'mpmath incomplete beta; fixed precision, independent numerical evidence',
'sage_exact_F_2_1_99_100': str(exact_example), 'nist_comparisons': {}}
for p in ('19/20', '99/100'):
comparisons = []
for row in nist:
if row['p'] != p:
continue
text, lo, hi = rows[(row['d1'], row['d2'], p)]
printed = Fraction(row['printed'])
unit = Fraction(10) ** Decimal(row['printed']).as_tuple().exponent
worst = max(abs(lo-printed), abs(hi-printed))/unit
comparisons.append({**row, 'stored': text,
'maximum_distance_in_printed_units': mp.nstr(m(worst), 25),
'within_one_printed_unit': worst <= 1})
bad = [r for r in comparisons if not r['within_one_printed_unit']]
assert len(comparisons) == 100
assert len(bad) == (0 if p == '19/20' else 7)
report['nist_comparisons'][p] = {'count': 100, 'outside_one_unit': len(bad),
'comparisons': comparisons}
print('NIST p =', p, ':', len(bad), 'of 100 outside one printed unit', flush=True)
assert Fraction(rows[(2,1,'99/100')][0]) == Fraction(9999,2)
example = next(r for r in nist if (r['d1'],r['d2'],r['p']) == (2,1,'99/100'))
assert (Fraction(example['printed'])-Fraction(9999,2))/Fraction(1,100) == 2
print('CHECK_RESULT_JSON=' + json.dumps(report, sort_keys=True))