back to table · edit · history · where entries came from · files · download
10980 bytes, as of the version from 2026-10-26 22:20 (current). Recorded here, not run.
"""Independent checks for T530; run with sage -python checks.py [--stored].
The resolvent method is derived in derivation.md. It shares neither the
K,Q,P,L recurrences nor the parity simplification with generate.py. Exact
Airy-series division provides a third coefficient computation. The repository
helpers in check.py supply exactness, agreement, ring and size checks.
--stored additionally fetches the draft and checks actual written intervals;
set NUMBERDB_API_KEY or NUMBERDB_KEY_FROM_STDIN=1 for a private draft.
"""
import os
import sys
import json
import hashlib
from pathlib import Path
from math import factorial
from fractions import Fraction
from decimal import Decimal
import numberdb.sage as numberdb
from numberdb._write import to_text
from numberdb._generate import _stored_entries
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
import mpmath as mp
from check import exactness, agrees_with, names_its_rings, measure
from one_value import moment
from generate import PositiveBrownianMoments, coefficients, rational_factor, MAX_ORDER
def derivative(poly):
return [i*poly[i] for i in range(1, len(poly))] or [QQ(0)]
def add(a, b):
return [(a[i] if i < len(a) else QQ(0))
+ (b[i] if i < len(b) else QQ(0))
for i in range(max(len(a),len(b)))]
def scale(poly, c):
return [c*x for x in poly]
def resolvent(last):
"""V_n(x)=P_n(x)+exp(-x)Q_n(x), for x>=0 and Exp(1/2) time."""
polynomial = [QQ(1)]
exponential = [QQ(0)]
constants = [QQ(1)]
for n in range(1,last+1):
rhs = [QQ(0)] + scale(polynomial, QQ(2*n))
solution = [QQ(0)]
term = rhs
while any(term):
solution = add(solution, term)
term = derivative(derivative(term))
assert add(solution, scale(derivative(derivative(solution)),QQ(-1))) == rhs
rhs_exp = [QQ(0)] + scale(exponential, QQ(2*n))
derivative_solution = [QQ(0)]
term = scale(rhs_exp, QQ(1)/2)
while any(term):
derivative_solution = add(derivative_solution, term)
term = scale(derivative(term), QQ(1)/2)
# Continuity of V_n and V_n' at zero matches C_n exp(x) at x<0.
constant = (derivative(solution)[0] + derivative_solution[0] - solution[0])/QQ(2)
exp_solution = [constant] + [x/QQ(i+1) for i,x in enumerate(derivative_solution)]
check_exp = add(scale(derivative(exp_solution),QQ(2)),
scale(derivative(derivative(exp_solution)),QQ(-1)))
# Trailing zero coefficients do not change a polynomial.
def normalized(a):
a = a[:]
while len(a)>1 and a[-1]==0: a.pop()
return a
assert normalized(check_exp) == normalized(rhs_exp)
assert derivative(solution)[0]+derivative(exp_solution)[0]-exp_solution[0] == solution[0]+exp_solution[0]
polynomial, exponential = solution, exp_solution
constants.append(polynomial[0]+exponential[0])
if n == 1:
assert normalized(polynomial) == [QQ(0),QQ(2)]
assert normalized(exponential) == [QQ(1)]
assert constants[:3] == [QQ(1),QQ(1),QQ(17)/2]
return constants
def airy_quotient(last):
alpha = [QQ(factorial(6*n))/(QQ(576)**n*factorial(3*n)*factorial(2*n))
for n in range(last+1)]
beta = [QQ(1)]
for n in range(1,last+1):
beta.append(alpha[n]+QQ(3*(2*n-1))/4*beta[n-1])
denominator = [QQ(2)*alpha[n]/QQ(1-6*n) for n in range(last+1)]
result = []
for n in range(last+1):
result.append((alpha[n]+beta[n]-sum(denominator[j]*result[n-j]
for j in range(1,n+1)))/denominator[0])
assert result[:3] == [QQ(1),QQ(1)/2,QQ(17)/16]
return result
def exact_decimal(text):
f = Fraction(text)
return QQ(f.numerator)/QQ(f.denominator)
def written_bounds(text):
centre = exact_decimal(text)
if '.' not in text and 'e' not in text.lower():
return centre, centre
unit = QQ(10)**int(Decimal(text).as_tuple().exponent)
return centre-unit, centre+unit
def check_written(text, expected, exact=None):
lo,hi = written_bounds(text)
if exact is not None:
assert lo==hi==exact
else:
assert lo<hi and expected.is_finite()
assert expected.accuracy() > 700
assert expected > lo and expected < hi
def gaussian_checks():
mp.mp.dps = 85
control = mp.quad(lambda t: mp.sin(t)**2*mp.cos(t)**2, [0,mp.pi/2])
assert abs(control-mp.pi/16) < mp.mpf('1e-80')
mean = mp.quad(lambda t: mp.sqrt(t/(2*mp.pi)),[0,1])
second = 2/(3*mp.pi)*mp.quad(
lambda t: mp.sin(t)**2*(mp.cos(t)+mp.sin(t)*(mp.pi/2+t))*mp.cos(t),
[0,mp.pi/2])
assert abs(mean-mp.sqrt(2/mp.pi)/3)<mp.mpf('1e-80')
assert abs(second-mp.mpf(17)/96)<mp.mpf('1e-80')
return {'working_decimal_digits':85,'nonzero_control':'integral sin^2(theta) cos^2(theta) = pi/16',
'control_error':mp.nstr(abs(control-mp.pi/16),6),
'mean':mp.nstr(mean,80),'second_moment':mp.nstr(second,80),
'method':'Gaussian marginal mean and joint positive-part product integral; agreement check, not certification'}
def run(stored=False):
base = Path(__file__).resolve().parent
gen = PositiveBrownianMoments()
params = list(gen.enumerate())
assert params == [{'n':n} for n in range(1,13)]
assert not names_its_rings(str(base/'generate.py'))
values = {p['n']:gen.value(p,100) for p in params}
assert not exactness(values)
gaussian = gaussian_checks() # nonzero control precedes this measurement
fk = resolvent(MAX_ORDER)
ell = coefficients(MAX_ORDER)[3]
assert not agrees_with(dict(enumerate(ell)),lambda n:airy_quotient(MAX_ORDER)[n])
assert not agrees_with(dict(enumerate(ell)),lambda n:fk[n]/(QQ(2)**n*factorial(n)))
assert rational_factor(0)==1
for n,c in {1:QQ(1)/3,2:QQ(17)/96,3:QQ(251)/1260,4:QQ(6989)/40960}.items():
assert rational_factor(n)==c
high = RealBallField(768)
independent = {n:high(fk[n])/(high(2)**(QQ(3*n)/2)*(high(QQ(3*n)/2+1)).gamma())
for n in range(MAX_ORDER+1)}
# Perman and Wellner, printed Table 1, positive-motion column. Odd
# coefficients multiply 1/sqrt(2*pi), not sqrt(2/pi).
printed = {1:QQ(2)/3, 2:QQ(17)/96, 3:QQ(251)/630,
4:QQ(6989)/40960, 5:QQ(188267)/360360,
6:QQ(37235311)/132120576, 7:QQ(1451280043)/1396755360,
8:QQ(522258818027)/797253304320, 9:QQ(8096107769)/2930256000,
10:QQ(1011667945773427)/515911471595520}
for n,c in printed.items():
assert rational_factor(n)==(c/2 if n%2 else c)
rows = []
accuracy = []
for n,v in values.items():
exact = fk[n]/(QQ(2)**(3*n//2)*factorial(3*n//2)) if n%2==0 else None
if n%2:
assert v.is_finite() and v.accuracy()>=500
accuracy.append(int(v.accuracy()))
assert v.overlaps(independent[n])
else:
assert v==exact and (v in independent[n])
program_value = moment(n)
if n%2:
assert program_value.is_finite() and program_value.accuracy()>490
assert program_value.overlaps(independent[n])
else:
assert program_value==exact
text = to_text(v,digits=100)
check_written(text,independent[n],exact)
rows.append({'params':{'n':str(n)},'number':text})
for n in range(1,MAX_ORDER):
assert independent[n-1]*independent[n+1]-independent[n]**2 > 0
all_recurrences = True
k,q,p,l = coefficients(MAX_ORDER)
for n in range(1,MAX_ORDER+1):
assert k[n]==QQ(3*n-4)/4*k[n-1]+sum(k[j]*k[n-j] for j in range(1,n))
assert q[n]==QQ(3*n-2)/2*q[n-1]+2*sum(k[j]*q[n-j] for j in range(1,n+1))
assert p[n]==sum(k[j]*p[n-j] for j in range(1,n+1))
assert 2*l[n]==p[n]+sum(p[j]*q[n-j] for j in range(n+1))
stored_report = None
if stored:
doc = numberdb.table('T530')
records = _stored_entries(doc)
assert len(records)==12
assert list(records)==[str(n) for n in range(1,13)]
balls = {0: high(1)}
for identity,text in records.items():
n=int(identity)
exact = fk[n]/(QQ(2)**(3*n//2)*factorial(3*n//2)) if n%2==0 else None
check_written(text,independent[n],exact)
lo,hi=written_bounds(text)
balls[n]=high((lo+hi)/2).add_error(high((hi-lo)/2))
# Invert the moment formula on the stored interval, then test
# the independently computed coefficient and convolution.
recovered = balls[n]*(high(QQ(3*n)/2+1)).gamma()*high(2)**(QQ(n)/2)/factorial(n)
assert (ell[n] in recovered)
assert (airy_quotient(MAX_ORDER)[n] in recovered)
assert recovered.accuracy()>320
for n in range(1,MAX_ORDER):
assert balls[n-1]*balls[n+1]-balls[n]**2 > 0
stored_report={'entries':12,'all_independent_enclosures_contained':True,
'inverted_moment_formula_and_convolution':True,'strict_log_convexity':True}
measured=measure({str(n):to_text(v,100) for n,v in values.items()})
evidence={'table':'T530','orders':[1,12],'entries':12,'exact_even_moments':6,
'certified_odd_moments':6,'stored_digits':100,'working_bits':numberdb.bits(100,losing=180),
'minimum_relative_accuracy_bits':min(accuracy),'per_odd_order_accuracy_bits':accuracy,
'independent_working_bits':768,'resolvent_coefficients':[str(x) for x in fk],
'rational_factors':[str(rational_factor(n)) for n in range(1,13)],
'airy_quotient_all_orders':True,'feynman_kac_all_orders':True,
'all_stated_recurrences':all_recurrences,'janson_table8_orders_1_to_4':True,
'perman_wellner_table1_orders_1_to_10':True,'program_checked_all_orders':True,
'written_interval_containment':True,'strict_log_convexity_orders_1_to_11':True,
'gaussian_checks':gaussian,'measurement':measured,'stored_readback':stored_report,
'instruction_manifest':'agents/runs/20261026T220819Z-build-skills/manifest.json',
'instruction_sha256':'b4d51ded0d52cfe3dcc0859de0845c3879b0c2363273ba7e809a7e06ded40b18',
'generator_sha256':hashlib.sha256((base/'generate.py').read_bytes()).hexdigest()}
from sage.version import version as sage_version
import importlib.metadata
evidence['versions']={'sage':sage_version,'mpmath':mp.__version__,
'numberdb':importlib.metadata.version('numberdb')}
print('EVIDENCE_JSON='+json.dumps(evidence,sort_keys=True))
print('ROWS_JSON='+json.dumps(rows))
return evidence
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
run(stored='--stored' in sys.argv)