"""Independent arithmetic checks for lambda_g Hodge integrals, numberdb.org/T462.

    sage -python checks.py             # local generated values
    sage -python checks.py --stored    # values read from T462; needs its key

The checker uses Python fractions and a Bernoulli recurrence, independently
of the generator's sine-series inversion and binomial products. It also uses
three published Macaulay2 HodgeIntegrals examples, genus-one boundary
reduction to genus-zero string recursion, the corpus Bernoulli numbers,
and a 512-bit Arb calculation from even zeta values. Files bernoulli.json
and T458.json are API snapshots dated 2026-10-16; the latter is a companion
draft and is used only as an additional check, not as mathematical authority.
"""
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from fractions import Fraction as F
from functools import lru_cache
from itertools import combinations_with_replacement
from math import factorial
from pathlib import Path
import json
import os
import sys
from check import exactness, measure, agrees_with, names_its_rings

ROOT = Path(__file__).resolve().parent


def bernoullis(limit):
    # B_0=1 and sum_{k=0}^{m} binomial(m+1,k) B_k = 0 for m>0.
    from math import comb
    bs = [F(1)]
    for m in range(1, limit+1):
        bs.append(-sum(F(comb(m+1,k))*bs[k] for k in range(m))/F(m+1))
    return bs


def independent(key):
    g, n, ds = key
    bs = bernoullis(2*g)
    if g == 0:
        base = F(1)
    else:
        base = (1-F(1, 2**(2*g-1))) * abs(bs[2*g]) / factorial(2*g)
    denominator = 1
    for d in ds:
        denominator *= factorial(d)
    return base * F(factorial(sum(ds)), denominator)


@lru_cache(None)
def genus_zero(ds):
    # Pure psi intersections from the string equation alone; no multinomial.
    if len(ds) < 3 or min(ds) < 0 or sum(ds) != len(ds)-3:
        return F(0)
    if ds == (0,0,0):
        return F(1)
    rest = list(ds)
    rest.remove(0)
    answer = F(0)
    for j, d in enumerate(rest):
        if d:
            lowered = rest.copy()
            lowered[j] -= 1
            answer += genus_zero(tuple(sorted(lowered)))
    return answer


def scalar(value):
    return value['number'] if isinstance(value,dict) else value


def flatten(doc):
    return {(int(g),int(n),tuple(map(int,d.split(',')))): QQ(scalar(value))
            for g, ns in doc['Numbers'].items() for n, profiles in ns.items()
            for d, value in profiles.items()}


def validate(vals):
    # Controls are evaluated before comparisons with the proposed data.
    bs = bernoullis(10)
    assert bs[:5] == [F(1),F(-1,2),F(1,6),F(0),F(-1,30)]
    assert genus_zero((0,0,0)) == 1
    assert genus_zero((0,0,0,1,1)) == 2
    assert genus_zero((0,0,1)) == 0
    assert independent((1,1,(0,))) == F(1,24)
    assert independent((2,2,(1,2))) == F(7,1920)
    assert not agrees_with({(1,1,(0,)):F(1,24)},independent)
    assert agrees_with({(1,1,(0,)):F(1,12)},independent)
    assert exactness({'float-control':0.5})

    expected = set()
    for g in range(1,6):
        for n in range(1,6):
            total = 2*g-3+n
            for ds in combinations_with_replacement(range(total+1),n):
                if sum(ds) == total:
                    expected.add((g,n,ds))
    assert set(vals) == expected
    assert not exactness(vals)
    measurement = measure(vals)
    vals = {k:F(int(v.numerator()), int(v.denominator())) for k,v in vals.items()}
    assert not names_its_rings(str(ROOT/'generate.py'))
    disagreements = agrees_with(vals, independent)
    assert not disagreements, disagreements[:5]
    for (g,n,ds), value in vals.items():
        assert len(ds) == n and ds == tuple(sorted(ds)) and min(ds) >= 0
        assert 2*g-2+n > 0 and sum(ds)+g == 3*g-3+n
        assert value > 0

    # Published outputs from the independent Macaulay2 HodgeIntegrals package.
    examples = {(1,1,(0,)):F(1,24), (2,2,(1,2)):F(7,1920),
                (3,3,(1,2,3)):F(31,16128)}
    for key, value in examples.items():
        assert vals[key] == value

    # b_g values tested against an external Bernoulli table, and via zeta.
    corpus = json.loads((ROOT/'bernoulli.json').read_text())['Numbers']
    R = RealBallField(512)
    def zeta_base(g):
        return 2*(1-R(2)**(1-2*g))*R(2*g).zeta()/(2*R.pi())**(2*g)
    control = zeta_base(1)
    assert control.is_finite() and R(control.rad()) < R(10)**(-140)
    assert control.overlaps(R(QQ(1)/QQ(24)))
    assert not control.overlaps(R(QQ(1)/QQ(12)))
    zeta_radii = []
    b = [QQ(1)]
    for g in range(1,6):
        B = F(str(scalar(corpus[str(2*g)])))
        assert B == bs[2*g]
        base = vals[(g,1,(2*g-2,))]
        assert base == (1-F(1,2**(2*g-1)))*abs(B)/factorial(2*g)
        ball = zeta_base(g)
        assert ball.is_finite() and R(ball.rad()) < R(10)**(-140)
        assert ball.overlaps(R(QQ(base.numerator)/QQ(base.denominator)))
        zeta_radii.append(str(ball.rad()))
        b.append(QQ(base.numerator)/QQ(base.denominator))
    for g in range(1,6):
        convolution = sum(QQ((-1)**k)/QQ(2**(2*k)*factorial(2*k+1))*b[g-k]
                          for k in range(g+1))
        assert convolution == 0

    # Genus one: lambda_1=delta_irr/12, with gluing degree two.
    # This reduces to 1/24 of a genus-zero pure-psi intersection with two
    # extra zero exponents. The genus-zero recursion above is independent.
    genus_one = 0
    companion = 0
    psi = json.loads((ROOT/'T458.json').read_text())['Numbers']['0']
    for (g,n,ds), value in vals.items():
        if g != 1:
            continue
        zero_profile = tuple(sorted(ds+(0,0)))
        assert value == genus_zero(zero_profile)/24
        genus_one += 1
        key = ','.join(map(str,zero_profile))
        if str(n+2) in psi and key in psi[str(n+2)]:
            assert value == F(str(scalar(psi[str(n+2)][key])))/24
            companion += 1

    string_count = dilaton_count = 0
    for (g,n,ds), value in vals.items():
        if n <= 1:
            continue
        if ds[0] == 0:
            rest = list(ds[1:])
            rhs = F(0)
            for j,d in enumerate(rest):
                if d:
                    lowered = rest.copy()
                    lowered[j] -= 1
                    rhs += vals[(g,n-1,tuple(sorted(lowered)))]
            assert value == rhs
            string_count += 1
        if 1 in ds:
            rest = list(ds)
            rest.remove(1)
            assert value == (2*g-3+n)*vals[(g,n-1,tuple(rest))]
            dilaton_count += 1

    return {'entries':len(vals), 'exactness':'pass', 'range':'all g=1..5, n=1..5, sorted top-degree profiles',
            'independent_bernoulli_formula_matches':len(vals),
            'published_macaulay2_examples':len(examples),
            'corpus_bernoulli_matches':5, 'arb_zeta_matches':5,
            'arb_zeta_working_bits':512, 'arb_zeta_absolute_radii':zeta_radii,
            'sine_series_coefficients':5, 'genus_one_boundary_checks':genus_one,
            'T458_genus_zero_matches':companion,
            'string_equations':string_count, 'dilaton_equations':dilaton_count,
            'controls':'known rational, nonzero genus-zero, wrong-degree, wrong normalization and float controls passed',
            'measurement':measurement}


if __name__ == '__main__':
    stored = '--stored' in sys.argv or os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1'
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
    if stored:
        vals = flatten(numberdb.table('T462'))
    else:
        from generate import LambdaGHodge
        gen = LambdaGHodge()
        vals = {(int(p['g']),int(p['n']),tuple(map(int,p['d'].split(',')))):
                scalar(gen.value(p,100)) for p in gen.enumerate()}
    result = validate(vals)
    result['source'] = 'database T462' if stored else 'local generator'
    print(json.dumps(result,indent=2,default=str))