checks.py

back to table · edit · history · where entries came from · files · download

8406 bytes, as of the version from 2026-10-26 23:00. Recorded here, not run.

"""Independent numerical checks for Airy zeta values, numberdb.org/T531.

Run: sage -python checks.py. With NUMBERDB_KEY_FROM_STDIN=1, read the
private draft and repeat the checks on the intervals its written digits claim.
No writes. The checked selection is n=2,...,20; there are no generator imports
in the Taylor-quotient, source-formula or spectral-sum implementations.
"""
import json
import os
import sys
from decimal import Decimal
from pathlib import Path

import numberdb.sage as numberdb
import numberdb as core
from numberdb._write import to_text
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
from sage.version import version as sage_version
import mpmath as mp

from generate import AiryZeta, evaluate, zeta_polynomials
from check import exactness, measure, names_its_rings, agrees_with


def require_ball(x, digits):
    assert x.is_finite() and x > 0
    B = RealBallField(700)
    assert B(x.rad()) < B(abs(x.center())) * B(10)**(-digits)


def overlap(a, b, digits=100):
    require_ball(a, digits)
    require_ball(b, digits)
    assert a.overlaps(b), (str(a), str(b))


def taylor_quotient(max_order, precision):
    """Ai Taylor coefficients followed by long division, not Riccati."""
    C = ComplexBallField(precision)
    B = RealBallField(precision)
    ai0 = C(0).airy_ai()
    aip0 = C(0).airy_ai_prime()
    assert ai0.imag().is_zero() and aip0.imag().is_zero()
    c = [ai0.real(), aip0.real(), B(0)]
    for j in range(3, max_order+1):
        c.append(c[j-3] / B(j*(j-1)))
    q = []
    for m in range(max_order):
        rhs = (m+1)*c[m+1]
        for j in range(1, m+1):
            rhs -= c[j]*q[m-j]
        q.append(rhs/c[0])
    return {n: (-1)**(n-1)*q[n-1] for n in range(2, max_order+1)}


def published_formulas(precision):
    """Crandall 4.12-4.16, divided by 2^n, and Voros 42.2/43.3."""
    B=RealBallField(precision)
    # Reflection-form kappa is separately evaluated, avoiding Gamma(1/3).
    r=B(3)**(QQ(5)/6)*B(QQ(2)/3).gamma()**2/(2*B.pi())
    return {2:r*r, 3:B(QQ(1)/2)-r**3, 4:r**4-r/3,
            5:B(QQ(5)/12)*r*r-r**5,
            6:r**6-r**3/2+B(QQ(1)/20)}


def decimal_rational(text):
    d=Decimal(text)
    sign, digits, exponent=d.as_tuple()
    num=int(''.join(map(str,digits))) * (-1 if sign else 1)
    return QQ(num) * QQ(10)**exponent


def written_ball(text, precision=700):
    """A decimal means its midpoint plus/minus one final-place unit."""
    d=Decimal(text)
    center=decimal_rational(text)
    unit=QQ(10)**d.as_tuple().exponent
    return RealBallField(precision)(center, unit)


def assert_decimal_encloses(text, ball):
    center=decimal_rational(text)
    unit=QQ(10)**Decimal(text).as_tuple().exponent
    assert ball > center-unit and ball < center+unit


def spectral_bounds(count=64):
    """Independent Airy-root sanity check, with an explicit infinite tail.

    Dirichlet-Neumann bracketing at x=L>lambda gives
    N(lambda) <= L*sqrt(lambda)/pi + 1/2. Let L decrease to lambda:
    lambda_k >= [pi*(k-1/2)]^(2/3). Integral comparison then bounds
    sum_{k>K} lambda_k^(-n) by
    pi^(-2n/3)*(K-1/2)^(1-2n/3)/(2n/3-1).

    Roots are located by mpmath and individually bracketed with Arb signs.
    Their index assignment is numerical, not a certified zero count. Thus
    this is independent validation, not the source of the 100-digit proof.
    """
    C=ComplexBallField(300)
    B=RealBallField(300)
    # A nonzero known-answer control for the same Airy evaluator.
    exact_ai0=1/(B(3)**(QQ(2)/3)*B(QQ(2)/3).gamma())
    control=C(0).airy_ai().real()
    overlap(control,exact_ai0,digits=80)
    mp.mp.dps=70
    roots=[]
    eps=QQ(1)/QQ(10)**55
    for j in range(1,count+1):
        center=decimal_rational(mp.nstr(mp.airyaizero(j),65))
        lo=C(center-eps).airy_ai().real()
        hi=C(center+eps).airy_ai().real()
        assert lo.is_finite() and hi.is_finite() and lo*hi < 0
        roots.append(-B(center,eps))
    assert all(roots[i+1] > roots[i] for i in range(count-1))
    bounds={}
    for n in range(2,21):
        partial=sum((r**(-n) for r in roots),B(0))
        power=QQ(2*n)/3
        tail=B.pi()**(-power)*B(QQ(2*count-1)/2)**(1-power)/B(power-1)
        bounds[n]=(partial,tail)
    return bounds, str(control)


def check_relations(values):
    # Positivity, strict decrease, and strict log convexity of this spectrum.
    for n,v in values.items(): require_ball(v,95)
    assert all(values[n]>values[n+1] for n in range(2,20))
    assert all(values[n]*values[n+2]>values[n+1]**2 for n in range(2,19))
    # Riccati recurrence evaluated on these actual values. q_0=-kappa is
    # independently obtained from Arb's Airy implementation.
    C=ComplexBallField(700)
    q={0:(C(0).airy_ai_prime()/C(0).airy_ai()).real()}
    q.update({n-1:(-1)**(n-1)*v for n,v in values.items()})
    B=RealBallField(700)
    for m in range(19):
        residual=(m+1)*q[m+1]+sum((q[j]*q[m-j] for j in range(m+1)),B(0))-(1 if m==1 else 0)
        assert residual.is_finite() and residual.contains_zero()
        assert B(residual.rad())<B(10)**-90


def main():
    stored=None
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        token=sys.stdin.read().strip()
        if token.startswith('NUMBERDB_') and '=' in token:
            token=token.split('=',1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY']=token
        stored=core.table('T531')
    g=AiryZeta()
    generated={int(p['n']):g.value(p,100) for p in g.enumerate()}
    assert not exactness(generated)
    assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
    assert not exactness(zeta_polynomials(20))
    # Control: polynomial arithmetic must reproduce the known n=2 identity.
    P=zeta_polynomials(20)
    assert P[2]==P[2].parent().gen()**2
    higher={n:evaluate(n,629) for n in generated}
    independent=taylor_quotient(20,700)
    source=published_formulas(700)
    source_at_low={n:P[n](QQ(0)) for n in range(2,7)}
    assert not agrees_with(source_at_low,lambda n:{2:QQ(0),3:QQ(1)/2,4:QQ(0),5:QQ(0),6:QQ(1)/20}[n])
    rows={}
    for n,v in generated.items():
        overlap(v,higher[n],115)
        overlap(v,independent[n],115)
        if n in source: overlap(v,source[n],115)
        text=to_text(v,digits=100)
        written=written_ball(text)
        assert_decimal_encloses(text,independent[n])
        rows[str(n)]={'number':text,'ball':str(v),'accuracy_bits':int(v.accuracy()),
                      'relative_radius':str(v.rad()/abs(v.center()))}
    check_relations(generated)
    bounds,control=spectral_bounds()
    tail_summary={}
    for n,v in generated.items():
        partial,tail=bounds[n]
        assert v>partial and v<partial+tail
        tail_summary[str(n)]={'partial':str(partial),'tail_bound':str(tail)}
    if stored is not None:
        raw=stored['Numbers']
        if isinstance(raw,dict):
            texts={int(n):r['number'] for n,r in raw.items()}
        else:
            texts={int(r['params']['n']):r['number'] for r in raw}
        read={n:written_ball(t) for n,t in texts.items()}
        assert list(read)==list(range(2,21))
        check_relations(read)
        for n,v in read.items():
            assert_decimal_encloses(texts[n],independent[n])
            if n in source: assert_decimal_encloses(texts[n],source[n])
            partial,tail=bounds[n]
            assert v>partial and v<partial+tail
        report=g.verify(sample=None)
        assert report.ok
        print('STORED VERIFY',report,flush=True)
    payload={'table':'T531','skill_manifest':'agents/runs/20261026T224836Z-build-skills/manifest.json',
      'skill_sha256':'b4d51ded0d52cfe3dcc0859de0845c3879b0c2363273ba7e809a7e06ded40b18',
      'versions':{'sage':sage_version,'numberdb':core.__version__,'mpmath':mp.__version__},
      'working_bits':numberdb.bits(100,losing=96),'rerun_bits':629,'quotient_bits':700,
      'minimum_accuracy_bits':min(r['accuracy_bits'] for r in rows.values()),
      'measurement':measure(generated),'serialized_measurement':measure({n:r['number'] for n,r in rows.items()}),
      'all_19_passed':True,'stored_checked':stored is not None,'known_answer_control':control,
      'entries':rows,'root_partial_sums_and_tails':tail_summary}
    print('EVIDENCE_JSON '+json.dumps(payload),flush=True)


if __name__=='__main__':
    main()