#!/usr/bin/env python3
"""Exact certificate for the sharp degree-seven signed Karlin H6 threshold."""

from __future__ import annotations

import hashlib
from collections import defaultdict
from fractions import Fraction
from itertools import combinations, permutations
from math import comb, factorial

import sympy as sp


RECIPROCAL_ORBITS = {
    (6,): 125,
    (4, 2): -45,
    (4, 1, 1): 60,
    (3, 3): -56,
    (3, 2, 1): 72,
    (3, 1, 1, 1): -66,
    (2, 2, 2): 99,
    (2, 2, 1, 1): -90,
    (2, 1, 1, 1, 1): 324,
    (1, 1, 1, 1, 1, 1): -2340,
}
VARIABLES = 6
TOTAL_DEGREE = 6
SIX_ROOT_CONTENT = 23_887_872_000
FIVE_ROOT_CONTENT = 2_985_984_000_000
EXPECTED_SOURCE_TERMS = 432
EXPECTED_CHAMBER_TERMS = 462
EXPECTED_MINIMUM_COEFFICIENT = 104
EXPECTED_MAXIMUM_COEFFICIENT = 426_060
EXPECTED_COEFFICIENT_SUM = 26_978_400
EXPECTED_CONSTANT_RAY_COEFFICIENT = 8_100
EXPECTED_CHAMBER_SHA256 = (
    "30aed3f2fc017e69d048180129a3af3455b96f1dc41f73107286469310cff85c"
)
EXPECTED_BOUNDARY_WITNESS = -1_246_946_918_400
EXPECTED_WITNESS = (
    -4_202_074_601_182_498_649_700_878_160_382_679_442_880_356_812_390_400
)


if not __debug__:
    raise SystemExit("verify.py requires assertions")


def orbit_dictionary() -> dict[tuple[int, ...], int]:
    result: dict[tuple[int, ...], int] = defaultdict(int)
    for partition, coefficient in RECIPROCAL_ORBITS.items():
        padded = partition + (0,) * (VARIABLES - len(partition))
        for exponents in set(permutations(padded)):
            result[exponents] += coefficient
    return dict(result)


def dictionary_expression(
    coefficients: dict[tuple[int, ...], int],
    variables: tuple[sp.Symbol, ...],
) -> sp.Expr:
    return sp.Add(
        *(
            coefficient
            * sp.prod(
                variable**exponent
                for variable, exponent in zip(
                    variables, exponents, strict=True
                )
            )
            for exponents, coefficient in coefficients.items()
        )
    )


def reciprocal_determinant(
    support: int,
) -> tuple[tuple[sp.Symbol, ...], sp.Expr]:
    elementary = sp.symbols(f"e1:{support + 1}")
    values = (sp.Integer(1),) + elementary
    matrix = sp.Matrix(
        [
            [
                factorial(row + column)
                * (
                    values[support - row - column]
                    if 0 <= support - row - column <= support
                    else 0
                )
                for column in range(6)
            ]
            for row in range(6)
        ]
    )
    return elementary, sp.expand(-matrix.det())


def verify_orbit_identity(
    source: dict[tuple[int, ...], int],
) -> None:
    beta = sp.symbols("beta0:6")
    expression = dictionary_expression(source, beta)
    symmetric, remainder, mapping = sp.symmetrize(
        expression, *beta, formal=True
    )
    assert remainder == 0

    elementary, determinant = reciprocal_determinant(6)
    determinant_poly = sp.Poly(determinant, *elementary)
    assert sp.gcd_list(determinant_poly.coeffs()) == SIX_ROOT_CONTENT
    normalized = sp.expand(determinant / SIX_ROOT_CONTENT)
    formal_to_elementary = {
        formal: elementary[index]
        for index, (formal, _) in enumerate(mapping)
    }
    assert sp.expand(
        symmetric.subs(formal_to_elementary) - normalized
    ) == 0


def verify_five_root_formula() -> None:
    elementary, determinant = reciprocal_determinant(5)
    assert sp.Poly(determinant, *elementary).total_degree() == 0
    assert determinant == FIVE_ROOT_CONTENT


def ordered_difference_transform(
    coefficients: dict[tuple[int, ...], int],
) -> dict[tuple[int, ...], int]:
    transformed = coefficients
    for index in range(VARIABLES - 1):
        updated: dict[tuple[int, ...], int] = defaultdict(int)
        for exponents, coefficient in transformed.items():
            power = exponents[index]
            carried = exponents[index + 1]
            for retained in range(power + 1):
                new_exponents = list(exponents)
                new_exponents[index] = retained
                new_exponents[index + 1] = carried + power - retained
                updated[tuple(new_exponents)] += (
                    coefficient * comb(power, retained)
                )
        transformed = {
            exponents: coefficient
            for exponents, coefficient in updated.items()
            if coefficient
        }
    return transformed


def dictionary_digest(coefficients: dict[tuple[int, ...], int]) -> str:
    payload = "".join(
        ",".join(map(str, exponents)) + ":" + str(coefficient) + "\n"
        for exponents, coefficient in sorted(coefficients.items())
    ).encode()
    return hashlib.sha256(payload).hexdigest()


def evaluate_dictionary(
    coefficients: dict[tuple[int, ...], int],
    values: tuple[int | Fraction, ...],
) -> int | Fraction:
    return sum(
        coefficient
        * sp.prod(
            value**exponent
            for value, exponent in zip(values, exponents, strict=True)
        )
        for exponents, coefficient in coefficients.items()
    )


def verify_transform(
    source: dict[tuple[int, ...], int],
    chamber: dict[tuple[int, ...], int],
) -> None:
    assert len(source) == EXPECTED_SOURCE_TERMS
    assert len(chamber) == EXPECTED_CHAMBER_TERMS
    assert EXPECTED_CHAMBER_TERMS == comb(
        TOTAL_DEGREE + VARIABLES - 1, VARIABLES - 1
    )
    assert all(
        sum(exponents) == TOTAL_DEGREE
        for exponents in chamber
    )
    assert min(chamber.values()) == EXPECTED_MINIMUM_COEFFICIENT
    assert max(chamber.values()) == EXPECTED_MAXIMUM_COEFFICIENT
    assert sum(chamber.values()) == EXPECTED_COEFFICIENT_SUM
    assert all(coefficient > 0 for coefficient in chamber.values())
    assert (
        chamber[(0, 0, 0, 0, 0, 6)]
        == EXPECTED_CONSTANT_RAY_COEFFICIENT
    )
    assert dictionary_digest(chamber) == EXPECTED_CHAMBER_SHA256

    for gaps in (
        (1, 1, 1, 1, 1, 1),
        (1, 2, 3, 4, 5, 6),
        (6, 1, 4, 2, 5, 3),
    ):
        beta = tuple(sum(gaps[index:]) for index in range(VARIABLES))
        assert evaluate_dictionary(source, beta) == evaluate_dictionary(
            chamber, gaps
        )


def elementary_symmetric(
    values: tuple[int | Fraction | sp.Expr, ...],
    degree: int,
) -> int | Fraction | sp.Expr:
    return sum(
        sp.prod(values[index] for index in subset)
        for subset in combinations(range(len(values)), degree)
    )


def centered_determinant(
    values: tuple[int | Fraction | sp.Expr, ...],
) -> sp.Expr:
    elementary = [
        elementary_symmetric(values, degree)
        for degree in range(11)
    ]
    derivatives = [
        factorial(degree) * elementary[degree]
        for degree in range(11)
    ]
    return sp.expand(
        sp.Matrix(
            [
                [
                    derivatives[5 + column - row]
                    for column in range(6)
                ]
                for row in range(6)
            ]
        ).det()
    )


def verify_reciprocal_relation(
    source: dict[tuple[int, ...], int],
) -> None:
    alpha = (1, 2, 3, 4, 5, 6)
    beta = tuple(Fraction(1, value) for value in alpha)
    product_alpha = sp.prod(alpha)
    reciprocal_value = evaluate_dictionary(source, beta)
    assert centered_determinant(alpha) == (
        SIX_ROOT_CONTENT
        * product_alpha**6
        * reciprocal_value
    )


def verify_polynomial_transfer() -> None:
    variable = sp.symbols("x")
    roots = (0, 1, 2, 3, 4, 5)
    polynomial = -sp.prod(variable + root for root in roots)
    point = 2
    derivatives = [
        sp.diff(polynomial, variable, degree).subs(variable, point)
        for degree in range(11)
    ]
    hankel = sp.Matrix(
        [
            [
                derivatives[row + column]
                for column in range(6)
            ]
            for row in range(6)
        ]
    ).det()
    reciprocal_roots = tuple(
        Fraction(1, point + root) for root in roots
    )
    assert -hankel == (
        polynomial.subs(variable, point) ** 6
        * centered_determinant(reciprocal_roots)
    )
    assert -hankel > 0


def verify_witness_family() -> None:
    scale = sp.symbols("A")
    determinant = sp.factor(
        centered_determinant((scale, scale, scale, scale, 1, 1, 1))
    )
    polynomial = (
        116 * scale**12
        - 5952 * scale**11
        - 5544 * scale**10
        - 14880 * scale**9
        - 36540 * scale**8
        - 44352 * scale**7
        - 76436 * scale**6
        - 73632 * scale**5
        - 71595 * scale**4
        - 28560 * scale**3
        - 47664 * scale**2
        - 10512 * scale
        - 4624
    )
    expected = -10_749_542_400 * scale**12 * polynomial
    assert sp.expand(determinant - expected) == 0
    assert determinant.subs(scale, 52) > 0
    assert determinant.subs(scale, 53) == EXPECTED_WITNESS


def verify_simple_root_witness() -> None:
    denominator = 10_000_000
    numerators = (
        53 * denominator - 3,
        53 * denominator - 1,
        53 * denominator + 1,
        53 * denominator + 3,
        denominator - 2,
        denominator,
        denominator + 2,
    )
    reciprocal_roots = tuple(
        Fraction(numerator, denominator)
        for numerator in numerators
    )
    assert len(set(reciprocal_roots)) == 7
    assert min(reciprocal_roots) > 0
    assert centered_determinant(reciprocal_roots) < 0


def verify_boundary_witness() -> None:
    # f(x)=x^4(x+1)^3, represented by ascending coefficients.
    coefficients = (0, 0, 0, 0) + tuple(
        comb(3, degree) for degree in range(4)
    )
    derivatives = [
        factorial(degree) * (
            coefficients[degree] if degree < len(coefficients) else 0
        )
        for degree in range(11)
    ]
    hankel = sp.Matrix(
        [
            [
                derivatives[row + column]
                for column in range(6)
            ]
            for row in range(6)
        ]
    ).det()
    assert -hankel == EXPECTED_BOUNDARY_WITNESS


def main() -> None:
    source = orbit_dictionary()
    verify_orbit_identity(source)
    verify_five_root_formula()
    chamber = ordered_difference_transform(source)
    verify_transform(source, chamber)
    verify_reciprocal_relation(source)
    verify_polynomial_transfer()
    verify_boundary_witness()
    verify_witness_family()
    verify_simple_root_witness()
    assert (-1) ** (6 * 5 // 2) == -1
    print("six-root reciprocal determinant and 10 orbit coefficients: PASS")
    print("five-root closed formula: PASS")
    print("ordered chamber source terms:", len(source))
    print("ordered chamber expanded terms:", len(chamber))
    print("ordered chamber minimum coefficient:", min(chamber.values()))
    print("ordered chamber SHA-256:", dictionary_digest(chamber))
    print("reciprocal complement identity: PASS")
    print("negative-leading polynomial transfer: PASS")
    print("degree-seven boundary witness: PASS")
    print("degree-seven witness family: PASS")
    print("degree-seven simple-root rational witness: PASS")


if __name__ == "__main__":
    main()
