back to table · edit · history · where entries came from · files · download
6810 bytes, as of the version from 2026-10-26 12:40. Recorded here, not run.
"""Independent checks for Scorer Hi -- numberdb.org/T515.
sage -python checks.py
sage -python checks.py --stored # authenticated access while it is a draft
Certified quadrature integrates exp(x*t-t^3/3)/pi over [0,12]. For t>=T,
the exponent has derivative x-t^2 <= x-T^2 <0. Its infinite tail is at most
exp(x*T-T^3/3)/(pi*(T^2-x)). Adding this bound produces a rigorous enclosure
independent of the Taylor coefficient recurrence. Quadrature is controlled
first with integrals of t^2 and exp(-t), whose answers are known exactly or
have elementary rigorous enclosures. mpmath comparisons are numerical only.
"""
import os
import sys
import json
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_mpfi import RealIntervalField
from sage.version import version as sage_version
from mpmath import mp
import mpmath
from generate import ScorerHi, ARGUMENTS, scorer_hi
from check import exactness, measure, names_its_rings
R = RealBallField(665)
C = ComplexBallField(665)
TOL = R(10)**-155
def narrow(v, relative_bits=450):
assert v.is_finite(), 'nonfinite comparison'
assert v != 0
assert v.accuracy() >= relative_bits, ('wide comparison', v.accuracy())
def quadrature(x):
x = QQ(x)
assert x < 144
q = C.integral(lambda t, _: (C(x)*t-t**3/3).exp()/C.pi(),
0, 12, rel_tol=TOL, abs_tol=TOL)
assert q.real().is_finite() and q.imag().is_finite() and q.imag().contains_zero()
tail = (R(x)*12-R(12)**3/3).exp()/(R.pi()*(144-R(x)))
result = q.real().add_error(tail)
narrow(result)
return result, tail
def controls():
a = C.integral(lambda t, _: t**2, 0, 3, rel_tol=TOL, abs_tol=TOL)
b = C.integral(lambda t, _: (-t).exp(), 0, 1,
rel_tol=TOL, abs_tol=TOL)
narrow(a.real()); narrow(b.real())
assert (9 in a.real()) and a.imag().contains_zero()
assert b.real().overlaps(1-R(-1).exp()) and b.imag().contains_zero()
assert not (8 in a.real()), 'quadrature control cannot reject wrong answer'
print('controls: integral t^2=9 and integral exp(-t)=1-exp(-1); wrong 8 rejected', flush=True)
def mp_value(x, dps, derivative=0):
with mp.workdps(dps):
xq = QQ(x)
xm = mp.mpf(int(xq.numerator()))/int(xq.denominator())
v = mp.scorerhi(xm) if derivative == 0 else mp.diff(mp.scorerhi, xm, derivative)
return mp.nstr(v, dps)
def close_mp(v, text, tol):
# This tolerance is a numerical cross-check, never a certified mp bound.
narrow(v)
return abs(v-R(text)) < tol*max(R(1),abs(v))
def main(read_stored=False):
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
numberdb.configure(api_key=sys.stdin.read().strip())
controls()
assert not names_its_rings(os.path.join(os.path.dirname(__file__),'generate.py'))
g=ScorerHi()
vals={p['x']:g.value(p,100) for p in g.enumerate()}
assert not exactness(vals)
texts={x:to_text(v,digits=100) for x,v in vals.items()}
print('versions',sage_version,numberdb.__version__,mpmath.__version__)
print('stored-size',measure(texts))
evidence=[]
quads={}
for x,v in vals.items():
q,tail=quadrature(x);quads[x]=q
assert v.overlaps(q), ('integral disagreement',x)
assert close_mp(v,mp_value(x,150),R(10)**-140), ('mp150',x)
assert close_mp(v,mp_value(x,200),R(10)**-140), ('mp200',x)
# Independently differentiated mpmath implementation checks forcing.
d2=R(mp_value(x,200,2))
ode=d2-R(QQ(x))*v-1/R.pi()
assert abs(ode) < R(10)**-140*max(R(1),abs(v)), ('ODE',x)
assert not abs(d2-R(QQ(x))*v+1/R.pi()) < R(10)**-100, ('wrong sign accepted',x)
with mp.workdps(200):
xq=QQ(x);xm=mp.mpf(int(xq.numerator()))/int(xq.denominator())
gi=R(mp.nstr(mp.scorergi(xm),200))
bi=C(QQ(x)).airy_bi().real();narrow(bi)
assert abs(v+gi-bi) < R(10)**-140*max(R(1),abs(v)), ('Gi+Hi=Bi',x)
row={'x':x,'number':texts[x],'series_accuracy_bits':int(v.accuracy()),
'quadrature_accuracy_bits':int(q.accuracy()),'quadrature_tail_bound':str(tail),
'series_ball':str(v),'quadrature_ball':str(q)}
evidence.append(row)
print(json.dumps(row),flush=True)
# Initial conditions checked using independent Airy values and mp derivative.
assert vals['0'].overlaps(2*C(0).airy_bi().real()/3)
derivative=2/(R(3)**(QQ(5)/6)*R(QQ(1)/3).gamma())
assert close_mp(derivative,mp_value('0',200,1),R(10)**-180)
# Positivity and strict monotonicity also follow by differentiating the
# positive integral. Check both properties on the whole chosen selection.
ordered=sorted(vals,key=QQ)
assert all(vals[x]>0 for x in ordered)
assert all(vals[b]>vals[a] for a,b in zip(ordered,ordered[1:]))
if read_stored:
document=numberdb.table('T515')
entries=document['Numbers']
print('stored-shape',type(entries).__name__)
by_x={str(x):(entry['number'] if isinstance(entry,dict) else entry)
for x,entry in entries.items()}
assert list(by_x)==list(ARGUMENTS)
for x,text in by_x.items():
# A decimal denotes centre +/- one unit in its last place.
from decimal import Decimal
dec=Decimal(text)
centre=QQ(str(dec.as_integer_ratio()[0]))/QQ(str(dec.as_integer_ratio()[1]))
ulp=QQ(10)**int(dec.as_tuple().exponent)
qi=RealIntervalField(665)(quads[x])
assert centre-ulp<=qi.lower() and centre+ulp>=qi.upper(), ('stored',x)
assert abs(R(centre)-R(mp_value(x,200))) < R(ulp), ('stored-mp',x)
# Check the differential equation and sum identity on the stored
# interval, allowing the independently measured numerical mp tolerance.
value=R(centre).add_error(R(ulp))
d2=R(mp_value(x,200,2))
tolerance=R(10)**-180*max(R(1),abs(value))
assert (d2-R(QQ(x))*value-1/R.pi()).add_error(tolerance).contains_zero()
with mp.workdps(200):
xq=QQ(x);xm=mp.mpf(int(xq.numerator()))/int(xq.denominator())
gi=R(mp.nstr(mp.scorergi(xm),200))
assert (value+gi-C(QQ(x)).airy_bi().real()).add_error(tolerance).contains_zero()
print('stored independent containment, ODE and sum',x,flush=True)
report=g.verify(sample=None);print(report);assert report.ok
print('PASS: all nine values, independent quadrature, mp150/mp200, ODE, sum, initial values',flush=True)
if __name__=='__main__':
main('--stored' in sys.argv or os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1')