"""Independent numerical checks for numberdb.org/T524. Download this file, generate.py and check.py; run sage -python checks.py. For the stored draft, run sage -python checks.py --stored with an API key. Prints JSON evidence and never writes to the website. Independent construction: let G_s(0,0) be the diagonal resolvent at energy 1/2 of -1/2 d^2/dx^2 + s|x|. On x>=0 its decaying Green function solves f''=(1+2sx)f, f'(0)=-1. Writing f=exp(-x) sum s^n P_n(x) gives P_0=1, P_n''-2P_n'=2x P_(n-1), P_n'(0)=P_n(0) for n>0. We solve these polynomial linear equations coefficient by coefficient over Python Fraction, without using the generator's quadratic recurrence. Feynman-Kac and Brownian scaling (Janson, Appendix C.2) imply P_n(0)=(-1)^n 2^(3n/2) Gamma((3n+1)/2) M_n/(n! sqrt(pi)). A second exact route uses the Airy asymptotic-series quotient, Janson (137), with factorial coefficients (104)-(107), again without the D recurrence. The independent marginal check uses E|N(0,v)|=sqrt(2v/pi) and integral_0^1 sqrt(t(1-t)) dt = B(3/2,3/2)=pi/8. """ 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.real_mpfi import RealIntervalField from fractions import Fraction as F from decimal import Decimal from math import factorial from pathlib import Path import hashlib import json import os import sys import check import generate R = RealBallField(768) I = RealIntervalField(768) def qq(x): return QQ(x.numerator)/QQ(x.denominator) def ode_polynomials(maximum): ps = [[F(1)]] for n in range(1, maximum+1): rhs = [F(0)] + [2*x for x in ps[-1]] q = [F(0)]*(len(rhs)+1) for k in range(len(rhs)-1, -1, -1): q[k] = ((k+1)*q[k+1]-rhs[k])/2 p = [q[0]] + [q[k]/(k+1) for k in range(len(rhs))] assert p[0] == p[1] for k, target in enumerate(rhs): second = (k+2)*(k+1)*p[k+2] if k+2 < len(p) else F(0) assert second-2*(k+1)*p[k+1] == target ps.append(p) return ps def airy_quotient(maximum): a = [F(factorial(6*k),576**k*factorial(3*k)*factorial(2*k)) for k in range(maximum+1)] b = [-F(6*k+1,6*k-1)*a[k] for k in range(maximum+1)] quotient = [] for n in range(maximum+1): quotient.append(a[n]-sum((b[j]*quotient[n-j] for j in range(1,n+1)),F(0))) assert a[0] == b[0] == quotient[0] == 1 assert a[1] == F(5,48) and b[1] == -F(7,48) return quotient def independent_moment(n, polynomial): return (R((-1)**n)*R(qq(polynomial[0]))*R(factorial(n))*R.pi().sqrt() / (R(2).sqrt()**(3*n) * R(QQ(3*n+1)/QQ(2)).gamma())) def bounds(text): """The precise interval a NumberDB spelling claims; no rounded parser.""" if '.' not in text and 'e' not in text.lower(): q = qq(F(text)) return q,q d = Decimal(text) unit = F(10)**d.as_tuple().exponent center = F(d) return qq(center-unit), qq(center+unit) def enclosed(text, expected, exact=None): low,high = bounds(text) assert expected.is_finite() and expected.accuracy() >= 700 iv = I(expected) if low == high: assert exact is not None and low == exact and exact in expected else: assert low < QQ(iv.lower()) and QQ(iv.upper()) < high, text assert len(Decimal(text).as_tuple().digits) == 100 def run(stored=False): # Controls first. Their answers are fixed independently of the generator. controls = ode_polynomials(1) assert controls == [[F(1)],[F(-1,2),F(-1,2),F(-1,2)]] assert R(QQ(1)/QQ(2)).gamma().overlaps(R.pi().sqrt()) assert 1 in independent_moment(0,controls[0]) assert R(1).gamma() == 1 and R(2).gamma() == 1 marginal_mean = (R(2)/R.pi()).sqrt()*R(QQ(3)/QQ(2)).gamma()**2/R(3).gamma() assert marginal_mean.is_finite() and marginal_mean.accuracy()>700 assert marginal_mean.overlaps((R.pi()/R(2)).sqrt()/R(4)) g = generate.AbsoluteBridgeAreaMoments() params = list(g.enumerate()) assert params == [{'n':n} for n in range(1,13)] values = {n:g.value({'n':n},100) for n in range(1,13)} for findings in [check.exactness(values),check.names_its_rings(str(Path(generate.__file__))), check.exact_in_disguise(values)]: assert not findings, findings d = generate.d_coefficients(12) ps = ode_polynomials(12) aq = airy_quotient(12) assert not check.agrees_with(dict(enumerate(d)),lambda n:qq(ps[n][0]/F((-2)**n))) assert not check.agrees_with(dict(enumerate(d)),lambda n:qq(aq[n])) assert [d[n] for n in range(5)] == [QQ(1),QQ(1)/QQ(4),QQ(7)/QQ(32),QQ(21)/QQ(64),QQ(1463)/QQ(2048)] known = [QQ(1),QQ(1)/QQ(4),QQ(7)/QQ(60),QQ(21)/QQ(512),QQ(19)/QQ(720)] assert [generate.rational_factor(n) for n in range(5)] == known assert R(values[1]).overlaps(marginal_mean) expected = {n:independent_moment(n,ps[n]) for n in range(13)} records=[] accuracy=[] for n in range(1,13): value = values[n] if n%2: assert value.is_finite() and value.accuracy()>=500 accuracy.append(int(value.accuracy())) else: assert value in QQ text = to_text(value,digits=100) enclosed(text,expected[n],generate.rational_factor(n) if n%2==0 else None) records.append({'params':{'n':str(n)},'number':text, 'rational_factor':str(generate.rational_factor(n)), 'exact':n%2==0}) # Positive, strictly log-convex raw moments, including M_0=1. for n in range(1,12): assert expected[n-1]*expected[n+1]-expected[n]**2 > 0 # Re-check identities and the independent values on the stored spellings. readback=[] if stored: report=g.verify(sample=None) assert report.ok, str(report) document=numberdb.table(g.table) nums=document['Numbers'] if isinstance(nums,dict): nums=[{'params':{'n':str(k)},**(v if isinstance(v,dict) else {'number':v})} for k,v in nums.items()] assert [int(rec['params']['n']) for rec in nums]==list(range(1,13)) from_stored={0:R(1)} for rec in nums: n=int(rec['params']['n']); text=str(rec['number']) enclosed(text,expected[n],generate.rational_factor(n) if n%2==0 else None) low,high=bounds(text) assert (low==high)==(n%2==0) from_stored[n]=R(I(low,high)) readback.append({'n':n,'enclosure_contains_independent_value':True}) # Recover D from the stored moments, then test the recurrence on intervals. dbd={n:from_stored[n]*R(2).sqrt()**n*R(QQ(3*n+1)/QQ(2)).gamma() /(R.pi().sqrt()*R(factorial(n))) for n in range(13)} for n in range(1,13): residual=dbd[n]-R(QQ(3*n-2)/QQ(4))*dbd[n-1]+sum( (dbd[j]*dbd[n-j] for j in range(1,n)),R(0))/R(2) assert residual.is_finite() and residual.contains_zero() assert R(residual.rad())0 return { 'table':'T524','stage':'stored' if stored else 'pre-fill', 'skill_manifest':'agents/runs/20261026T191343Z-build-skills/manifest.json', 'skill_sha256':'b4d51ded0d52cfe3dcc0859de0845c3879b0c2363273ba7e809a7e06ded40b18', 'environment':g.environment(),'working_bits':numberdb.bits(100,losing=generate.WORKING_GUARD), 'minimum_odd_relative_accuracy_bits':min(accuracy), 'independent_working_bits':768,'ode_and_airy_orders_checked':list(range(13)), 'outside_family_check':'Folded-normal marginal mean and beta integral, independent of the moment recurrence', 'published_table2_positive_moments_checked':[1,2,3,4], 'controls':'P0=1, P1=-(1+x+x^2)/2; M0=1; Gamma(1)=Gamma(2)=1 and Gamma(1/2)=sqrt(pi)', 'exactness':True,'serialized_enclosures_checked':12,'strict_log_convexity_checks':11, 'raw_measurement':check.measure(values), 'serialized_measurement':check.measure({n:r['number'] for n,r in enumerate(records,1)}), 'rows':records,'stored_checks':readback, 'generator_sha256':hashlib.sha256(Path(generate.__file__).read_bytes()).hexdigest(), 'checks_sha256':hashlib.sha256(Path(__file__).read_bytes()).hexdigest() } if __name__=='__main__': if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1': os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip() print(json.dumps(run('--stored' in sys.argv),indent=2))