"""Independent checks for numberdb.org/T461.

Run: sage -python checks.py [--stored]
Requires generate.py and the attached check.py beside this file.

The exact comparison uses Python Fraction and the Stirling-number sum for
Bernoulli numbers, followed by the closed Harer-Zagier formula. It does not
call the generator's recurrence or Sage's Bernoulli implementation. A second
comparison uses positive even zeta values and the zeta functional equation
in 512-bit Arb arithmetic. This relates the table to a different family.
The known genus-one group Euler characteristic is an additional control.
"""
import json
import os
import sys
from fractions import Fraction
from functools import lru_cache
from math import factorial
from pathlib import Path

import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from check import agrees_with, exactness, measure, names_its_rings


@lru_cache(None)
def bernoulli_by_partitions(m):
    row = [1]
    for order in range(1, m + 1):
        old = row
        row = [0] + [(old[k - 1] if k - 1 < len(old) else 0)
                     + k * (old[k] if k < len(old) else 0)
                     for k in range(1, order + 1)]
    return sum((Fraction((-1)**k * factorial(k) * row[k], k + 1)
                for k in range(m + 1)), Fraction(0))


def independent_exact(key):
    g, n = key
    if g == 0:
        return Fraction((-1)**(n - 3) * factorial(n - 3))
    return (Fraction((-1)**n * (2*g - 1) * factorial(2*g + n - 3),
                     factorial(2*g)) * bernoulli_by_partitions(2*g))


def check_values(values):
    field = RealBallField(512)
    assert bernoulli_by_partitions(0) == 1
    assert bernoulli_by_partitions(1) == Fraction(-1, 2)
    assert bernoulli_by_partitions(2) == Fraction(1, 6)
    assert bernoulli_by_partitions(4) == Fraction(-1, 30)
    assert field(2).zeta().is_finite()
    assert (field(2).zeta() - field.pi()**2 / 6).contains_zero()
    elliptic = QQ(1)/4 + QQ(1)/6 - QQ(1)/2
    assert elliptic == -QQ(1)/12
    genus_one_ball = -2 * field(2).zeta() / (2*field.pi())**2
    assert (elliptic in genus_one_ball)
    assert field(genus_one_ball.rad()) < field(2)**(-490)
    assert not (QQ(1)/12 in genus_one_ball)
    assert exactness({'deliberate bad float': 1.0})
    assert agrees_with({(1,1): QQ(1)}, lambda k: QQ(independent_exact(k)))

    expected_keys = {(g,n) for n in range(13) for g in range(1,13)
                     if 2*g - 2 + n > 0}
    assert set(values) == expected_keys
    assert len(values) == 155
    assert not exactness(values)
    mismatches = agrees_with(values, lambda k: QQ(independent_exact(k)))
    assert not mismatches, mismatches
    assert values[1,1] == elliptic
    anchors = {(2,0): QQ(-1)/240, (3,0): QQ(1)/1008,
               (2,1): QQ(1)/120, (3,1): QQ(-1)/252,
               (4,1): QQ(1)/240, (5,1): QQ(-1)/132,
               (6,1): QQ(691)/32760}
    for key, value in anchors.items():
        assert values[key] == value

    seeds = {}
    for g in range(1,13):
        z = field(2*g).zeta()
        seed = ((-1)**g * 2 * factorial(2*g-1) * z
                / (2*field.pi())**(2*g))
        assert z.is_finite() and seed.is_finite()
        assert (values[g,1] in seed)
        seeds[g] = seed
    accuracy = []
    for (g,n), value in values.items():
        # This factorial quotient scales the zeta seed without a recurrence.
        factor = (QQ(1) / QQ(2-2*g) if n == 0 else
                  QQ((-1)**(n-1) * factorial(2*g+n-3)) / QQ(factorial(2*g-2)))
        enclosure = seeds[g] * factor
        assert enclosure.is_finite()
        assert (value in enclosure), (g,n,value,enclosure)
        assert field(enclosure.rad()) < field(2)**(-450) * abs(value)
        accuracy.append(next(bits for bits in range(512,449,-1)
                             if field(enclosure.rad()) < field(2)**(-bits)*abs(value)))

    recurrence_count = 0
    for (g,n), value in values.items():
        if (g,n+1) in values:
            assert values[g,n+1] == (2-2*g-n)*value
            recurrence_count += 1
    for n in range(1,13):
        assert values[1,n] == QQ((-1)**n * factorial(n-1)) / QQ(12)
    for g in range(2,13):
        assert values[g,0] == QQ(bernoulli_by_partitions(2*g)) / QQ(2*g*(2*g-2))

    return {'entries':len(values), 'exact_independent_matches':len(values),
            'zeta_ball_matches':len(values), 'zeta_working_bits':512,
            'minimum_zeta_relative_accuracy_bits':min(accuracy),
            'forgetful_recurrence_pairs':recurrence_count,
            'genus_one_group_control':'1/4 + 1/6 - 1/2 = -1/12',
            'negative_controls':'wrong sign and float rejected',
            'named_rings':names_its_rings(str(Path(__file__).with_name('generate.py'))),
            'measurement':measure(values),
            'values':[{ 'g':g,'n':n,'number':str(v)} for (g,n),v in values.items()]}


def main(stored=False):
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        numberdb.configure(api_key=sys.stdin.read().strip())
    if stored:
        doc=numberdb.table('T461')
        rows=doc['Numbers']
        if isinstance(rows, list):
            values={(int(e['params']['g']),int(e['params']['n'])):QQ(e['number']) for e in rows}
            assert len(rows)==len(values)
        else:
            values={(int(g),int(n)):QQ(e['number'] if isinstance(e,dict) else e)
                    for g,group in rows.items() for n,e in group.items()}
    else:
        from generate import OrbifoldEulerCharacteristics, chi
        gen=OrbifoldEulerCharacteristics()
        values={(p['g'],p['n']):gen.value(p,100)['number'] for p in gen.enumerate()}
        for n in range(3,13):
            assert chi(0,n)==QQ(independent_exact((0,n)))
            assert chi(0,n+1)==(2-n)*chi(0,n)
    report=check_values(values)
    report['source']='stored table' if stored else 'generator'
    print(json.dumps(report,indent=2))
    return report


if __name__=='__main__':
    main('--stored' in sys.argv)