"""Independent repair checks for convex-position probabilities, numberdb.org/T468.

With numberdb installed, run `sage -python repair_checks.py` beside generate.py,
T264-document.json, and T468-document-before.json (or T468-document-after.json).
In this campaign run through agents/sage.sh, mounting those files as well.
No API writes and no credentials are needed. The generator is loaded from its
website attachment; the independent formulas are from Marckert's preprint,
https://arxiv.org/pdf/1402.3512, equations (7), (9), and (10).
Equation (7) gives the complement, which is subtracted from 1 here.
"""
import json
from decimal import Decimal
from math import comb, 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 sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing

from generate import ConvexPosition, GUARD_BITS

ROOT = Path(__file__).resolve().parent
AFTER = ROOT / 'T468-document-after.json'
doc = json.loads((AFTER if AFTER.exists() else ROOT / 'T468-document-before.json').read_text())
catalan_doc = json.loads((ROOT / 'T264-document.json').read_text())
R = RealBallField(768)
pi = R.pi()
disk = {
    3: QQ(1),
    4: 1 - 35/(12*pi**2),
    5: 1 - 305/(48*pi**2),
    6: 1 - (146400*pi**2-473473)/(11520*pi**4),
    7: 1 - (512400*pi**2-2900611)/(23040*pi**4),
    8: 1 - (62664108221+1721664000*pi**4-18670881600*pi**2)/(48384000*pi**6),
}


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


def decimal_interval(text):
    d = Decimal(text)
    numerator, denominator = d.as_integer_ratio()
    midpoint = QQ(numerator) / QQ(denominator)
    ulp = QQ(10)**int(d.as_tuple().exponent)
    assert len(d.as_tuple().digits) == 100
    return midpoint, ulp


def dyck_count(m):
    # Counts nonnegative walks by height, independently of the Catalan formula.
    counts = {0: 1}
    for _ in range(2*m):
        next_counts = {}
        for h, count in counts.items():
            for next_h in (h-1, h+1):
                if next_h >= 0:
                    next_counts[next_h] = next_counts.get(next_h, 0) + count
        counts = next_counts
    return counts.get(0, 0)


assert list(doc['Parameters']) == ['region', 'n']
gen = ConvexPosition()
assert gen.digits == 100 and GUARD_BITS == 80
rows = []
exact_count = 0
decimal_count = 0
for params in gen.enumerate():
    region, n = params['region'], int(params['n'])
    text = number(doc['Numbers'][region][str(n)])
    if region == 'equilateral-triangle':
        independent = QQ(2**n * factorial(3*n-3)) / QQ(factorial(n-1)**3 * factorial(2*n))
    elif region == 'unit-square':
        independent = (QQ(comb(2*n-2, n-1)) / QQ(factorial(n)))**2
    else:
        independent = disk[n]
    computed = gen.value(params, gen.digits)['number']
    if '.' not in text:
        assert QQ(text) == independent == computed
        exact_count += 1
        rows.append({'region': region, 'n': n, 'exact': text, 'passed': True})
    else:
        midpoint, ulp = decimal_interval(text)
        assert independent.is_finite() and computed.is_finite()
        assert R(independent.rad()) < R(10)**(-220)
        assert R(computed.rad()) < R(2)*R(10)**(-123)
        assert independent > R(midpoint-ulp) and independent < R(midpoint+ulp)
        assert computed > R(midpoint-ulp) and computed < R(midpoint+ulp)
        assert computed.overlaps(independent)
        difference_ulps = abs((independent - R(midpoint)) / R(ulp))
        assert difference_ulps < R(465)/1000
        decimal_count += 1
        rows.append({'region': region, 'n': n, 'passed': True,
                     'absolute_discrepancy_ulps': str(difference_ulps),
                     'independent_ball_radius': str(independent.rad()),
                     'generator_ball_radius': str(computed.rad())})

P = PolynomialRing(QQ, 'q')
catalan_checks = []
for n in range(3, 9):
    m = n-1
    paths = dyck_count(m)
    assert paths == comb(2*m, m)//(m+1)
    if str(m) in catalan_doc['Numbers']:
        polynomial = P(number(catalan_doc['Numbers'][str(m)]))
        assert polynomial(1) == paths
        source = f'T264 row {m}'
    else:
        assert m == 2
        # T264 explicitly gives this omitted small case in Formulas.
        assert '$C_0(q)=C_1(q)=1$ and $C_2(q)=q+1$.' == catalan_doc['Formulas']['formula-small']
        assert (P.gen()+1)(1) == paths
        source = 'T264 formula-small: C_2(q)=q+1'
    stored = QQ(number(doc['Numbers']['unit-square'][str(n)]))
    assert stored == (QQ(paths) / QQ(factorial(m)))**2
    catalan_checks.append({'n': n, 'catalan': paths, 'source': source, 'passed': True})

assert exact_count == 13 and decimal_count == 5 and len(rows) == 18
print(json.dumps({'table': 'T468', 'checked_document': 'after' if AFTER.exists() else 'before',
                  'exact_entries_passed': exact_count, 'decimal_entries_passed': decimal_count,
                  'working_bits': 768, 'rows': rows, 'catalan_checks': catalan_checks,
                  'scope': 'All stored values, the unchanged attached generator, and the linked Catalan identity. '
                           'Closed forms check transcription and precision; they do not re-prove the disk formulas.'}, indent=2))