#!/usr/bin/env python3
"""Exact certificate for the sharp degree-eight Karlin H5 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 = {
    (5, 5, 5): 96,
    (5, 5, 3, 2): -30,
    (5, 5, 3, 1, 1): 60,
    (5, 5, 2, 2, 1): 70,
    (5, 5, 2, 1, 1, 1): -100,
    (5, 5, 1, 1, 1, 1, 1): 720,
    (5, 3, 3, 3, 1): -30,
    (5, 3, 3, 2, 2): 40,
    (5, 3, 3, 2, 1, 1): 50,
    (5, 3, 3, 1, 1, 1, 1): 120,
    (5, 3, 2, 2, 2, 1): -30,
    (5, 3, 2, 2, 1, 1, 1): 50,
    (5, 2, 2, 2, 2, 2): 540,
    (5, 2, 2, 2, 2, 1, 1): -120,
    (3, 3, 3, 3, 3): 40,
    (3, 3, 3, 3, 1, 1, 1): -45,
    (3, 3, 3, 2, 2, 2): 90,
    (3, 3, 3, 2, 2, 1, 1): 210,
    (3, 3, 2, 2, 2, 2, 1): -210,
    (3, 2, 2, 2, 2, 2, 2): 900,
}
VARIABLES = 7
DETERMINANT_CONTENT = 82_944
EXPECTED_SOURCE_TERMS = 6_496
EXPECTED_CHAMBER_TERMS = 48_503
EXPECTED_MINIMUM_COEFFICIENT = 24
EXPECTED_MAXIMUM_COEFFICIENT = 452_485_502_400
EXPECTED_COEFFICIENT_SUM = 370_753_878_476_800
EXPECTED_CONSTANT_RAY_COEFFICIENT = 180_075
EXPECTED_CHAMBER_SHA256 = (
    "d4ae7b6db47d4491756a677bf1e84a5d3b497f5548e2f4a34fb972790a94e343"
)
EXPECTED_BOUNDARY_WITNESS = -9_953_280
EXPECTED_DEGREE_EIGHT_WITNESS = -52_168_558_894_579_833_105_369_392_640_000


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() -> tuple[tuple[sp.Symbol, ...], sp.Expr]:
    elementary = sp.symbols("e1:8")
    values = (sp.Integer(1),) + elementary
    matrix = sp.Matrix(
        [
            [
                factorial(row + column)
                * (
                    values[7 - row - column]
                    if row + column <= 7
                    else 0
                )
                for column in range(5)
            ]
            for row in range(5)
        ]
    )
    determinant = sp.Poly(sp.expand(matrix.det()), *elementary)
    assert sp.gcd_list(determinant.coeffs()) == DETERMINANT_CONTENT
    return elementary, sp.expand(
        determinant.as_expr() / DETERMINANT_CONTENT
    )


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

    elementary, determinant = reciprocal_determinant()
    formal_to_elementary = {
        formal: elementary[index]
        for index, (formal, _) in enumerate(mapping)
    }
    assert sp.expand(
        symmetric.subs(formal_to_elementary) - determinant
    ) == 0


def ordered_difference_transform(
    coefficients: dict[tuple[int, ...], int],
) -> dict[tuple[int, ...], int]:
    transformed = coefficients
    # Substitute beta_i=y_i+beta_(i+1) successively. After six steps,
    # beta_i=y_i+...+y_7 on the ordered chamber.
    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 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, 0, 15)]
        == EXPECTED_CONSTANT_RAY_COEFFICIENT
    )
    assert dictionary_digest(chamber) == EXPECTED_CHAMBER_SHA256

    # Independent pointwise checks of the triangular substitution.
    for gaps in (
        (1, 1, 1, 1, 1, 1, 1),
        (1, 2, 3, 4, 5, 6, 7),
        (7, 1, 4, 2, 5, 3, 6),
    ):
        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, ...],
    degree: int,
) -> int | Fraction:
    return sum(
        sp.prod(values[index] for index in subset)
        for subset in combinations(range(len(values)), degree)
    )


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


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


def witness_determinant(scale: int) -> int:
    values = (scale, scale, scale, 1, 1, 1, 1, 1)
    return centered_determinant(values)


def boundary_witness_determinant() -> int:
    # f(x)=x^3(x+1)^5, represented by ascending coefficients.
    coefficients = (0, 0, 0) + tuple(comb(5, degree) for degree in range(6))
    derivatives = [
        factorial(degree) * coefficients[degree]
        for degree in range(9)
    ]
    return int(
        sp.Matrix(
            [
                [
                    derivatives[row + column]
                    for column in range(5)
                ]
                for row in range(5)
            ]
        ).det()
    )


def main() -> None:
    source = orbit_dictionary()
    verify_orbit_identity(source)
    chamber = ordered_difference_transform(source)
    verify_transform(source, chamber)
    verify_reciprocal_relation(source)
    assert boundary_witness_determinant() == EXPECTED_BOUNDARY_WITNESS
    assert witness_determinant(50) == EXPECTED_DEGREE_EIGHT_WITNESS
    assert (-1) ** (5 * 4 // 2) == 1
    print("reciprocal determinant and 20 orbit coefficients: 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("degree-eight boundary witness: PASS")
    print("degree-eight strict-root witness: PASS")


if __name__ == "__main__":
    main()
