back to table · edit · history · where entries came from · files · download
8732 bytes, as of the version from 2026-10-26 19:49 (current). Recorded here, not run.
"""Independent checks for T525; run sage -python checks.py [--stored].
Offline by default. --stored reads the private draft using NUMBERDB_API_KEY.
Reports JSON to stdout. Feynman-Kac polynomial construction and its controls
are explained in derivation.md. The generator is used only for comparison.
"""
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from math import factorial
from fractions import Fraction
from decimal import Decimal
import importlib.metadata
import json
import os
import sys
import hashlib
from pathlib import Path
from numberdb._write import to_text
from numberdb._compare import digits_of
from numberdb._generate import _stored_entries
from check import exactness, measure, agrees_with, names_its_rings
from generate import AbsoluteBrownianAreaMoments, coefficients, rational_factor, moment
LIMIT = 12
R = RealBallField(768)
def resolvent(limit):
ring = PolynomialRing(QQ, 'y')
y = ring.gen()
p, q = ring(1), ring(0)
result = [QQ(1)]
for n in range(1, limit + 1):
f, h = y*p, y*q
pnext, term = ring(0), f
while term:
pnext += term
term = term.derivative().derivative()
v, term = ring(0), h/QQ(2)
while term:
v += term
term = term.derivative()/QQ(2)
# Integrate coefficient by coefficient; no lazy symbolic machinery.
qnext = ring([QQ(0)] + [v[k]/QQ(k+1) for k in range(v.degree()+1)]) if v else ring(0)
qnext += pnext.derivative()(0) + v(0)
assert pnext - pnext.derivative().derivative() == f
assert 2*qnext.derivative() - qnext.derivative().derivative() == h
assert pnext.derivative()(0) + qnext.derivative()(0) - qnext(0) == 0
if n == 1:
assert pnext == y and qnext == 1
p, q = pnext, qnext
result.append(p(0)+q(0))
return result
def independent_moment(n, origins):
return R(factorial(n))*R(origins[n]) / (R(2)**(R(n)/R(2)) * (R(3*n)/R(2)+1).gamma())
def decimal_claim(text):
value = Decimal(text)
numerator, denominator = Fraction(value).as_integer_ratio()
center = QQ(numerator)/QQ(denominator)
power = int(value.as_tuple().exponent)
ulp = QQ(10)**power
return center, ulp
def written_contains(text, expected, exact_expected=None):
if '/' in text or ('.' not in text and 'e' not in text.lower()):
if exact_expected is not None:
return QQ(text) == exact_expected
return QQ(text) in expected
center, ulp = decimal_claim(text)
return bool(expected > R(center-ulp) and expected < R(center+ulp))
def written_ball(text):
if '/' in text or ('.' not in text and 'e' not in text.lower()):
return R(QQ(text))
center, ulp = decimal_claim(text)
return R(center).add_error(R(ulp))
def main(stored=False):
# Positive and nonzero controls precede every subsequent comparison.
origins = resolvent(LIMIT)
assert origins[:3] == [QQ(1), QQ(1), QQ(9)/QQ(4)]
assert moment(0) == 1
independent = {n:independent_moment(n, origins) for n in range(LIMIT+1)}
assert 1 in independent[0]
mean = R(QQ(2)/QQ(3))*(R(2)/R.pi()).sqrt()
assert independent[1].overlaps(mean)
assert QQ(3)/QQ(8) in independent[2]
assert independent[3].overlaps(R(QQ(263)/QQ(630))*(R(2)/R.pi()).sqrt())
assert QQ(903)/QQ(2560) in independent[4]
# Independent analytic Gaussian-product integral control (derivation.md).
second = R(8)/(R(3)*R.pi())*(R.pi()/16 + R(5)*R.pi()/64)
assert QQ(3)/QQ(8) in second
assert all(x.is_finite() and x > 0 and x.accuracy() > 740
for n,x in independent.items() if n > 0)
d, ell = coefficients(LIMIT)
assert not agrees_with(dict(enumerate(ell)), lambda n:origins[n])
assert not exactness(d + ell)
assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
gen = AbsoluteBrownianAreaMoments()
params = list(gen.enumerate())
assert params == [{'n':n} for n in range(1,LIMIT+1)]
values = {p['n']:gen.value(p,100) for p in params}
assert not exactness(values)
texts = {n:to_text(value,100) for n,value in values.items()}
exact = {2:QQ(3)/8,4:QQ(903)/2560}
# For all even orders, rational simplification is checked separately
# against direct Arb Gamma, while exact L matches the resolvent.
for n in range(2,LIMIT+1,2):
exact[n] = QQ(factorial(n))*origins[n]/QQ(2**(n//2)*factorial(3*n//2))
assert values[n] == exact[n]
for n in range(1,LIMIT+1):
assert written_contains(texts[n],independent[n],exact.get(n)), n
if n%2:
assert digits_of(texts[n]) == 100
assert values[n].accuracy() > 490
assert values[n].overlaps(independent[n])
# Negative controls: the enclosure checker must reject perturbed answers.
center, ulp = decimal_claim(texts[1])
assert not bool(independent[1] > R(center+9*ulp) and independent[1] < R(center+11*ulp))
assert not written_contains('4/8',independent[2],exact[2])
assert not agrees_with({0:QQ(1)},lambda n:QQ(1))
assert agrees_with({0:QQ(1)},lambda n:QQ(2))
checked = texts
if stored:
doc = numberdb.table('T525')
assert doc['Title'] == 'Moments of the absolute Brownian motion area'
entries = _stored_entries(doc)
assert [int(n) for n in entries] == list(range(1,LIMIT+1))
checked = {int(n):text for n,text in entries.items()}
for n,text in checked.items():
assert written_contains(text,independent[n],exact.get(n)), (n,text)
if n%2: assert digits_of(text)==100
else: assert '/' in text or text.isdigit()
intervals = {0:R(1), **{n:written_ball(text) for n,text in checked.items()}}
reconstructed = {n:intervals[n]*(R(2)**(R(n)/2))
*(R(3*n)/2+1).gamma()/R(factorial(n)) for n in range(LIMIT+1)}
for n in range(1,LIMIT+1):
residual = reconstructed[n] - R(QQ(3*n-1)/2)*reconstructed[n-1]
residual += sum((R(d[n-j])*reconstructed[j] for j in range(1,n)),R(0))
assert residual.is_finite() and residual.contains_zero(), n
assert residual.rad_as_ball() < R(QQ(1)/10**80), n
for n in range(1,LIMIT):
gap = intervals[n-1]*intervals[n+1]-intervals[n]**2
assert gap > 0, n
# At T=4, Brownian scaling is multiplication of area by 8.
for n in range(1,LIMIT+1):
scaled_direct = independent[n]*R(4)**(R(3*n)/2)
assert scaled_direct.overlaps(independent[n]*R(8)**n)
report = {
'table':'T525','phase':'stored' if stored else 'offline',
'instruction_manifest':'agents/runs/20261026T193829Z-build-skills/manifest.json',
'skill_sha256':'b4d51ded0d52cfe3dcc0859de0845c3879b0c2363273ba7e809a7e06ded40b18',
'orders':list(range(1,LIMIT+1)),
'source':'Janson 2007, Section 21, Table 3, and Appendix C.1',
'all_resolvent_coefficients_match':True,
'polynomial_ode_and_boundary_residuals_exact':True,
'source_first_four_match':True,
'independent_gaussian_mean_and_second_moment':True,
'negative_controls_reject':True,
'written_intervals_contain_independent_enclosures':True,
'recurrence_residuals_contain_zero':True,
'strict_log_convexity_orders':list(range(1,LIMIT)),
'working_bits':numberdb.bits(100,losing=180),
'min_odd_relative_accuracy_bits':min(int(v.accuracy()) for n,v in values.items() if n%2),
'independent_working_bits':768,
'min_independent_relative_accuracy_bits':min(int(v.accuracy()) for n,v in independent.items() if n>0),
'measurement':measure(checked),
'values':checked,
'rational_factors':{n:str(rational_factor(n,ell[n])) for n in range(1,LIMIT+1)},
'source_sha256':{name:hashlib.sha256(Path(__file__).with_name(name).read_bytes()).hexdigest()
for name in ('generate.py','checks.py','check.py','argument-selection.md','derivation.md','skill-manifest.json')},
'versions':{'python':sys.version.split()[0],'numberdb':numberdb.__version__},
}
try:
from sage.env import SAGE_VERSION
report['versions']['sage'] = SAGE_VERSION
except ImportError: pass
for package in ['python-flint','passagemath-standard','passagemath-symbolics']:
try: report['versions'][package] = importlib.metadata.version(package)
except importlib.metadata.PackageNotFoundError: pass
print(json.dumps(report,indent=2))
return report
if __name__ == '__main__':
main('--stored' in sys.argv)