"""Check the Gaussian-limit correction in T606 without changing the table.

Run with SageMath, numberdb installed, and the unchanged T606 generate.py nearby:
    sage -python repair-limit-check.py

A private draft requires NUMBERDB_API_KEY or the key on stdin. This reads T606,
checks its original generator and independent integer moment recurrences, and
corroborates the limiting-norm claim. The accompanying repair-verification.md
contains the proof by the central limit theorem and uniform integrability;
finite numerical checks alone do not prove convergence.
"""
import sys,json,os,re,math,hashlib,platform
from pathlib import Path
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 sage.env import SAGE_VERSION
from generate import UpperKhintchine

key=os.environ.get('NUMBERDB_API_KEY') or sys.stdin.read().strip()
if re.match(r'^[A-Z_][A-Z0-9_]*\s*=',key):
    key=key.split('=',1)[1].strip().strip('\"\'')
numberdb.configure(api_key=key)
doc=numberdb.table('T606')
gen=UpperKhintchine()
rows=doc['Numbers']
assert isinstance(rows,dict)
assert list(rows)==[v['p'] for v in gen.enumerate()]

def decimal_exact(s):
    s=str(s).lower(); mant,sep,exp=s.partition('e'); exp=int(exp) if sep else 0
    whole,dot,frac=mant.partition('.')
    return QQ(int(whole+frac))*QQ(10)**(exp-len(frac)),QQ(10)**(exp-len(frac))

def stored(s):
    if isinstance(s,dict):s=s['number']
    return decimal_exact(s)

R=RealBallField(700)
max_radius=R(0)
checks=0
independent=0
for ps,row in rows.items():
    p=QQ(ps)
    value=gen.value({'p':ps},100)
    if p<=2:
        assert value==1 and str(row.get('number') if isinstance(row,dict) else row)=='1'
        continue
    center,ulp=stored(row)
    assert value.is_finite() and R(R(value).rad())<R(10)**-120
    assert R(value)>R(center-ulp) and R(value)<R(center+ulp)
    checks+=1
    if p.denominator()==1:
        order=int(p)
        if order%2==0:
            moment=R(math.prod(range(1,order,2)))
        else:
            moment=R(2/R.pi()).sqrt()*math.prod(range(2,order,2))
        v=(moment.log()/R(p)).exp()
        assert v.is_finite() and R(v.rad())<R(10)**-180
        assert v>R(center-ulp) and v<R(center+ulp)
        independent+=1

# Exact coefficientwise comparison: 1/(2k)! <= 1/(2^k k!).
# For a given p, the proof uses any integer m with 2m>p; 2m=34 covers
# every currently selected p. Multiplication preserves coefficient ordering.
assert max(QQ(x) for x in rows)<34
for k in range(18):
    assert QQ(1)/math.factorial(2*k)<=QQ(1)/(2**k*math.factorial(k))

# Gaussian moments by integration by parts, corroborated against the bound.
# Symbolic binomial identities for the two known finite-sum controls.
P=PolynomialRing(QQ,'n');n=P.gen()
assert n/n==1
assert (n+3*n*(n-1))==3*n*n-2*n

# Exact binomial probabilities; all roots/powers and uncertainty propagated
# through Arb. The normalization is S_n=(2K-n)/sqrt(n).
exponents=[QQ(13)/6,QQ(3),QQ(4),QQ(19)/2,QQ(32)]
limits={str(p):R(2).sqrt()*((R((p+1)/2).gamma()/R.pi().sqrt()).log()/R(p)).exp() for p in exponents}
results=[];last={}
for N in [16,64,256,1024]:
    probs=[QQ(math.comb(N,k))/2**N for k in range(N+1)]
    assert sum(probs)==1
    exact2=sum(prob*QQ((2*k-N)**2)/N for k,prob in enumerate(probs))
    exact4=sum(prob*QQ((2*k-N)**4)/N**2 for k,prob in enumerate(probs))
    assert exact2==1 and exact4==3-QQ(2)/N
    xs=[R(abs(2*k-N))/R(N).sqrt() for k in range(N+1)]
    for p in exponents:
        moment=R(0)
        for x,prob in zip(xs,probs):
            if not x.is_zero():moment+=R(prob)*(x.log()*R(p)).exp()
        norm=(moment.log()/R(p)).exp();limit=limits[str(p)];gap=limit-norm
        assert norm.is_finite() and gap>0 and R(norm.rad())<R(10)**-180
        if str(p) in last:assert gap<last[str(p)]
        last[str(p)]=gap
        results.append({'n':N,'p':str(p),'norm':str(norm),'gaussian_norm':str(limit),'gap':str(gap)})
print(json.dumps({'table':'T606','skill_manifest':'agents/runs/20261029T095925Z-repair-skills/manifest.json','skill_sha256':{'local':'28ecbeeccb0c5e5d6eb541b5a5b0f88aec47b154374d909085116d34f5629199','live':'28ecbeeccb0c5e5d6eb541b5a5b0f88aec47b154374d909085116d34f5629199'},'versions':{'Sage':SAGE_VERSION,'numberdb':numberdb.__version__,'Python':platform.python_version()},'generator_decimal_enclosures':checks,'independent_integer_enclosures':independent,'exact_unit_entries':2,'coefficientwise_bounds_checked_through':34,'finite_sum_controls':'E S_n^2=1 and E S_n^4=3-2/n at n=16,64,256,1024','finite_sum_norm_checks':results,'limitations':'Finite norm comparisons corroborate but do not prove convergence. See repair-verification.md for the proof. No values were changed.'},indent=2))