"""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)