back to table · edit · history · where entries came from · files · download
8586 bytes, as of the version from 2026-10-25 15:31. Recorded here, not run.
"""Independent checks for Struve H -- numberdb.org/T492.
Download beside generate.py and check.py, then run:
sage -python checks.py
sage -python checks.py --stored
The stored mode needs NUMBERDB_API_KEY while the table is a draft.
No writes are made. Independent certified quadrature uses DLMF 11.5.1,
with t=cos(theta), so its integrand is entire for the five held orders.
Controls run first: integrals of theta and sin(theta), H_{1/2}(1),
and the integral of J_0 at 0 and 1. Every comparison checks finite bounds
and a strict radius threshold before allowing overlap or containment.
"""
import json
import os
import sys
from decimal import Decimal
from fractions import Fraction
from functools import lru_cache
from pathlib import Path
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_arb import RealBallField
from sage.rings.rational_field import QQ
import mpmath as mp
from generate import StruveH, ORDERS, arguments, struve_h
from check import exactness, measure, names_its_rings
C = ComplexBallField(480)
R = RealBallField(480)
HALF = QQ(1) / 2
def narrow(v, digits=110):
assert v.is_finite(), 'non-finite enclosure'
assert R(v.rad()) < R(10)**(-digits) * max(R(1), abs(v)), 'wide enclosure'
return v
def real(v):
assert v.real().is_finite() and v.imag().is_finite()
assert v.imag().contains_zero()
return narrow(v.real())
def decimal_parts(text):
d = Decimal(text)
f = Fraction(d)
mid = QQ(f.numerator) / f.denominator
ulp = QQ(10) ** int(d.as_tuple().exponent)
return mid, ulp
def written_ball(text):
mid, ulp = decimal_parts(text)
return narrow(R(mid).add_error(R(ulp)), 97)
def contains_written(text, value):
narrow(value)
mid, ulp = decimal_parts(text)
assert value > R(mid - ulp) and value < R(mid + ulp), ('outside decimal', text)
@lru_cache(None)
def integral_h(nu, x):
assert nu in ORDERS and x >= 0
n, z = C(nu), C(x)
power = int(2 * nu)
v = C.integral(lambda t, analytic: (z*t.cos()).sin()*t.sin()**power,
C(0), C.pi()/2)
return real(2 * (z/2)**n * v / (C.pi().sqrt()*(n+HALF).gamma()))
def half_elementary(nu, x):
z=R(x)
factor=(2/(R.pi()*z)).sqrt()
if nu==HALF:
return narrow(factor*(1-z.cos()))
assert nu==3*HALF
return narrow(factor*(z/2+(1-z.cos())/z-z.sin()))
def bessel_integral(x):
# Integral of J_0's power series: exact rational sum through k=249.
# Consecutive absolute terms decrease after k=249; ratios are bounded
# by x^2/(4*251^2) for all terms after the first omitted term.
x=QQ(x)
term=x
total=QQ(0)
for k in range(250):
total += term
term *= -x*x*QQ(2*k+1)/(4*(k+1)**2*(2*k+3))
ratio=x*x/(4*QQ(251)**2)
assert 0 <= ratio < 1
tail=abs(term)/(1-ratio)
return narrow(R(total).add_error(R(tail)))
def controls():
v=real(C.integral(lambda t,a:t,0,1))
assert (HALF in v) and not v.contains_zero()
v=real(C.integral(lambda t,a:t.sin(),0,C.pi()/2))
assert (1 in v) and not v.contains_zero()
v=integral_h(HALF,QQ(1))
assert v.overlaps(half_elementary(HALF,QQ(1))) and not v.contains_zero()
assert bessel_integral(0).contains_zero()
j=real(C.integral(lambda t,a:t.bessel_J(0),0,1))
assert j.overlaps(bessel_integral(1)) and not j.contains_zero()
for nu in ORDERS:
assert struve_h(nu,QQ(0)).contains_zero()
print('Controls passed: nonzero known integrals, half order, origin, J0 integral.',flush=True)
def run(stored=False):
import platform
from sage.env import SAGE_VERSION
assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
controls()
mp.mp.dps=160
xs=arguments()
gen=StruveH()
rows={}
if stored:
doc=numberdb.table(gen.table)
assert set(doc['Parameters'])=={'nu','x'}
rows={(QQ(n),QQ(x)): val['number'] if isinstance(val,dict) else val
for n, group in doc['Numbers'].items() for x,val in group.items()}
assert set(rows)=={(n,x) for n in ORDERS for x in xs}
vals={}
texts={}
evidence=[]
max_rad=R(0)
max_rel_rad=R(0)
max_quad_rad=R(0)
max_mp_relative=mp.mpf(0)
for i,params in enumerate(gen.enumerate(),1):
n,x=QQ(params['nu']),QQ(params['x'])
v=narrow(gen.value(params,100))
text=rows[(n,x)] if stored else to_text(v,digits=100)
contains_written(text,v)
q=integral_h(n,x)
assert v.overlaps(q), ('integral disagreement',params)
contains_written(text,q)
if n in (HALF,3*HALF):
e=half_elementary(n,x)
assert v.overlaps(e)
contains_written(text,e)
# A different implementation at 160 decimal digits.
m=mp.struveh(mp.mpf(str(n.numerator()))/int(n.denominator()),
mp.mpf(str(x.numerator()))/int(x.denominator()))
mid,ulp=decimal_parts(text)
mm=mp.mpf(str(mid.numerator()))/int(mid.denominator())
uu=mp.mpf(str(ulp.numerator()))/int(ulp.denominator())
assert abs(m-mm)<uu
# Compare mpmath with a 130-digit midpoint of the integral enclosure.
qmid=to_text(q,digits=130)
rel=abs(mp.mpf(qmid)-m)/max(mp.mpf(1),abs(m))
assert rel<mp.mpf('1e-128')
max_mp_relative=max(max_mp_relative,rel)
max_rad=max(max_rad,R(v.rad()))
max_rel_rad=max(max_rel_rad,R(v.rad())/abs(R(v)))
max_quad_rad=max(max_quad_rad,R(q.rad()))
vals[(n,x)]=written_ball(text) if stored else R(v)
texts[(str(n),str(x))]=text
evidence.append({'nu':str(n),'x':str(x),'number':text,
'generator_radius':str(v.rad()),
'integral_radius':str(q.rad()),
'mpmath_relative_error':mp.nstr(rel,12)})
if i%100==0:
print('Checked %d/600 %s rows.'%(i,'stored' if stored else 'generated'),flush=True)
assert not exactness(vals)
# All five orders at all arguments, using stored values whenever present.
for n in ORDERS:
for x in xs:
prev=vals.get((n-1,x))
if prev is None: prev=R(narrow(struve_h(n-1,x,140)))
nxt=vals.get((n+1,x))
if nxt is None: nxt=R(narrow(struve_h(n+1,x,140)))
rhs=2*R(n)/R(x)*vals[n,x]+(R(x)/2)**R(n)/(R.pi().sqrt()*R(n+3*HALF).gamma())
assert narrow(prev+nxt,96).overlaps(narrow(rhs,96)),('recurrence',n,x)
# H_-1 and the corrected J_0 integral identity at every held argument.
bessel_doc=numberdb.table('T187') if stored else None
stored_j_pairs=0
for x in xs:
hm=R(narrow(struve_h(QQ(-1),x,140)))
assert narrow(2/R.pi()-vals[QQ(1),x],96).overlaps(hm)
if bessel_doc and str(x) in bessel_doc['Numbers']['0']:
j0=written_ball(bessel_doc['Numbers']['0'][str(x)])
j1=written_ball(bessel_doc['Numbers']['1'][str(x)])
stored_j_pairs+=1
else:
j0,j1=real(C(x).bessel_J(0)),real(C(x).bessel_J(1))
bracket=j0*(2/R.pi()-vals[QQ(1),x])+j1*vals[QQ(0),x]
rhs=R.pi()/2*R(x)*bracket
exact=bessel_integral(x)
assert narrow(rhs,96).overlaps(exact),('Bessel integral',x)
if x==1:
assert not narrow(R(x)*bracket,96).overlaps(exact)
summary={'table':'T492','mode':'stored' if stored else 'generated',
'entries':len(vals),'orders':[str(n) for n in ORDERS],
'arguments':len(xs),'precision_bits':numberdb.bits(100,losing=64),
'integral_precision_bits':480,'max_generator_radius':str(max_rad),
'max_generator_relative_radius':str(max_rel_rad),
'max_integral_radius':str(max_quad_rad),
'max_mpmath_relative_error':mp.nstr(max_mp_relative,16),
'measure':measure(texts),'recurrence_cases':600,
'half_integer_cases':240,'negative_order_cases':120,
'bessel_integral_cases':120,'stored_bessel_pairs':stored_j_pairs,
'missing_pi_over_2_rejected':True,
'versions':{'sage':SAGE_VERSION,'python':platform.python_version(),
'mpmath':mp.__version__}}
print('SUMMARY '+json.dumps(summary),flush=True)
print('EVIDENCE '+json.dumps(evidence),flush=True)
return summary,evidence
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN'):
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
run(stored='--stored' in sys.argv)