back to table · edit · history · where entries came from · files · download
4230 bytes, as of the version from 2026-10-16 07:34 (current). Recorded here, not run.
"""Quantiles of the Kolmogorov distribution -- numberdb.org/T456.
Run with SageMath (install the client once with ``sage -pip install numberdb``):
sage -python generate.py # verify every stored entry
sage -python checks.py # independent numerical checks
sage -python generate.py --fill # write entries to the existing draft
Private drafts require NUMBERDB_API_KEY, or NUMBERDB_KEY_FROM_STDIN=1 with
that key piped on stdin. --fill only writes entries and attachments; it never
publishes or offers the draft for review.
The lower-tail inverse is computed with exact rational bisection endpoints.
At positive x, F(x) = 1 - 2*sum((-1)**(j-1)*exp(-2*j*j*x*x), j>=1).
The omitted tail after N terms has absolute value at most
2*exp(-2*(N+1)**2*x*x); this bound is added to the Arb ball, including all
rounding error. A bisection step is taken only when F(mid)-p has a certified
sign. The result encloses both endpoints of the final bracket.
The working precision is bits(digits, losing=128). The bracket width target
is 10**(-digits-12). Checks measure the actual retained accuracy rather than
assuming the guard suffices. The measured minimum was 110 decimal digits.
The conventional 14 probabilities come from the
shared quantile family, not a uniform grid. See checks.py for an independent
MPFI evaluation of the transformed positive theta series and its tail bound.
"""
import os
import sys
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
WORKING_GUARD = 128
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')
def cdf(x, field):
"""Certified alternating series, with the first omitted term bound."""
x = field(x)
if x <= 0:
return field(0)
if not x > 0:
raise ValueError('the argument must be certainly positive')
total = field(1)
scale = -2*x*x
threshold = field(2)**(-field.precision()+24)
for j in range(1, 100000):
term = 2*(scale*(j*j)).exp()
total = total - term if j % 2 else total + term
tail = 2*(scale*((j+1)*(j+1))).exp()
if tail < threshold:
result = total.add_error(tail.upper())
assert result.is_finite()
return result
raise ArithmeticError('series did not converge within the term limit')
def quantile(p, digits=100):
"""Return a ball enclosing the whole certified rational bracket."""
p = QQ(p)
if not 0 < p < 1:
raise ValueError('p must be between zero and one')
field = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
lo, hi = QQ(1)/4, QQ(4)
while not cdf(lo, field) < field(p):
lo /= 2
while not cdf(hi, field) > field(p):
hi *= 2
target = QQ(1)/10**(digits+12)
while hi-lo > target:
mid = (lo+hi)/2
delta = cdf(mid, field)-field(p)
if delta < 0:
lo = mid
elif delta > 0:
hi = mid
else:
# Never interpret an ambiguous ball comparison as a sign.
field = RealBallField(field.precision()+64)
result = field(lo).union(field(hi))
assert result.is_finite() and result.rad() > 0
assert field(result.rad()) < field(target)
return result
class KolmogorovQuantiles(numberdb.Generator):
table = 'T456'
parameters = ('p',)
type = 'R'
digits = 100
rigour = 'proven'
files = ('generate.py', 'checks.py', 'check.py', 'numerical-evidence.json',
'build-notes.md')
def enumerate(self):
for p in PROBABILITIES:
yield {'p': p}
def value(self, params, digits):
return quantile(params['p'], digits)
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
numberdb.configure(api_key=sys.stdin.read().strip())
generator = KolmogorovQuantiles()
if '--fill' in sys.argv:
print(generator.publish(message='add certified Kolmogorov quantiles'))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)