checks.py

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

10630 bytes, as of the version from 2026-10-27 21:27 (current). Recorded here, not run.

"""Numerical checks for Area Siegel–Veech constants -- numberdb.org/T561.

Download generate.py, check.py and table_records.py alongside this file, then run:
    sage -python checks.py
To check a private stored draft, use a NumberDB client with its author's key:
    check_stored(numberdb.table('T561', client=client))

The EMZ checks use separately transcribed configuration constants and the number
of cylinders, not EKZ sums. EKZ itself used EMZ, so these are independent code
paths and source consistency checks, not an independent proof of those values.
The Chen–Möller checks use their geometric divisor slopes, independently of EMZ
cylinder counting. Hyperelliptic checks use the genus-zero quadratic differential
formula and the orienting double cover, independently of the tabulated sums.
"""
import json
from fractions import Fraction
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.version import version as sage_version
from sage.rings.rational_field import QQ
from sage.rings.real_mpfi import RealIntervalField
from sage.rings.real_arb import RealBallField
from generate import AreaSiegelVeech, SUMS, SOURCE_ROWS, coefficient, area_constant
from check import exactness, measure, names_its_rings
from table_records import to_records

# Each pair is (number of cylinders, configuration constant times zeta(2)).
# Zero-cylinder configurations contribute zero and are omitted explicitly.
# EMZ v2 Examples 14.12–14.13, pp. 86–87 and Appendix B.2, pp. 98–100.
EMZ_CYLINDERS = {
    (3, (4,), 'odd'): [(1,'81/32'), (1,'15/32')],
    (3, (2,2), 'odd'): [(1,'6/5'), (1,'32/15'), (2,'1/6')],
    (4, (6,), 'hyp'): [(1,'6')],
    (4, (6,), 'even'): [(1,'10125/8192'), (1,'875/256'), (1,'2835/8192')],
    (4, (6,), 'odd'): [(1,'175/81'), (1,'25/32'), (1,'175/162'),
                         (1,'7/32'), (1,'7/32'), (1,'35/864')],
    (4, (3,3), 'hyp'): [(1,'15/2')],
    (4, (3,3), 'nonhyp'): [(1,'3699/1120'), (1,'64/35'),
                            (1,'27/80'), (1,'1/32')],
    (4, (4,2), 'even'): [(1,'4/21'), (1,'256/315'), (1,'8/9'),
                          (1,'32768/8505'), (2,'2/45'), (1,'256/1215')],
    (4, (4,2), 'odd'): [(1,'7/18'), (1,'21/16'), (1,'3/10'),
                         (1,'128/45'), (1,'7/80'), (2,'21/320'),
                         (1,'7/45'), (2,'7/576')],
    (4, (2,2,2), 'even'): [(1,'45/74'), (1,'225/37'), (2,'5/37')],
    (4, (2,2,2), 'odd'): [(1,'63/62'), (1,'144/31'), (2,'9/155'),
                           (2,'16/155'), (3,'1/186')],
}

# Chen–Möller §§5.1–5.2: the hyperelliptic divisor has slope 9 and is
# disjoint from these Teichmüller curves. §6.1: C.(30 lambda+60 psi-4 delta)=0,
# psi=(lambda-delta/12)/4, so delta/lambda=60/7. Their Proposition 4.5 gives
# L=12*kappa/(12-slope); their area constant is pi**2/3 times ours.
CM_DIVISOR_SLOPES = {
    (3, (4,), 'odd'): QQ(9),
    (3, (3,1), 'connected'): QQ(9),
    (4, (6,), 'even'): QQ(60)/QQ(7),
}


def kappa_independent(mu):
    # m(m+2)/(m+1) = (m+1)-1/(m+1).
    return sum((QQ(m+1) - QQ(1)/QQ(m+1) for m in mu), QQ(0))/QQ(12)


def expected_emz(key):
    g, mu, _ = key
    return QQ(6)*sum((QQ(q)*QQ(c) for q,c in EMZ_CYLINDERS[key]), QQ(0))/QQ(2*g+len(mu)-2)


def expected_cm(key):
    kappa=kappa_independent(key[1]); slope=CM_DIVISOR_SLOPES[key]
    total=QQ(12)*kappa/(QQ(12)-slope)
    return QQ(3)*(total-kappa)


def expected_hyperelliptic(g, mu):
    # EKZ Theorem 3(a) and Lemma 1.1: twice the genus-zero area constant.
    if mu == (2*g-2,):
        quadratic_orders=(2*g-3,) + (-1,)*(2*g+1)
    elif mu == (g-1,g-1):
        quadratic_orders=(2*g-2,) + (-1,)*(2*g+2)
    else:
        raise ValueError('Not a hyperelliptic component profile.')
    assert sum(quadratic_orders)==-4
    return -sum((QQ(d*(d+4))/QQ(d+2) for d in quadratic_orders),QQ(0))/QQ(4)


def rational_from_decimal(text):
    f=Fraction(text)
    return QQ(f.numerator)/QQ(f.denominator)


def written_interval(text):
    """Exact interval claimed by a plain NumberDB decimal (one last-place unit)."""
    text=str(text)
    assert '.' in text
    mantissa, _, exponent=text.lower().partition('e')
    places=len(mantissa.partition('.')[2])
    power=(int(exponent) if exponent else 0)-places
    center=rational_from_decimal(text)
    ulp=QQ(10)**power
    return center-ulp, center+ulp


def controls():
    field=RealIntervalField(512)
    third=field(1)/field(3)
    assert QQ(third.lower()) <= QQ(1)/QQ(3) <= QQ(third.upper())
    assert QQ(third.upper())-QQ(third.lower()) < QQ(2)**(-500)
    assert written_interval('0.125') == (QQ(124)/1000,QQ(126)/1000)
    assert expected_hyperelliptic(2,(2,))==QQ(10)/3
    assert expected_hyperelliptic(2,(1,1))==QQ(15)/4
    assert expected_emz((3,(4,),'odd'))==QQ(18)/5
    # Independent exact solution of the CM divisor equation.
    slope=QQ(60)/7
    assert QQ(30)+QQ(60)*(QQ(1)-slope/12)/4-QQ(4)*slope==0
    return 'All exact, interval, and divisor controls passed.'


def check_written(key, text):
    g, mu, component=key
    lo, hi=written_interval(text)
    assert 0 < lo < hi
    assert (hi-lo)/lo < QQ(10)**(-98)
    field=RealIntervalField(512)
    # A separately written MPFI calculation from the source sum.
    q=QQ(3)*(SUMS[key]-kappa_independent(mu))
    independent=field(q)/field.pi()**2
    assert lo <= QQ(independent.lower()) <= QQ(independent.upper()) <= hi
    # The EKZ identity is also checked on the whole written uncertainty interval.
    reconstructed=kappa_independent(mu) + field.pi()**2*field(lo,hi)/field(3)
    assert QQ(reconstructed.lower()) <= SUMS[key] <= QQ(reconstructed.upper())
    if component=='hyp':
        expected=field(expected_hyperelliptic(g,mu))/field.pi()**2
        assert lo <= QQ(expected.lower()) <= QQ(expected.upper()) <= hi
    if key in EMZ_CYLINDERS:
        expected=field(expected_emz(key))/field.pi()**2
        assert lo <= QQ(expected.lower()) <= QQ(expected.upper()) <= hi
    if key in CM_DIVISOR_SLOPES:
        expected=field(expected_cm(key))/field.pi()**2
        assert lo <= QQ(expected.lower()) <= QQ(expected.upper()) <= hi


def check_stored(document):
    """Check values read from the API, including every displayed formula."""
    controls()
    namespace={}
    exec(document['Programs']['sage-hyperelliptic']['code'],namespace)
    entries=to_records(document)
    assert isinstance(entries,list) and len(entries)==25
    generator=AreaSiegelVeech()
    expected_params=[{k:str(v) for k,v in p.items()} for p in generator.enumerate()]
    actual=[]
    for row in entries:
        params=row['params'];actual.append({k:str(v) for k,v in params.items()})
        key=(int(params['g']),tuple(map(int,str(params['mu']).split(','))),params['component'])
        check_written(key,row['number'])
        if key[2]=='hyp':
            ball=namespace['hyperelliptic_area'](key[0],zeros=len(key[1]),digits=100)
            assert ball.is_finite() and relative_bits(ball)>=330
            lo,hi=written_interval(row['number'])
            assert lo<=QQ(ball.lower())<=QQ(ball.upper())<=hi
    assert actual==expected_params
    return {'stored_entries':25,'all_written_intervals_contain_independent_mpfi':True,
            'ekz_identity':25,'hyperelliptic_formula':6,'emz_cylinder_checks':11,
            'cm_divisor_checks':3,'enumeration_matches':True}


def relative_bits(ball):
    lo,hi=QQ(ball.lower()),QQ(ball.upper())
    ratio=lo/(hi-lo)
    bits=int(ratio.numerator().nbits()-ratio.denominator().nbits())
    return bits if QQ(2)**bits<=ratio else bits-1


def run():
    import importlib.metadata
    from pathlib import Path
    assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
    generator=AreaSiegelVeech(); rows=list(generator.enumerate())
    actual=[(int(p['g']),tuple(map(int,p['mu'].split(','))),p['component']) for p in rows]
    assert len(actual)==len(set(actual))==len(SUMS)==25
    assert set(actual)==set(SUMS)
    assert len({(g,mu) for g,mu,c in actual})==18
    assert [sum(g==j for g,mu,c in actual) for j in (2,3,4)]==[2,7,16]
    values={};serialized=[];results=[];accuracy=[]
    control=controls()
    for params,key in zip(rows,actual):
        g,mu,c=key
        q=coefficient(*key)
        assert q>0
        value=generator.value(params,100)
        ball=value['number'] if isinstance(value,dict) else value
        assert ball.is_finite() and ball.lower()>0
        accuracy.append(relative_bits(ball))
        assert accuracy[-1]>=330
        higher=area_constant(*key,digits=160)
        assert higher.is_finite() and relative_bits(higher)>=530
        assert ball.overlaps(higher)
        if c=='hyp':
            assert q==expected_hyperelliptic(g,mu)
            closed = QQ(g*(2*g+1))/QQ(2*g-1) if len(mu)==1 else QQ((g+1)*(2*g+1))/QQ(2*g)
            assert q==closed
        if key in EMZ_CYLINDERS: assert q==expected_emz(key),(key,q,expected_emz(key))
        if key in CM_DIVISOR_SLOPES: assert q==expected_cm(key)
        text=to_text(ball,digits=100)
        check_written(key,text)
        row={'params':params,'number':text}
        if isinstance(value,dict): row.update({k:v for k,v in value.items() if k!='number'})
        serialized.append(row);values[str(key)]=value
        results.append({'g':g,'mu':list(mu),'component':c,'lyapunov_sum':str(SUMS[key]),
                        'pi_squared_times_area_constant':str(q),'number':text,
                        'relative_accuracy_bits':accuracy[-1]})
    assert not exactness(values)
    assert not exactness({str(k):coefficient(*k) for k in actual})
    record={'table':'T561','controls':control,'entries':25,'profiles':18,
            'min_relative_accuracy_bits':min(accuracy),'working_guard_bits':64,
            'digits':100,'rigour':'proven','precision_agreement_entries':25,
            'mpfi_512_bit_containment_entries':25,'hyperelliptic_checks':6,
            'emz_cylinder_checks':11,'cm_divisor_checks':3,
            'distinct_entries_with_additional_mathematical_checks':16,
            'measure':measure(values),'rows':results,
            'skill_manifest':'agents/runs/20261027T205535Z-build-skills/manifest.json',
            'skill_local_sha256':'e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c',
            'skill_live_sha256':'e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c',
            'versions':{'sage':sage_version,'numberdb':getattr(numberdb, '__version__', 'source checkout')}}
    # Stdout permits use of the campaign's read-only Sage wrapper.
    print(json.dumps({'evidence':record,'serialized_entries':serialized},indent=2))


if __name__=='__main__':
    run()