#!/usr/bin/env python3
"""Exact finite validation of the new recurrence and bounds, not an infinite proof."""
import argparse
from fractions import Fraction as F
from hashlib import sha256
import importlib.util
from itertools import product
import json
from math import comb, factorial, prod
from pathlib import Path

HERE = Path(__file__).resolve().parent
PREVIOUS = HERE.parent / '0024' / 'check_hardy_sign.py'
spec = importlib.util.spec_from_file_location('cycle24', PREVIOUS)
old = importlib.util.module_from_spec(spec)
spec.loader.exec_module(old)
ZERO = (0, 0, 0, 0)


def multinomial(a):
    return factorial(sum(a)) // prod(factorial(n) for n in a)


def weight(a):
    if a == ZERO:
        return 0
    if all(n % 2 == 0 for n in a):
        return 8 * (multinomial(a) - multinomial(tuple(n // 2 for n in a)))
    if all(n % 2 == 1 for n in a):
        return 8 * multinomial(a)
    return 0


def base_coefficient(a):
    odds = sum(n % 2 for n in a)
    half = tuple(n // 2 for n in a)
    if odds == 0:
        return comb(sum(half) + 2, 2) * multinomial(half)
    if odds == 2:
        return 2 * comb(sum(half) + 3, 3) * multinomial(half)
    return 0


def lower_coefficient(a):
    if all(n % 2 == 1 for n in a):
        half = tuple(n // 2 for n in a)
        return 48 * comb(sum(half) + 2, 2) * multinomial(half)
    return base_coefficient(a)


def within(a):
    return product(*(range(n + 1) for n in a))


def subtract(a, b):
    return tuple(x - y for x, y in zip(a, b))


def compute(degree):
    coefficients = {ZERO: F(1)}
    weights = {a: weight(a) for s in range(1, degree + 1) for a in old.indices(s)}
    assert all(w >= 0 for w in weights.values())
    for s in range(1, degree + 1):
        for a in old.indices(s):
            c = sum((weights.get(b, 0) * coefficients.get(subtract(a, b), 0)
                     for b in within(a) if b != ZERO), F(0)) / s
            assert c >= 0 and c.denominator == 1, (a, c)
            if c:
                coefficients[a] = c
    return coefficients


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--degree', type=int, default=12)
    parser.add_argument('--output', type=Path, default=HERE / 'positive-recurrence-checks.json')
    args = parser.parse_args()
    if args.degree < 4:
        parser.error('--degree must be at least 4 to cover every parity case')
    positive = compute(args.degree)
    # Compare to direct integer rational inversion, not another use of the recurrence.
    r = old.inverse_series(old.denominator(old.FORMS), args.degree)
    numerator = {ZERO: 1}
    one_minus_q = {ZERO: 1, **{tuple(2 * int(i == j) for i in range(4)): -1 for j in range(4)}}
    for _ in range(4):
        numerator = old.mul(numerator, one_minus_q)
    direct_c = old.mul(numerator, r)
    direct_m = old.inverse_series(old.denominator(old.FORMS[1:]), args.degree)
    counts = {'all_indices': 0, 'even_indices': 0, 'odd_indices': 0,
              'C_identity_comparisons': 0, 'MGF_comparisons': 0,
              'interval_integral_comparisons': 0, 'lower_bound_comparisons': 0,
              'axis_equality_comparisons': 0}
    examples = []
    for s in range(args.degree + 1):
        for a in old.indices(s):
            c = positive.get(a, 0)
            assert c == direct_c.get(a, 0), ('C identity', a, c, direct_c.get(a, 0))
            counts['all_indices'] += 1
            counts['C_identity_comparisons'] += 1
            g = sum((base_coefficient(subtract(a, b)) * positive.get(b, 0)
                     for b in within(a)), F(0))
            odd_entries = sum(n % 2 for n in a)
            if s % 2:
                assert g == 0 and direct_m.get(a, 0) == 0
                counts['odd_indices'] += 1
                continue
            counts['even_indices'] += 1
            assert g == (-1)**(odd_entries // 2) * direct_m.get(a, 0), ('MGF', a, g)
            counts['MGF_comparisons'] += 1
            lower = lower_coefficient(a)
            assert g >= lower > 0, ('lower bound', a, g, lower)
            counts['lower_bound_comparisons'] += 1
            scale = F(12 * prod(factorial(n) for n in a), 2**s * factorial(s + 4))
            sign = (-1)**sum(n // 2 for n in a)
            hardy = sign * scale * g
            assert hardy == old.integral_hardy(a), ('interval integral', a, hardy)
            counts['interval_integral_comparisons'] += 1
            if sum(n > 0 for n in a) <= 1:
                N = s // 2
                exact = F((-1)**N * 3, 2 * 4**N * (2*N + 1) * (2*N + 3))
                assert g == lower and hardy == exact
                counts['axis_equality_comparisons'] += 1
            if a in (ZERO, (1, 1, 0, 0), (1, 1, 1, 1), (2, 2, 0, 0), (3, 1, 0, 0), (2, 2, 2, 2)):
                examples.append({'alpha': a, 'C': str(c), 'G': str(g), 'HARDY': str(hardy),
                                 'magnitude_lower_bound': str(scale * lower)})
    assert positive[(1, 1, 1, 1)] == 48
    degree_two_axis = (2, 0, 0, 0)
    lambda_five = old.mul(one_minus_q, direct_c)
    assert lambda_five[degree_two_axis] == -1
    report = {
        'status': 'pass', 'degree': args.degree, 'counts': counts,
        'failures': 0, 'arithmetic': 'Python integers and fractions.Fraction only',
        'examples': examples,
        'sharpness_witness': {'lambda': '5',
                              'coefficient_t1_squared': str(lambda_five[degree_two_axis]),
                              'method': 'Direct rational expansion of (1-Q)^5 R'},
        'source_sha256': {str(p.relative_to(HERE.parents[3])): sha256(p.read_bytes()).hexdigest()
                          for p in (Path(__file__), PREVIOUS, HERE / 'sharp-factorization.md')},
        'scope': 'Finite tests of new recurrence, rational identity, magnitude bounds and boundary cases. The interval formula is inherited from cycle24, not the literal source-integral checker. All-orders validity requires analytic review; no full Lean proof or originality claim.'
    }
    args.output.write_text(json.dumps(report, indent=2) + '\n')
    print(json.dumps(report, indent=2))


if __name__ == '__main__':
    main()