checks.py

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

10598 bytes, as of the version from 2026-10-28 15:38 (current). Recorded here, not run.

"""Reproduce T584 build checks without writing to NumberDB.

sage -python checks.py
NUMBERDB_CHECK_STORED=1 sage -python checks.py  # authenticated for a draft

The certified computation is independently compared with mpmath roots and
siegelz derivatives at 125 decimal digits. Every output is also recomputed at
120 requested digits in FLINT. Counts with multiplicity certify the numbering;
nonzero first derivatives certify simplicity. Additional controls, symmetries,
zeta derivative identities, source values and corpus checks are recorded.
Prints one EVIDENCE_JSON record; no file is silently overwritten.
"""
import os
import sys
import json
import time
import hashlib
import platform
from decimal import Decimal
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.complex_arb import ComplexBallField
from sage.env import SAGE_VERSION
from numberdb._write import to_text
from numberdb._generate import _stored_entries
import mpmath as mp
from generate import HardyZDerivatives, compute_zero, bridge, ZERO_INDICES, DERIVATIVE_ORDERS
from check import exactness, measure, names_its_rings


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


def decimal_interval(text):
    if '+/-' in text:
        center, radius = map(decimal_q, text.split('+/-'))
        return center-radius, center+radius
    d = Decimal(text)
    center = decimal_q(text)
    ulp = QQ(10)**int(d.as_tuple().exponent)
    return center - ulp, center + ulp


def inside_decimal(ball, text):
    lo, hi = decimal_interval(text)
    return bool(ball > lo and ball < hi)


def finite_complex(z):
    return z.real().is_finite() and z.imag().is_finite()


def main():
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        key = sys.stdin.read().strip()
        if '=' in key: key = key.split('=', 1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY'] = key
    start = time.monotonic()
    folder = Path(__file__).parent
    b = bridge()
    R = RealBallField(493)
    C = ComplexBallField(493)
    mp.mp.dps = 125
    manifest = json.loads((folder/'skill-manifest.json').read_text())
    assert manifest['status'] == 'ready'
    assert len({x['sha256'] for x in manifest['versions']}) == 1
    evidence = {
        'table':'T584', 'skill_manifest': 'agents/runs/20261028T151417Z-build-skills/manifest.json',
        'skill_versions':manifest['versions'],
        'versions':{'python':platform.python_version(),'sage':SAGE_VERSION,
                    'numberdb':numberdb.__version__,'mpmath':mp.__version__,'flint':b.flint_version()},
        'working_bits':493,'working_guard_bits':160,'stored_significant_digits':100,
        'independent_mpmath_dps':125,
        'n_groups':{'discrete_index':[str(n) for n in range(1,101)],
                    'documented_application':['6709','6710']},
        'k_values':['1','2','3'], 'controls':[], 'rows':[], 'counts':[], 'checks':{}}
    # Nonzero control for both implementations; the odd terms are zero at origin.
    control = b.jet(R(0), 4, 493)
    assert control[0].real().overlaps(R(QQ(1)/2).zeta())
    assert control[0].real() < -1 and control[0].real() > -2
    assert control[1].contains_zero() and control[3].contains_zero()
    assert abs(mp.siegelz(0) - mp.zeta(mp.mpf('0.5'))) < mp.mpf('1e-120')
    assert b.count(0,128) == 0 and b.count(20,128) == 1
    evidence['controls'] += ['Z(0)=zeta(1/2) in (-2,-1), both implementations',
                             'odd derivatives vanish at zero', 'N(0)=0 and N(20)=1']
    # Distinct critical-line roots account for the entire zero count up to
    # 7006, proving the critical-line indexing even before the distant pair.
    print('certifying all 6710 preceding roots at 96 bits',file=sys.stderr,flush=True)
    prefix = b.zero_prefix(6710,96)
    previous = R(0)
    min_gap = R(100)
    for root in prefix:
        assert finite_complex(root) and root.real() == QQ(1)/2
        assert root.imag() > previous
        min_gap = min(min_gap,root.imag()-previous)
        previous = root.imag()
    assert previous < 7006 and b.count(7006,128) == 6710
    evidence['finite_indexing_proof'] = {
        'certified_distinct_critical_line_roots':6710, 'precision_bits':96,
        'count_at_7006':6710, 'first':str(prefix[0]),'last':str(prefix[-1]),
        'minimum_gap_ball':str(min_gap),
        'root_list_sha256':hashlib.sha256('\n'.join(map(str,prefix)).encode()).hexdigest(),
        'conclusion':'All nontrivial zeros with 0 < Im(s) <= 7006 are simple and on the critical line.'}
    print('finite indexing proof complete',file=sys.stderr,flush=True)
    values = {}; higher = {}; texts = {}; gammas = {}
    stored = None
    if os.environ.get('NUMBERDB_CHECK_STORED') == '1':
        stored = _stored_entries(numberdb.table('T584'))
        assert list(stored) == [f'{n},{k}' for n in ZERO_INDICES for k in DERIVATIVE_ORDERS]
    t3 = _stored_entries(numberdb.table('T3'))
    mp_worst = mp.mpf(0)
    min_accuracy = 100000
    paper = {34:['-1.590846','3.8401','1.4834'],35:['1.361151','2.2575','-4.6657'],
             6709:['0.414558','-21.3028','-57.1590'],6710:['-0.427037','-23.2882','-47.9737']}
    for n in ZERO_INDICES:
        gamma, derivatives = compute_zero(n,100)
        gh, dh = compute_zero(n,120)
        gammas[n] = gamma
        assert gamma.overlaps(gh) and gamma.overlaps(prefix[n-1].imag())
        assert abs(gh-gamma) < R(10)**(-135)
        # Counting on either side is well separated from the root. It proves
        # the ordinal among all nontrivial zeros, counted with multiplicity.
        center = gamma.center().exact_rational()
        left = center - QQ(1)/1000000
        right = center + QQ(1)/1000000
        a, z = b.count(left,128), b.count(right,128)
        assert a == n-1 and z == n
        evidence['counts'].append({'n':n,'left_count':n-1,'right_count':n,
                                  'offset':'1/1000000','gamma_ball':str(gamma)})
        if n <= 100:
            assert inside_decimal(gamma,t3[str(n)])
        mz = mp.zetazero(n)
        assert mp.re(mz) == mp.mpf('0.5')
        mg = mp.im(mz)
        assert abs(decimal_q(mp.nstr(mg,125))-gamma) < R(10)**(-115)
        negative = b.jet(-gamma,4,493)
        # First derivative identity independently through Sage's zeta derivative.
        s = C(QQ(1)/2,gamma)
        theta = C(QQ(1)/4,gamma/2).log_gamma().imag() - gamma*R.pi().log()/2
        rhs = C(0,1)*C(0,theta).exp()*s.zetaderiv(1)
        assert finite_complex(rhs) and rhs.imag().contains_zero()
        assert rhs.real().accuracy() > 400 and rhs.real().overlaps(derivatives[0])
        assert derivatives[0] * (-1)**(n+1) > 0
        for k in DERIVATIVE_ORDERS:
            key = f'{n},{k}'
            value = derivatives[k-1]; high = dh[k-1]
            assert value.overlaps(high)
            assert value.is_finite() and high.is_finite()
            min_accuracy = min(min_accuracy,int(value.accuracy()))
            par = negative[k].real() * (1,1,2,6)[k]
            assert negative[k].imag().contains_zero() and par.is_finite()
            assert par.accuracy() > 400 and par.overlaps((-1)**k * value)
            v = mp.siegelz(mg,derivative=k)
            vq = decimal_q(mp.nstr(v,125))
            error = abs(vq-value)
            assert error < R(10)**(-110)*max(R(1),abs(value))
            mp_error = abs(v-mp.mpf(str(value.center()))) / max(mp.mpf(1),abs(v))
            mp_worst = max(mp_worst,mp_error)
            text = stored[key] if stored is not None else to_text(value,100)
            assert inside_decimal(value,text) and inside_decimal(high,text)
            lo,hi = decimal_interval(text)
            assert lo < vq < hi
            if n in paper: assert inside_decimal(value,paper[n][k-1])
            values[key]=value; higher[key]=high; texts[key]=text
            evidence['rows'].append({'n':n,'k':k,'number':text,
                                     'ball':str(value),'accuracy_bits':int(value.accuracy()),
                                     'mpmath_125dps':mp.nstr(v,125)})
        if n % 20 == 0 or n>100:
            print(f'checked n={n}; {time.monotonic()-start:.1f}s',file=sys.stderr,flush=True)
    # Independent higher precision mpmath controls at both ordinary and close zeros.
    with mp.workdps(155):
        for n in [1,34,35,100,6709,6710]:
            t=mp.im(mp.zetazero(n))
            for k in DERIVATIVE_ORDERS:
                v=decimal_q(mp.nstr(mp.siegelz(t,derivative=k),155))
                lo,hi=decimal_interval(texts[f'{n},{k}'])
                assert lo < v < hi
    assert not inside_decimal(values['1,1']+QQ(1)/10**90,texts['1,1'])
    evidence['controls'].append('adding 1e-90 to the first derivative is rejected by its stored 100-digit interval')
    assert not exactness(values) and not names_its_rings(str(folder/'generate.py'))
    # Supplemental cross-family diagnostic from the paper, not a stored quantity.
    metrics={}
    for n,target in [(34,'0.57'),(6709,'0.007')]:
        fs=[]
        for j in [n,n+1]:
            d1,d2,d3=[values[f'{j},{k}'] for k in DERIVATIVE_ORDERS]
            fs.append(-d3/d1+QQ(3)/4*(d2/d1)**2)
        ghat=(gammas[n+1]-gammas[n])**2*(fs[0]+fs[1])/3-2
        assert ghat.is_finite() and ghat.accuracy()>400
        assert inside_decimal(ghat,target)
        metrics[str(n)]=str(ghat)
    evidence['checks']={'all_306_certified':True,'all_306_mpmath_125dps':True,
       'all_306_higher_precision_recomputed':True,'all_306_decimal_intervals_contain_balls':True,
       'all_306_parity':True,'all_102_zeta_derivative_identities':True,
       'all_102_root_counts_and_simplicity':True,'first_100_roots_match_T3':True,
       'all_12_Simonic_Table_1_values':True,'Simonic_Table_2_ghat':metrics,
       '18_mpmath_155dps_checks':True,'stored_values_read_back':stored is not None}
    evidence['min_accuracy_bits']=min_accuracy
    evidence['max_mpmath_relative_error']=mp.nstr(mp_worst,8)
    evidence['size_ball_representation']=measure(values)
    evidence['size_decimal_representation']=measure(texts)
    evidence['seconds']=time.monotonic()-start
    evidence['source_sha256']={p:hashlib.sha256((folder/p).read_bytes()).hexdigest()
                              for p in ['generate.py','hardy_arb.pyx','checks.py','argument-selection.md','check.py']}
    print('EVIDENCE_JSON '+json.dumps(evidence,sort_keys=True))

if __name__ == '__main__':
    main()