#!/usr/bin/env python3
"""Dependency-free hostile verifier for C-0059 and its primary sources."""

from __future__ import annotations

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


SIZE = 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_FAMILY_WITNESS = (
    -4_202_074_601_182_498_649_700_878_160_382_679_442_880_356_812_390_400
)
EXPECTED_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,
}
EXPECTED_SOURCE_SHA256 = {
    "cc2002-page16.png":
        "6c7b50158ff09e43285c05ad6e3ecead08c8b708d2bd624d2e4d33884713438f",
    "cc2002-page17.png":
        "bf1e3c33ab0130949222ab1ad6b445d9c07f4ffa1a864926c68c5977dcc4a993",
    "cc2005-page64.png":
        "96eb61475e398a8150f0cd681f4d36b888c4485be99aabbfa6308c6e907afbec",
    "cc2005-page67.png":
        "6413d21d1938a751cacd93b3552fde39270099c5c501a5ec790acc09d7b9d4b7",
    "cc2005-page68.png":
        "de8e54c4fc8b7451906b7a7b90c9173bd8cc117a5cce6bb9009be8cb2689ccb9",
    "craven-csordas-2002-iterated-laguerre-turan.pdf":
        "945310b929c47aee17cbc5f0a8bf3fb25cd6c52f14da88c5117dd9bed3b59714",
    "craven-csordas-2002-iterated-laguerre-turan.txt":
        "db7f2cda0b0e3ca995c42f7854dbeaa4c4940368dce2e357821c3eeb4e04a318",
    "craven-csordas-2005-karlin-conjecture.pdf":
        "bd9fc2f503621a371c2cbda72a7130df5ba2ba56b776260d865ed96f2e6c6746",
    "craven-csordas-2005-karlin-conjecture.txt":
        "49603f4dd6fe6e5c3da0fe84004b3512476d9b3d244258a3f1934b946499dac9",
    "original-verify.py":
        "5406185dff46e2aae06a5a55ecd3d8afb4dd4a4c9a204cbc74bc4f44f29115d3",
}

Exponent = tuple[int, ...]
Polynomial = dict[Exponent, int]
UniPolynomial = tuple[int, ...]


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


def add(left: Polynomial, right: Polynomial) -> Polynomial:
    result: defaultdict[Exponent, int] = defaultdict(int)
    for exponents, coefficient in left.items():
        result[exponents] += coefficient
    for exponents, coefficient in right.items():
        result[exponents] += coefficient
    return {key: value for key, value in result.items() if value}


def scale(polynomial: Polynomial, scalar: int) -> Polynomial:
    return {
        exponents: scalar * coefficient
        for exponents, coefficient in polynomial.items()
        if scalar * coefficient
    }


def multiply(left: Polynomial, right: Polynomial) -> Polynomial:
    result: defaultdict[Exponent, int] = defaultdict(int)
    for left_exp, left_coefficient in left.items():
        for right_exp, right_coefficient in right.items():
            exponents = tuple(
                a + b
                for a, b in zip(left_exp, right_exp, strict=True)
            )
            result[exponents] += left_coefficient * right_coefficient
    return {key: value for key, value in result.items() if value}


def permutation_sign(permutation: tuple[int, ...]) -> int:
    inversions = sum(
        permutation[i] > permutation[j]
        for i in range(len(permutation))
        for j in range(i + 1, len(permutation))
    )
    return -1 if inversions % 2 else 1


def determinant(matrix: list[list[Polynomial]], variables: int) -> Polynomial:
    result: Polynomial = {}
    one = {(0,) * variables: 1}
    for permutation in permutations(range(len(matrix))):
        term = one
        for row, column in enumerate(permutation):
            term = multiply(term, matrix[row][column])
            if not term:
                break
        result = add(result, scale(term, permutation_sign(permutation)))
    return result


def elementary_polynomials(variables: int) -> list[Polynomial]:
    result: list[Polynomial] = []
    for degree in range(variables + 1):
        polynomial: Polynomial = {}
        for subset in combinations(range(variables), degree):
            exponents = [0] * variables
            for index in subset:
                exponents[index] = 1
            polynomial[tuple(exponents)] = 1
        result.append(polynomial)
    return result


def reciprocal_determinant(support: int) -> Polynomial:
    elementary = elementary_polynomials(support)
    zero: Polynomial = {}
    matrix = [
        [
            scale(
                elementary[support - row - column]
                if 0 <= support - row - column <= support
                else zero,
                factorial(row + column),
            )
            for column in range(SIZE)
        ]
        for row in range(SIZE)
    ]
    return scale(determinant(matrix, support), -1)


def derive_six_root_source() -> Polynomial:
    raw = reciprocal_determinant(6)
    content = 0
    for coefficient in raw.values():
        content = gcd(content, abs(coefficient))
    assert content == SIX_ROOT_CONTENT
    assert all(coefficient % content == 0 for coefficient in raw.values())
    source = {
        exponents: coefficient // content
        for exponents, coefficient in raw.items()
    }
    assert len(source) == EXPECTED_SOURCE_TERMS
    assert all(sum(exponents) == 6 for exponents in source)
    return source


def orbit_coefficients(source: Polynomial) -> dict[Exponent, int]:
    coefficients: dict[Exponent, int] = {}
    multiplicities: defaultdict[Exponent, int] = defaultdict(int)
    for exponents, coefficient in source.items():
        partition = tuple(
            sorted((power for power in exponents if power), reverse=True)
        )
        multiplicities[partition] += 1
        previous = coefficients.setdefault(partition, coefficient)
        assert previous == coefficient
    for partition, multiplicity in multiplicities.items():
        padded = partition + (0,) * (SIZE - len(partition))
        assert multiplicity == len(set(permutations(padded)))
    return coefficients


def ordered_difference_transform(source: Polynomial) -> Polynomial:
    transformed = source
    for index in range(SIZE - 1):
        updated: defaultdict[Exponent, int] = defaultdict(int)
        for exponents, coefficient in transformed.items():
            power = exponents[index]
            for retained in range(power + 1):
                next_exponents = list(exponents)
                next_exponents[index] = retained
                next_exponents[index + 1] += power - retained
                updated[tuple(next_exponents)] += (
                    coefficient * comb(power, retained)
                )
        transformed = {key: value for key, value in updated.items() if value}
    return transformed


def dictionary_digest(polynomial: Polynomial) -> str:
    payload = "".join(
        ",".join(map(str, exponents)) + ":" + str(coefficient) + "\n"
        for exponents, coefficient in sorted(polynomial.items())
    ).encode()
    return hashlib.sha256(payload).hexdigest()


def evaluate_dictionary(
    polynomial: Polynomial,
    values: tuple[int | Fraction, ...],
) -> int | Fraction:
    return sum(
        coefficient
        * product(
            value**power
            for value, power in zip(values, exponents, strict=True)
        )
        for exponents, coefficient in polynomial.items()
    )


def verify_chamber(source: Polynomial) -> Polynomial:
    chamber = ordered_difference_transform(source)
    assert len(chamber) == EXPECTED_CHAMBER_TERMS
    assert EXPECTED_CHAMBER_TERMS == comb(11, 5)
    assert all(sum(exponents) == 6 for exponents in chamber)
    assert all(coefficient > 0 for coefficient in chamber.values())
    assert min(chamber.values()) == EXPECTED_MINIMUM_COEFFICIENT
    assert max(chamber.values()) == EXPECTED_MAXIMUM_COEFFICIENT
    assert sum(chamber.values()) == EXPECTED_COEFFICIENT_SUM
    assert (
        chamber[(0, 0, 0, 0, 0, 6)]
        == EXPECTED_CONSTANT_RAY_COEFFICIENT
    )
    assert dictionary_digest(chamber) == EXPECTED_CHAMBER_SHA256
    for gaps in ((1, 2, 3, 4, 5, 6), (6, 1, 4, 2, 5, 3)):
        beta = tuple(sum(gaps[index:]) for index in range(SIZE))
        assert evaluate_dictionary(source, beta) == evaluate_dictionary(
            chamber, gaps
        )
    return chamber


def product(values) -> int | Fraction:
    result = 1
    for value in values:
        result *= value
    return result


def elementary_values(values: tuple[int | Fraction, ...]) -> list[int | Fraction]:
    result: list[int | Fraction] = [1] + [0] * len(values)
    for value in values:
        for degree in range(len(values), 0, -1):
            result[degree] += value * result[degree - 1]
    return result


def scalar_determinant(matrix: list[list[int | Fraction]]) -> int | Fraction:
    total: int | Fraction = 0
    for permutation in permutations(range(len(matrix))):
        total += permutation_sign(permutation) * product(
            matrix[row][column]
            for row, column in enumerate(permutation)
        )
    return total


def centered_determinant(values: tuple[int | Fraction, ...]) -> int | Fraction:
    elementary = elementary_values(values)
    derivatives = [
        factorial(degree) * (
            elementary[degree] if degree < len(elementary) else 0
        )
        for degree in range(11)
    ]
    return scalar_determinant(
        [
            [derivatives[5 + column - row] for column in range(SIZE)]
            for row in range(SIZE)
        ]
    )


def verify_lower_support_and_reciprocal(source: Polynomial) -> None:
    five = reciprocal_determinant(5)
    assert five == {(0, 0, 0, 0, 0): FIVE_ROOT_CONTENT}
    for support in range(5):
        assert reciprocal_determinant(support) == {}

    alpha = (1, 2, 3, 4, 5, 6)
    beta = tuple(Fraction(1, value) for value in alpha)
    assert centered_determinant(alpha) == (
        SIX_ROOT_CONTENT
        * product(alpha) ** 6
        * evaluate_dictionary(source, beta)
    )


def polynomial_multiply(left: list[int], right: list[int]) -> list[int]:
    result = [0] * (len(left) + len(right) - 1)
    for left_degree, left_coefficient in enumerate(left):
        for right_degree, right_coefficient in enumerate(right):
            result[left_degree + right_degree] += (
                left_coefficient * right_coefficient
            )
    while len(result) > 1 and result[-1] == 0:
        result.pop()
    return result


def derivative_value(coefficients: list[int], degree: int, point: int) -> int:
    return sum(
        coefficient
        * factorial(power) // factorial(power - degree)
        * point ** (power - degree)
        for power, coefficient in enumerate(coefficients)
        if power >= degree
    )


def verify_transfer_and_boundary() -> None:
    coefficients = [-1]
    for root in (0, 1, 2, 3, 4, 5):
        coefficients = polynomial_multiply(coefficients, [root, 1])
    point = 2
    derivatives = [
        derivative_value(coefficients, degree, point)
        for degree in range(11)
    ]
    hankel = scalar_determinant(
        [
            [derivatives[row + column] for column in range(SIZE)]
            for row in range(SIZE)
        ]
    )
    reciprocal_roots = tuple(
        Fraction(1, point + root) for root in (0, 1, 2, 3, 4, 5)
    )
    value = sum(
        coefficient * point**power
        for power, coefficient in enumerate(coefficients)
    )
    assert -hankel == value**6 * centered_determinant(reciprocal_roots)
    assert -hankel > 0

    boundary_coefficients = [0, 0, 0, 0] + [
        comb(3, degree) for degree in range(4)
    ]
    boundary_derivatives = [
        factorial(degree) * (
            boundary_coefficients[degree]
            if degree < len(boundary_coefficients)
            else 0
        )
        for degree in range(11)
    ]
    boundary_hankel = scalar_determinant(
        [
            [
                boundary_derivatives[row + column]
                for column in range(SIZE)
            ]
            for row in range(SIZE)
        ]
    )
    assert -boundary_hankel == EXPECTED_BOUNDARY_WITNESS


def uni_add(left: UniPolynomial, right: UniPolynomial) -> UniPolynomial:
    size = max(len(left), len(right))
    result = tuple(
        (left[index] if index < len(left) else 0)
        + (right[index] if index < len(right) else 0)
        for index in range(size)
    )
    return uni_trim(result)


def uni_scale(polynomial: UniPolynomial, scalar: int) -> UniPolynomial:
    return uni_trim(tuple(scalar * coefficient for coefficient in polynomial))


def uni_multiply(
    left: UniPolynomial,
    right: UniPolynomial,
) -> UniPolynomial:
    result = [0] * (len(left) + len(right) - 1)
    for left_degree, left_coefficient in enumerate(left):
        for right_degree, right_coefficient in enumerate(right):
            result[left_degree + right_degree] += (
                left_coefficient * right_coefficient
            )
    return uni_trim(tuple(result))


def uni_trim(polynomial: UniPolynomial) -> UniPolynomial:
    end = len(polynomial)
    while end > 1 and polynomial[end - 1] == 0:
        end -= 1
    return polynomial[:end]


def uni_determinant(matrix: list[list[UniPolynomial]]) -> UniPolynomial:
    result: UniPolynomial = (0,)
    for permutation in permutations(range(len(matrix))):
        term: UniPolynomial = (1,)
        for row, column in enumerate(permutation):
            term = uni_multiply(term, matrix[row][column])
        result = uni_add(
            result,
            uni_scale(term, permutation_sign(permutation)),
        )
    return result


def uni_evaluate(polynomial: UniPolynomial, value: int) -> int:
    result = 0
    for coefficient in reversed(polynomial):
        result = result * value + coefficient
    return result


def verify_witness_family() -> None:
    elementary: list[UniPolynomial] = []
    for degree in range(8):
        coefficients = [0] * 5
        for powers_of_scale in range(5):
            powers_of_one = degree - powers_of_scale
            if 0 <= powers_of_one <= 3:
                coefficients[powers_of_scale] += (
                    comb(4, powers_of_scale) * comb(3, powers_of_one)
                )
        elementary.append(uni_trim(tuple(coefficients)))
    derivatives = [
        uni_scale(elementary[degree], factorial(degree))
        if degree < len(elementary)
        else (0,)
        for degree in range(11)
    ]
    determinant = uni_determinant(
        [
            [derivatives[5 + column - row] for column in range(SIZE)]
            for row in range(SIZE)
        ]
    )
    factor = -10_749_542_400
    factor_polynomial = (
        -4624, -10512, -47664, -28560, -71595, -73632, -76436,
        -44352, -36540, -14880, -5544, -5952, 116,
    )
    expected = (0,) * 12 + tuple(
        factor * coefficient for coefficient in factor_polynomial
    )
    assert determinant == expected
    assert uni_evaluate(determinant, 52) > 0
    assert uni_evaluate(determinant, 53) == EXPECTED_FAMILY_WITNESS

    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 sha256(path: Path) -> str:
    return hashlib.sha256(path.read_bytes()).hexdigest()


def verify_sources() -> None:
    source_dir = Path(__file__).resolve().parent / "source"
    assert {path.name for path in source_dir.iterdir()} == set(
        EXPECTED_SOURCE_SHA256
    )
    for name, expected in EXPECTED_SOURCE_SHA256.items():
        assert sha256(source_dir / name) == expected

    text_2005 = (
        source_dir / "craven-csordas-2005-karlin-conjecture.txt"
    ).read_text(errors="strict")
    assert "Karlin’s Conjecture. Let φ ∈ L-P + . Then for all n ≥ 1" in text_2005
    assert "(1.7)   (−1)n(n+1)/2 W (f (x), f" in text_2005
    assert "Example 3.1. Let f (x) := x3 (x + 1)5 . Then the 5 × 5 Wronskian" in text_2005
    assert "x4 (x + 1)5 is not a counterexample in the 5 × 5 Wronskian" in text_2005
    assert (-1) ** (5 * 6 // 2) == -1

    text_2002 = (
        source_dir / "craven-csordas-2002-iterated-laguerre-turan.txt"
    ).read_text(errors="strict")
    assert "In connection with the investigations of a conjecture of" in text_2002
    assert "S. Karlin, additional examples are considered in [7] and [1]." in text_2002


def main() -> None:
    verify_sources()
    source = derive_six_root_source()
    assert orbit_coefficients(source) == EXPECTED_ORBITS
    chamber = verify_chamber(source)
    verify_lower_support_and_reciprocal(source)
    verify_transfer_and_boundary()
    verify_witness_family()
    print("primary-source sign and prior-example scope: PASS")
    print("six-root determinant content:", SIX_ROOT_CONTENT)
    print("six-root source monomials:", len(source))
    print("six-root orbit coefficients:", len(EXPECTED_ORBITS))
    print("ordered chamber monomials:", len(chamber))
    print("ordered chamber minimum coefficient:", min(chamber.values()))
    print("ordered chamber SHA-256:", dictionary_digest(chamber))
    print("five-root formula and support exhaustion: PASS")
    print("polynomial transfer and degree-seven boundary: PASS")
    print("degree-seven family and simple-root witness: PASS")


if __name__ == "__main__":
    main()
