#!/usr/bin/env python3
"""Exact finite checks, not a proof of the infinite Hardy sign theorem.

Compare independently expanded integral moments with a rational MGF; verify
parity, normalization and the sign dictionary. No third-party packages.
"""
from collections import defaultdict
from fractions import Fraction as F
from functools import lru_cache
from hashlib import sha256
from itertools import product
from math import comb, factorial
from pathlib import Path
import json

ZERO = (0, 0, 0, 0)
FORMS = ((1, 1, 1, 1), (1, 1, -1, -1),
         (1, -1, 1, -1), (1, -1, -1, 1))


def indices(total):
    for a in range(total + 1):
        for b in range(total - a + 1):
            for c in range(total - a - b + 1):
                yield (a, b, c, total - a - b - c)


def mul(p, q):
    out = defaultdict(int)
    for a, x in p.items():
        for b, y in q.items():
            out[tuple(u + v for u, v in zip(a, b))] += x * y
    return {a: x for a, x in out.items() if x}


def denominator(forms):
    out = {ZERO: 1}
    for signs in forms:
        linear = {tuple(int(i == j) for j in range(4)): signs[i]
                  for i in range(4)}
        factor = {a: -c for a, c in mul(linear, linear).items()}
        factor[ZERO] = 1
        out = mul(out, factor)
    return out


def inverse_series(poly, degree):
    assert poly[ZERO] == 1
    out = {ZERO: 1}
    nonconstant = [(a, c) for a, c in poly.items() if a != ZERO]
    for total in range(1, degree + 1):
        for a in indices(total):
            value = 0
            for b, c in nonconstant:
                if all(x >= y for x, y in zip(a, b)):
                    value -= c * out.get(tuple(x-y for x, y in zip(a, b)), 0)
            if value:
                out[a] = value
    return out


@lru_cache(None)
def interval_coeff(k, l):
    return {j: F(sum(comb(k, a) * comb(l, j-a) * (-1)**(j-a)
                         for a in range(max(0, j-l), min(k, j)+1)), j+1)
            for j in range(0, k+l+1, 2)}


def integral_hardy(alpha):
    # J_kl in center/radius coordinates, followed by exact simplex integrals.
    k, l, m, n = alpha
    s = sum(alpha)
    if s % 2:
        return F(0)
    value = sum(v*w*factorial(s-j-h)*factorial(j+h+2)
                for j, v in interval_coeff(k, l).items()
                for h, w in interval_coeff(m, n).items())
    return value * F(6 * (-1)**(m+n+s//2), 2**s * factorial(s+4))


def source_negative(alpha):
    units = ((1, 0), (0, 1), (-1, 0), (0, -1))
    re = sum(units[a % 4][0] for a in alpha)
    im = sum(units[a % 4][1] for a in alpha)
    return re*re + im*im == 4


def main():
    degree = 12
    mgf = inverse_series(denominator(FORMS[1:]), degree)
    full = inverse_series(denominator(FORMS), degree)
    comparisons = 0
    for total in range(degree + 1):
        for alpha in indices(total):
            c = mgf.get(alpha, 0)
            r = sum(a % 2 for a in alpha)
            assert full.get(alpha, 0) >= 0
            if r not in (0, 4):
                assert full.get(alpha, 0) == 0
            if total % 2:
                assert c == 0
                continue
            moment = c
            for a in alpha:
                moment *= factorial(a)
            value = F(12*(-1)**(total//2)*moment,
                      2**total*factorial(total+4))
            assert value == integral_hardy(alpha), (alpha, value)
            assert (-1)**(r//2) * c > 0, alpha
            expected = (-1)**sum(a//2 for a in alpha)
            assert expected * value > 0
            assert (expected < 0) == source_negative(alpha)
            comparisons += 1
    extended = 0
    for total in range(0, 41, 2):
        for alpha in indices(total):
            if tuple(sorted(alpha, reverse=True)) != alpha:
                continue
            value = integral_hardy(alpha)
            assert (-1)**sum(a//2 for a in alpha) * value > 0, alpha
            extended += 1
    samples = {(0,0,0,0): F(1,2), (6,0,0,0): F(-1,2688),
               (3,1,1,1): F(-19,201600), (3,3,1,1): F(19,1774080)}
    for alpha, expected in samples.items():
        assert integral_hardy(alpha) == expected
    print(json.dumps({
        'status': 'pass', 'scope': 'finite exact checks only',
        'ordered_even_quadruples_through_total_12': comparisons,
        'sorted_even_quadruples_through_total_40': extended,
        'sample_values': {str(a): str(v) for a, v in samples.items()},
        'script_sha256': sha256(Path(__file__).read_bytes()).hexdigest(),
        'proof_sha256': sha256(Path(__file__).with_name('hardy-sign-proof.md').read_bytes()).hexdigest()
    }, indent=2))


if __name__ == '__main__':
    main()