back to table · edit · history · where entries came from · files · download
8712 bytes, as of the version from 2026-10-26 20:49 (current). Recorded here, not run.
"""Independent computations and stored-value checks for numberdb.org/T527.
Run alongside generate.py and the unmodified campaign check.py:
sage -python checks.py
sage -python checks.py --stored # NUMBERDB_API_KEY for the private draft
The polynomial resolvent computation and factorial-series quotient share no
coefficient code with generate.py. See derivation.md for the ODE derivation.
The script prints JSON evidence; redirect stdout to preserve it. No writes.
"""
import os
import sys
import json
from decimal import Decimal
from fractions import Fraction
from math import factorial, comb
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.version import version as sage_version
import generate
from check import exactness, measure, agrees_with, names_its_rings
MANIFEST = 'agents/runs/20261026T204019Z-build-skills/manifest.json'
SKILL_SHA = 'b4d51ded0d52cfe3dcc0859de0845c3879b0c2363273ba7e809a7e06ded40b18'
def derivative(a):
return [j*a[j] for j in range(1,len(a))] or [QQ(0)]
def add(a,b):
out=[QQ(0)]*max(len(a),len(b))
for j,x in enumerate(a):out[j]+=x
for j,x in enumerate(b):out[j]+=x
while len(out)>1 and out[-1]==0:out.pop()
return out
def scaled(a,c):
return [x*c for x in a]
def resolvent_q(count):
# v_n = P_n + exp(-x)R_n; (d²-1)v_n = x v_(n-1), v_n(0)=0.
p,r=[QQ(1)],[-QQ(1)]
out=[QQ(1)]
for n in range(1,count+1):
a=[QQ(0)]+p
pnew=[QQ(0)]
while any(a):
pnew=add(pnew,scaled(a,-1))
a=derivative(derivative(a))
a=[QQ(0)]+r
sprime=[QQ(0)]
power=QQ(2)
while any(a):
sprime=add(sprime,scaled(a,-QQ(1)/power))
a=derivative(a);power*=2
rnew=[-pnew[0]]+[sprime[j]/QQ(j+1) for j in range(len(sprime))]
# Verify the polynomial identities, independently of the recurrence extraction.
assert add(derivative(derivative(pnew)),scaled(pnew,-1))==add([QQ(0)]+p,[QQ(0)])
assert add(derivative(derivative(rnew)),scaled(derivative(rnew),-2))==add([QQ(0)]+r,[QQ(0)])
assert pnew[0]+rnew[0]==0
slope=derivative(pnew)[0]+derivative(rnew)[0]-rnew[0]
out.append((-1)**n*slope)
p,r=pnew,rnew
return out
def quotient_q(count):
alpha=[QQ(factorial(6*n))/QQ(576**n*factorial(3*n)*factorial(2*n))
for n in range(count+1)]
beta=[QQ(1)]
q=[QQ(1)]
for n in range(1,count+1):
beta.append(alpha[n]+QQ(6*n-3)/QQ(4)*beta[-1])
q.append(beta[n]-sum((alpha[j]*q[n-j] for j in range(1,n+1)),QQ(0)))
return q
def direct_moment(n,q,field):
w=sum((q[j]*q[n-j] for j in range(n+1)),QQ(0))
return field(2)**(-field(n)/2)*field(factorial(n))*field(w)/field(QQ(3*n+2)/QQ(2)).gamma()
def written_interval(text):
if '.' not in text and 'e' not in text.lower():
v=QQ(text);return v,v
d=Decimal(text);f=Fraction(d);center=QQ(f.numerator)/QQ(f.denominator)
unit=QQ(10)**int(d.as_tuple().exponent)
return center-unit,center+unit
def check_inside(text,ball):
assert ball.is_finite() and ball.accuracy()>=700
lo,hi=written_interval(text)
field=ball.parent()
if lo==hi:
assert lo in ball
else:
assert ball>field(lo) and ball<field(hi),text
def run(stored=False):
count=generate.MAX_ORDER
high=RealBallField(768)
q=resolvent_q(count)
# Controls, checked before comparisons: nonzero values, both boundary terms.
assert q[:3]==[QQ(1),QQ(3)/4,QQ(59)/32]
assert quotient_q(2)==q[:3]
assert 1 in direct_moment(0,q,high)
assert high(QQ(1)/2).gamma().overlaps(high.pi().sqrt())
# Reflection principle: -min B has the folded-standard-normal mean.
mean=high(2/high.pi()).sqrt()
assert direct_moment(1,q,high).overlaps(mean)
assert mean.accuracy()>=750
# Published values: Majumdar–Comtet (74); Janson Table 5 adds order 4.
controls={0:QQ(1),1:QQ(1),2:QQ(17)/24,3:QQ(123)/140,4:QQ(2963)/3840}
assert all(generate.rational_factor(n)==x for n,x in controls.items())
k,gq=generate.coefficients(count)
assert not agrees_with(dict(enumerate(gq)),lambda n:q[n])
assert not agrees_with(dict(enumerate(gq)),lambda n:quotient_q(count)[n])
assert not exactness(dict(enumerate(k)))
assert not exactness(dict(enumerate(gq)))
assert not names_its_rings(generate.__file__)
gen=generate.DoubleMeanderMoments()
values={int(p['n']):gen.value(p,100) for p in gen.enumerate()}
assert not exactness(values)
written={n:to_text(v.get('number') if isinstance(v,dict) else v,100) for n,v in values.items()}
fresh={n:direct_moment(n,q,high) for n in range(count+1)}
for n in range(1,count):
assert fresh[n-1]*fresh[n+1]-fresh[n]**2>0
# Probabilistic arcsine splitting, with independently computed meander moments.
me=[high.pi().sqrt()*high(2)**(-high(n)/2)*factorial(n)*q[n]
/high(QQ(3*n+1)/2).gamma() for n in range(count+1)]
for n in range(count+1):
conv=high(0)
for j in range(n+1):
a,b=high(QQ(3*j+1)/2),high(QQ(3*(n-j)+1)/2)
conv+=comb(n,j)*me[j]*me[n-j]*a.gamma()*b.gamma()/(a+b).gamma()/high.pi()
assert conv.is_finite() and conv.accuracy()>700
assert conv.overlaps(fresh[n])
rows=[]
for n,text in written.items():
check_inside(text,fresh[n])
v=generate.moment(n)
rows.append({'n':n,'exact_rational_factor':str(generate.rational_factor(n)),
'factor_multiplies_sqrt_2_over_pi':bool(n%2),'written':text,
'relative_accuracy_bits':int(v.accuracy()) if n%2 else None,
'independent_relative_accuracy_bits':int(fresh[n].accuracy()),
'written_interval_contains_independent_enclosure':True})
result={'table':'T527','skill_manifest':MANIFEST,'skill_sha256':SKILL_SHA,
'sage_version':sage_version,'working_bits':numberdb.bits(100,losing=generate.WORKING_GUARD),
'independent_bits':768,'minimum_odd_relative_accuracy_bits':min(x['relative_accuracy_bits'] for x in rows if x['n']%2),
'checks':{'resolvent_polynomial_identities':'all orders 1..12',
'coefficient_agreement_resolvent':'all orders 0..12',
'coefficient_agreement_factorial_series_quotient':'all orders 0..12',
'gamma_normalization_controls':True,'reflection_principle_mean':True,
'published_low_moments':'orders 0..4','arcsine_meander_convolution':'orders 0..12',
'strict_log_convexity':'orders 1..11, with M0=1',
'exactness_and_named_imports':True},
'serialized_measurement':measure(written),'rows':rows}
if stored:
report=gen.verify(sample=None)
assert report.ok,str(report)
doc=numberdb.table('T527')
entries=doc['Numbers']
if isinstance(entries,list):
texts={int(x['params']['n']):str(x['number']) for x in entries}
else:
texts={int(n):str(v['number'] if isinstance(v,dict) else v) for n,v in entries.items()}
assert list(texts)==list(range(1,count+1))
intervals={0:high(1)}
for n,text in texts.items():
assert text==written[n]
check_inside(text,fresh[n])
lo,hi=written_interval(text)
intervals[n]=high((lo+hi)/2, (hi-lo)/2)
# Reconstruct W_n from the database value, testing every displayed formula.
wread=intervals[n]*high(2)**(high(n)/2)*high(QQ(3*n+2)/2).gamma()/factorial(n)
wexact=sum((q[j]*q[n-j] for j in range(n+1)),QQ(0))
assert wread.is_finite() and wread.accuracy()>300
assert 0 in (wread-high(wexact))
for n in range(1,count):
assert intervals[n-1]*intervals[n+1]-intervals[n]**2>0
result['stored_checks']={'generator_verify_all':True,'count':len(texts),
'independent_enclosures_contained':True,
'reconstructed_W_matches_independent_convolution':True,
'strict_log_convexity':True,'enumeration_matches':True}
return result
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
key=sys.stdin.read().strip()
if key.startswith('NUMBERDB_API_KEY='):key=key.split('=',1)[1].strip().strip('\"\'')
os.environ['NUMBERDB_API_KEY']=key
use_stored='--stored' in sys.argv or os.environ.get('NUMBERDB_RESTATING')=='1'
print(json.dumps(run(use_stored),indent=2))