#!/usr/bin/env python3
"""Independent verification of the Bui-Hall sign-conjecture proof.

Unlike the in-repo check (which starts from the reduced interval integral,
equation (3)), this starts from the LITERAL source integral of Bui-Hall
Theorem 3 as quoted in the proof:

    I_alpha = int_{[0,1]^4} (u1-u2)^2 A^k B^l C^m D^n du
    A = 1/2 - u1 + (u1-u2) u3     B = 1/2 - u2 - (u1-u2) u3
    C = 1/2 - u1 + (u1-u2) u4     D = 1/2 - u2 - (u1-u2) u4
    HARDY = 3 (-1)^(m+n) i^S I_alpha

and checks the proof's closed form

    HARDY = 12 (-1)^(S/2) / (2^S (S+4)!) * E[X1^k X2^l X3^m X4^n]

where X_i are the four linear forms in three iid symmetric Laplace variables.
Everything is exact rational arithmetic; no third-party packages.

This tests section 1 of the proof (the reduction), which the existing
repository check does not cover.
"""
from collections import defaultdict
from fractions import Fraction as F
from itertools import product
from math import factorial

# ---------------------------------------------------------------- polynomials
# sparse dict: exponent tuple -> Fraction coefficient


def pmul(p, q):
    out = defaultdict(F)
    for ea, ca in p.items():
        for eb, cb in q.items():
            out[tuple(x + y for x, y in zip(ea, eb))] += ca * cb
    return {e: c for e, c in out.items() if c}


def padd(p, q):
    out = defaultdict(F, p)
    for e, c in q.items():
        out[e] += c
    return {e: c for e, c in out.items() if c}


def ppow(p, n, nvars):
    r = {(0,) * nvars: F(1)}
    for _ in range(n):
        r = pmul(r, p)
    return r


# ------------------------------------------------- side A: the source integral
# variables (u1, u2, u3, u4)

def mono4(c, e1, e2, e3, e4):
    return {(e1, e2, e3, e4): F(c)}


def source_forms():
    half = mono4(1, 0, 0, 0, 0)
    half = {k: F(1, 2) for k in half}
    u1 = mono4(1, 1, 0, 0, 0)
    u2 = mono4(1, 0, 1, 0, 0)
    u1u3 = mono4(1, 1, 0, 1, 0)
    u2u3 = mono4(-1, 0, 1, 1, 0)
    u1u4 = mono4(1, 1, 0, 0, 1)
    u2u4 = mono4(-1, 0, 1, 0, 1)
    neg = lambda p: {e: -c for e, c in p.items()}
    A = padd(padd(half, neg(u1)), padd(u1u3, u2u3))
    B = padd(padd(half, neg(u2)), neg(padd(u1u3, u2u3)))
    C = padd(padd(half, neg(u1)), padd(u1u4, u2u4))
    D = padd(padd(half, neg(u2)), neg(padd(u1u4, u2u4)))
    return A, B, C, D


def integrate_unit_cube4(p):
    """int_{[0,1]^4} p du  -- monomialwise, exact."""
    tot = F(0)
    for (a, b, c, d), coef in p.items():
        tot += coef / F((a + 1) * (b + 1) * (c + 1) * (d + 1))
    return tot


def I_alpha_from_source(k, l, m, n):
    A, B, C, D = source_forms()
    u1mu2 = padd(mono4(1, 1, 0, 0, 0), mono4(-1, 0, 1, 0, 0))
    p = pmul(u1mu2, u1mu2)
    for form, e in ((A, k), (B, l), (C, m), (D, n)):
        p = pmul(p, ppow(form, e, 4))
    return integrate_unit_cube4(p)


def hardy_from_source(k, l, m, n):
    """HARDY = 3 (-1)^(m+n) i^S I.  Even S only, so i^S = (-1)^(S/2)."""
    S = k + l + m + n
    assert S % 2 == 0
    I = I_alpha_from_source(k, l, m, n)
    return F(3) * F((-1) ** (m + n)) * F((-1) ** (S // 2)) * I


# ------------------------------------- side B: the proof's Laplace expectation
# X1 = U+Z+H, X2 = U-Z-H, X3 = -U+Z-H, X4 = -U-Z+H, variables (U,Z,H)

LINFORMS = ((1, 1, 1), (1, -1, -1), (-1, 1, -1), (-1, -1, 1))


def laplace_expect_monomial(e):
    """E[U^a Z^b H^c] for iid symmetric Laplace: a! b! c! if all even else 0."""
    a, b, c = e
    if a % 2 or b % 2 or c % 2:
        return F(0)
    return F(factorial(a) * factorial(b) * factorial(c))


def moment_from_proof(k, l, m, n):
    p = {(0, 0, 0): F(1)}
    for (cu, cz, ch), e in zip(LINFORMS, (k, l, m, n)):
        form = {}
        if cu:
            form[(1, 0, 0)] = F(cu)
        if cz:
            form[(0, 1, 0)] = F(cz)
        if ch:
            form[(0, 0, 1)] = F(ch)
        p = pmul(p, ppow(form, e, 3))
    return sum((c * laplace_expect_monomial(e) for e, c in p.items()), F(0))


def hardy_from_proof(k, l, m, n):
    S = k + l + m + n
    assert S % 2 == 0
    E = moment_from_proof(k, l, m, n)
    return F(12) * F((-1) ** (S // 2)) * E / (F(2 ** S) * F(factorial(S + 4)))


# --------------------------------------------------------------- the sign law

def predicted_sign(k, l, m, n):
    """Proof's claim: (-1)^(sum floor(a_j/2)) * HARDY > 0."""
    return (-1) ** (k // 2 + l // 2 + m // 2 + n // 2)


def source_says_negative(k, l, m, n):
    """Conjecture 1 condition: negative exactly when |i^k+i^l+i^m+i^n| = 2.

    Exact integer test. i^e cycles through (1,0), (0,1), (-1,0), (0,-1) as
    (real, imag), so the sum is a Gaussian integer and |sum| = 2 is the
    integer condition re^2 + im^2 == 4. No floating point, no tolerance.
    """
    table = ((1, 0), (0, 1), (-1, 0), (0, -1))
    re = im = 0
    for e in (k, l, m, n):
        a, b = table[e % 4]
        re += a
        im += b
    return re * re + im * im == 4


# ------------------------------------------------------------------- the runs

def quadruples(maxdeg):
    for S in range(0, maxdeg + 1, 2):
        for k in range(S + 1):
            for l in range(S - k + 1):
                for m in range(S - k - l + 1):
                    yield (k, l, m, S - k - l - m)


def main():
    import sys; MAXDEG = int(sys.argv[1]) if len(sys.argv)>1 else 8
    n_red = n_sign = n_equiv = 0
    red_fail, sign_fail, equiv_fail = [], [], []

    for a in quadruples(MAXDEG):
        hs = hardy_from_source(*a)
        hp = hardy_from_proof(*a)
        n_red += 1
        if hs != hp:
            red_fail.append((a, hs, hp))

        n_sign += 1
        if not (predicted_sign(*a) * hs > 0):
            sign_fail.append((a, hs))

        n_equiv += 1
        if (hs < 0) != source_says_negative(*a):
            equiv_fail.append((a, hs))

    print("INDEPENDENT VERIFICATION -- Bui-Hall sign conjecture proof")
    print("Start point: the LITERAL source integral over [0,1]^4 (Thm 3).")
    print("=" * 64)
    print(f"quadruples tested (even total degree <= {MAXDEG}): {n_red}")
    print()
    print("(1) REDUCTION  source integral  ==  proof's Laplace closed form")
    print(f"    checked {n_red}, mismatches {len(red_fail)}")
    for f in red_fail[:5]:
        print("      FAIL", f)
    print()
    print("(2) SIGN LAW   (-1)^(sum floor(a_j/2)) * HARDY > 0")
    print(f"    checked {n_sign}, violations {len(sign_fail)}")
    for f in sign_fail[:5]:
        print("      FAIL", f)
    print()
    print("(3) EQUIV      HARDY < 0  <==>  |i^k+i^l+i^m+i^n| = 2")
    print(f"    checked {n_equiv}, mismatches {len(equiv_fail)}")
    for f in equiv_fail[:5]:
        print("      FAIL", f)
    print()
    print("normalization check  HARDY(0,0,0,0) =", hardy_from_source(0, 0, 0, 0),
          " (proof states 1/2)")
    print("odd-degree vanishing E[X^(1,0,0,0)] =", moment_from_proof(1, 0, 0, 0))
    print()
    sample = [(0, 0, 0, 0), (2, 0, 0, 0), (1, 1, 0, 0), (2, 2, 0, 0),
              (1, 1, 1, 1), (4, 0, 0, 0), (2, 1, 1, 0), (3, 1, 0, 0)]
    print("sample exact values (from the SOURCE integral):")
    for a in sample:
        if sum(a) % 2 == 0 and sum(a) <= MAXDEG:
            print(f"   HARDY{a} = {hardy_from_source(*a)}")

    failures = len(red_fail) + len(sign_fail) + len(equiv_fail)
    print()
    print("RESULT:", "ALL CHECKS PASSED" if failures == 0
          else f"{failures} FAILURE(S)")
    return 1 if failures else 0


if __name__ == "__main__":
    raise SystemExit(main())