back to table · edit · history · where entries came from · files · download
8478 bytes, as of the version from 2026-10-26 18:42. Recorded here, not run.
"""Independent checks for T523; run: sage -python checks.py.
No generator helper is used to construct the expected moments. Compute the
formal logarithm of the factorial series for ALL graphs, then use Spencer's
connection to excursion moments: Janson (2007), equations (34), (36), (63)-(64).
Coefficient convolution is explicit, so no Sage power-series initialization is
required. Compare with the generator in exact QQ after stripping sqrt(pi/2).
Also use Arb Gamma directly for the displayed moment formula, known Table 1
values, Maxwell/Beta mean, and positive log-convexity margins. This is build
validation, not the separate campaign critique.
"""
import json
import os
import sys
from fractions import Fraction
from decimal import Decimal
from math import factorial
from pathlib import Path
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
import generate
def convolution(a, b, degree):
return [sum((a[j] * b[k-j] for j in range(k+1)), QQ(0))
for k in range(degree+1)]
def formal_log(series):
degree = len(series) - 1
assert series[0] == 1
u = [QQ(0)] + series[1:]
power = [QQ(1)] + [QQ(0)] * degree
out = [QQ(0)] * (degree+1)
for j in range(1, degree+1):
power = convolution(power, u, degree)
for k in range(degree+1):
out[k] += QQ((-1)**(j+1)) / QQ(j) * power[k]
return out
def independent_factors(max_order):
degree = max_order - 1
# Control: log(1+z) = z-z²/2+z³/3-... .
control = formal_log([QQ(1), QQ(1)] + [QQ(0)]*(degree-1))
assert control == [QQ(0)] + [QQ((-1)**(k+1))/QQ(k) for k in range(1, degree+1)]
# Control: log(exp(z)) = z through the full truncation order.
control = formal_log([QQ(1)/QQ(factorial(k)) for k in range(degree+1)])
assert control == [QQ(0), QQ(1)] + [QQ(0)]*(degree-1)
e = [QQ(factorial(6*r)) / QQ(2**(5*r)*3**(2*r)*factorial(3*r)*factorial(2*r))
for r in range(degree+1)]
c = formal_log(e)
assert c[1:4] == [QQ(5)/QQ(24), QQ(5)/QQ(16), QQ(1105)/QQ(1152)]
result = {0: QQ(1), 1: QQ(1)/QQ(2)}
# M_n = n! 2^(2-3n/2) sqrt(pi) c_(n-1) / Gamma(3(n-1)/2).
for n in range(2, max_order+1):
m = n // 2
if n % 2:
result[n] = (QQ(factorial(n)) * QQ(2)**(1-3*m) * c[n-1]
/ QQ(factorial(3*m-1)))
else:
g = QQ(1)
for j in range(3*m-2):
g *= QQ(2*j+1)/QQ(2)
result[n] = QQ(factorial(n)) * QQ(2)**(2-3*m) * c[n-1] / g
return result
def assert_tight(value, bits=450):
assert value.is_finite() and value > 0
assert value.rad_as_ball() < abs(value) * value.parent()(2)**(-bits)
def written_bounds(text):
"""Exact endpoints of the interval a decimal claims, one last place."""
decimal = Decimal(text)
center = Fraction(decimal)
unit = Fraction(10) ** decimal.as_tuple().exponent
def rational(x):
return QQ(x.numerator)/QQ(x.denominator)
return rational(center-unit), rational(center+unit)
def stored_checks():
"""Read the private draft and check its written intervals independently."""
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
key=sys.stdin.read().strip()
if '=' in key and key.split('=',1)[0].strip().isidentifier():
key=key.split('=',1)[1].strip().strip('"\'')
os.environ['NUMBERDB_API_KEY']=key
doc=numberdb.table('T523')
block=doc['Numbers']
if isinstance(block, dict):
rows=[{'params':{'n':str(n)},**v} for n,v in block.items()]
else:
rows=block
assert [r['params']['n'] for r in rows] == [str(n) for n in range(1,13)]
q=independent_factors(12)
field=RealBallField(768)
vals={0:field(1)}
for r in rows:
n=int(r['params']['n']); text=r['number']
if n%2 == 0:
assert '.' not in text and 'e' not in text.lower()
assert QQ(text) == q[n]
vals[n]=field(q[n])
else:
assert len(Decimal(text).as_tuple().digits) == 100
lo,hi=written_bounds(text)
expected=field(q[n])*(field.pi()/field(2)).sqrt()
assert_tight(expected,700)
assert field(lo) < expected < field(hi)
vals[n]=field((lo+hi)/QQ(2), (hi-lo)/QQ(2))
# Recurrence recovered from the stored moments, including exact seed K_0.
kval={0:-field(1)/field(2)}
for n in range(1,13):
kval[n]=(vals[n]*field(QQ(3*n-1)/QQ(2)).gamma()
*field(2)**(QQ(n)/QQ(2))
/ (field(4)*field.pi().sqrt()*factorial(n)))
residual=kval[n]-(field(QQ(3*n-4)/QQ(4))*kval[n-1]
+sum((kval[j]*kval[n-j] for j in range(1,n)),field(0)))
assert 0 in residual
for n in range(1,12):
assert vals[n-1]*vals[n+1]-vals[n]**2 > 0
report=generate.ExcursionMoments().verify(sample=None)
assert report.ok
print(report)
print(json.dumps({'stored_rows_checked':12,'exact_rationals':6,
'100_digit_enclosures_contain_independent_values':6,
'stored_recurrence_residuals_contain_zero':12,
'stored_log_convexity_margins_positive':11}))
def run():
limit = generate.MAX_ORDER
q = independent_factors(limit)
# Repository's comparison helper detects exact disagreement, not ball overlap.
from check import agrees_with, exactness, measure, names_its_rings
actual_q = {n:generate.rational_factor(n) for n in range(limit+1)}
assert not agrees_with(actual_q, lambda n:q[n])
assert not exactness(actual_q)
assert not names_its_rings(str(Path(generate.__file__)))
assert q[0] == 1 and q[1] == QQ(1)/QQ(2)
assert q[2] == QQ(5)/QQ(12)
assert q[3] == QQ(15)/QQ(64)
assert q[4] == QQ(221)/QQ(1008)
f = RealBallField(768)
# Maxwell radial mean times the Beta integral across bridge time.
# Gamma controls catch a bad special-function coercion.
assert 1 in f(1).gamma() and 1 in f(2).gamma()
chi_mean = f(2).sqrt() * f(2).gamma() / f(QQ(3)/QQ(2)).gamma()
beta = f(QQ(3)/QQ(2)).gamma()**2 / f(3).gamma()
mean = chi_mean * beta
assert_tight(mean, 700)
assert mean.overlaps((f.pi()/f(8)).sqrt())
g = generate.ExcursionMoments()
assert list(g.enumerate()) == [{'n':str(n)} for n in range(1,limit+1)]
values = {n:g.value({'n':str(n)},100) for n in range(1,limit+1)}
assert not exactness(values)
accuracy=[]
rows=[]
for n,v in values.items():
# Uncancelled equation (6), evaluated with Arb Gamma, independently
# checks the generator's half-integer simplification over every order.
rhs=(f(4)*f.pi().sqrt()*f(2)**(-QQ(n)/QQ(2))*factorial(n)
/ f(QQ(3*n-1)/QQ(2)).gamma()*f(generate.takacs(n)))
assert_tight(rhs, 700)
assert rhs.overlaps(f(v))
if n % 2:
assert_tight(v)
accuracy.append(int(v.accuracy()))
expected=f(q[n])*(f.pi()/f(2)).sqrt()
assert_tight(expected,700)
assert expected.overlaps(f(v))
else:
assert type(v).__name__ == 'Rational' and v == q[n]
text=to_text(v,digits=100)
if n % 2:
assert len(Decimal(text).as_tuple().digits) == 100
lo, hi = written_bounds(text)
assert f(lo) < expected < f(hi)
rows.append({'params':{'n':str(n)},'number':text})
# Strict log-convexity, including M_0 = 1.
all_values={0:f(1),**{n:f(v) for n,v in values.items()}}
for n in range(1,limit):
assert all_values[n-1]*all_values[n+1]-all_values[n]**2 > 0
print(json.dumps({'independent_exact_factors':{str(n):str(x) for n,x in q.items()},
'controls':['log(1+z)','log(exp(z))','graph c1,c2,c3','Gamma(1)=Gamma(2)=1','M0=1','Table 1 moments n=1..4','Maxwell/Beta mean'],
'independent_graph_comparisons':limit-1,
'gamma_formula_comparisons':limit,
'strict_log_convexity':limit-1,
'minimum_odd_accuracy_bits':min(accuracy),
'working_precision_bits':int(values[1].parent().precision()),
'serialized_measure':measure({r['params']['n']:r['number'] for r in rows}),
'rows':rows},indent=2))
return rows
if __name__ == '__main__':
if '--stored' in sys.argv or os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
stored_checks()
else:
run()