#!/usr/bin/env python3
"""Dependency-free hostile verifier for C-0043 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
from pathlib import Path


EXPECTED_SOURCE_SHA256 = {
    "cc2002-page13.png":
        "3e2b8ded781cc5f8a6f9e6cdf44a168991b9c9fb23a7ab01580f3aecad637c84",
    "cc2002-page16.png":
        "6c7b50158ff09e43285c05ad6e3ecead08c8b708d2bd624d2e4d33884713438f",
    "cc2002-page17.png":
        "bf1e3c33ab0130949222ab1ad6b445d9c07f4ffa1a864926c68c5977dcc4a993",
    "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":
        "5049666042b70c35fec02bddf6a22deb1ea797239ae7f8bfffc43ffcce915e09",
}
EXPECTED_ORBITS = {
    (4, 4, 4): 108,
    (3, 3, 3, 3): -288,
    (4, 3, 3, 2): 96,
    (4, 4, 2, 2): -24,
    (3, 3, 2, 2, 2): 96,
    (4, 2, 2, 2, 2): 72,
    (2, 2, 2, 2, 2, 2): 720,
}
EXPECTED_Q_ORBITS = {
    (2, 2, 2): 9,
    (2, 2, 1, 1): -2,
    (2, 1, 1, 1, 1): 6,
    (1, 1, 1, 1, 1, 1): 60,
}
EXPECTED_STRICT_WITNESS = -527_889_869_252_540_784

Q_VARIABLES = 6
ROOT_VARIABLES = 6
ZERO_Q_EXPONENT = (0,) * Q_VARIABLES
ZERO_ROOT_EXPONENT = (0,) * ROOT_VARIABLES

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


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


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


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


def multiply_polynomials(left: Polynomial, right: Polynomial) -> Polynomial:
    result: defaultdict[tuple[int, ...], int] = defaultdict(int)
    for left_exponents, left_coefficient in left.items():
        for right_exponents, right_coefficient in right.items():
            exponents = tuple(
                left_power + right_power
                for left_power, right_power in zip(
                    left_exponents, right_exponents, strict=True
                )
            )
            result[exponents] += left_coefficient * right_coefficient
    return {
        exponents: coefficient
        for exponents, coefficient in result.items()
        if coefficient
    }


def power_polynomial(polynomial: Polynomial, exponent: int) -> Polynomial:
    zero = (0,) * len(next(iter(polynomial)))
    result = {zero: 1}
    base = polynomial
    while exponent:
        if exponent & 1:
            result = multiply_polynomials(result, base)
        base = multiply_polynomials(base, base)
        exponent //= 2
    return result


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


def determinant_polynomial(
    matrix: list[list[Polynomial]],
) -> Polynomial:
    size = len(matrix)
    result: Polynomial = {}
    for permutation in permutations(range(size)):
        term = {(0,) * len(next(iter(matrix[0][0]))): 1}
        for row, column in enumerate(permutation):
            term = multiply_polynomials(term, matrix[row][column])
        result = add_polynomials(
            result,
            scale_polynomial(term, permutation_sign(permutation)),
        )
    return result


def normalized_power_sum_determinant() -> Polynomial:
    q_variables = []
    for index in range(Q_VARIABLES):
        exponents = [0] * Q_VARIABLES
        exponents[index] = 1
        q_variables.append({tuple(exponents): 1})

    derivatives: list[Polynomial] = [{ZERO_Q_EXPONENT: 1}]
    for order in range(6):
        value: Polynomial = {}
        for index in range(order + 1):
            term = multiply_polynomials(
                derivatives[order - index],
                q_variables[index],
            )
            scalar = comb(order, index) * (-1) ** index * factorial(index)
            value = add_polynomials(value, scale_polynomial(term, scalar))
        derivatives.append(value)

    matrix = [
        [derivatives[3 + column - row] for column in range(4)]
        for row in range(4)
    ]
    determinant = determinant_polynomial(matrix)
    assert all(exponents[0] == 0 for exponents in determinant)
    assert all(coefficient % 12 == 0 for coefficient in determinant.values())
    return {
        exponents: coefficient // 12
        for exponents, coefficient in determinant.items()
    }


def substitute_power_sums(polynomial: Polynomial) -> Polynomial:
    root_power_sums: list[Polynomial] = []
    for degree in range(1, Q_VARIABLES + 1):
        power_sum: Polynomial = {}
        for root_index in range(ROOT_VARIABLES):
            exponents = [0] * ROOT_VARIABLES
            exponents[root_index] = degree
            power_sum[tuple(exponents)] = 1
        root_power_sums.append(power_sum)

    result: Polynomial = {}
    for q_exponents, coefficient in polynomial.items():
        term = {ZERO_ROOT_EXPONENT: coefficient}
        for index, exponent in enumerate(q_exponents):
            if exponent:
                term = multiply_polynomials(
                    term,
                    power_polynomial(root_power_sums[index], exponent),
                )
        result = add_polynomials(result, term)
    return result


def orbit_coefficients(polynomial: Polynomial) -> dict[tuple[int, ...], int]:
    result: dict[tuple[int, ...], int] = {}
    counts: defaultdict[tuple[int, ...], int] = defaultdict(int)
    for exponents, coefficient in polynomial.items():
        partition = tuple(
            sorted((power for power in exponents if power), reverse=True)
        )
        previous = result.setdefault(partition, coefficient)
        assert previous == coefficient
        counts[partition] += 1
    for partition, count in counts.items():
        padded = partition + (0,) * (ROOT_VARIABLES - len(partition))
        assert count == len(set(permutations(padded)))
    return result


def verify_am_gm() -> None:
    orbit = set(permutations((4, 3, 3, 2)))
    assert len(orbit) == 12
    assert [
        sum(exponents[index] for exponents in orbit)
        for index in range(4)
    ] == [36, 36, 36, 36]


def verify_johnson_certificate() -> None:
    triples = list(combinations(range(12), 3))
    scaled_l_values = {3: 252, 2: -28, 1: 7, 0: -3}
    scaled_k_values = {3: 252, 2: -28, 1: 28, 0: 84}

    scaled_l = [
        [
            scaled_l_values[len(set(left).intersection(right))]
            for right in triples
        ]
        for left in triples
    ]
    assert all(
        scaled_l[row][column] == scaled_l[column][row]
        for row in range(len(triples))
        for column in range(len(triples))
    )

    for row in range(len(triples)):
        for column in range(len(triples)):
            square_entry = sum(
                scaled_l[row][middle] * scaled_l[middle][column]
                for middle in range(len(triples))
            )
            assert square_entry == 360 * scaled_l[row][column]

    assert sum(scaled_l[index][index] for index in range(len(triples))) == (
        154 * 360
    )
    assert all(
        scaled_k_values[intersection] - scaled_l_values[intersection] >= 0
        for intersection in range(4)
    )

    q_coefficients: defaultdict[tuple[int, ...], int] = defaultdict(int)
    for left in triples:
        left_set = set(left)
        for right in triples:
            exponents = [0] * 12
            for index in left:
                exponents[index] += 1
            for index in right:
                exponents[index] += 1
            q_coefficients[tuple(exponents)] += {
                3: 9,
                2: -1,
                1: 1,
                0: 3,
            }[len(left_set.intersection(right))]

    orbit_values: dict[tuple[int, ...], int] = {}
    for exponents, coefficient in q_coefficients.items():
        partition = tuple(
            sorted((power for power in exponents if power), reverse=True)
        )
        previous = orbit_values.setdefault(partition, coefficient)
        assert previous == coefficient
    assert orbit_values == EXPECTED_Q_ORBITS


def bareiss_determinant(matrix: list[list[int | Fraction]]) -> int | Fraction:
    work = [
        [Fraction(entry) for entry in row]
        for row in matrix
    ]
    size = len(work)
    sign = 1
    previous = Fraction(1)
    for pivot_index in range(size - 1):
        if work[pivot_index][pivot_index] == 0:
            swap_index = next(
                (
                    index
                    for index in range(pivot_index + 1, size)
                    if work[index][pivot_index] != 0
                ),
                None,
            )
            if swap_index is None:
                return 0
            work[pivot_index], work[swap_index] = (
                work[swap_index],
                work[pivot_index],
            )
            sign *= -1
        pivot = work[pivot_index][pivot_index]
        for row in range(pivot_index + 1, size):
            for column in range(pivot_index + 1, size):
                work[row][column] = (
                    work[row][column] * pivot
                    - work[row][pivot_index] * work[pivot_index][column]
                ) / previous
        previous = pivot
    return sign * work[-1][-1]


def convolve(
    left: list[int | Fraction],
    right: list[int | Fraction],
) -> list[int | Fraction]:
    result: list[int | Fraction] = [0] * (len(left) + len(right) - 1)
    for left_index, left_coefficient in enumerate(left):
        for right_index, right_coefficient in enumerate(right):
            result[left_index + right_index] += (
                left_coefficient * right_coefficient
            )
    return result


def product_coefficients(
    reciprocal_roots: list[int | Fraction],
) -> list[int | Fraction]:
    coefficients: list[int | Fraction] = [1]
    for root in reciprocal_roots:
        coefficients = convolve(coefficients, [1, root])
    return coefficients


def centered_product_determinant(
    reciprocal_roots: list[int | Fraction],
) -> int | Fraction:
    coefficients = product_coefficients(reciprocal_roots)
    derivatives = [
        factorial(degree) * (
            coefficients[degree] if degree < len(coefficients) else 0
        )
        for degree in range(7)
    ]
    matrix = [
        [
            derivatives[3 + column - row]
            if 3 + column - row >= 0
            else 0
            for column in range(4)
        ]
        for row in range(4)
    ]
    return bareiss_determinant(matrix)


def polynomial_wronskian_at_zero(
    coefficients: list[int | Fraction],
) -> int | Fraction:
    derivatives = [
        factorial(degree) * (
            coefficients[degree] if degree < len(coefficients) else 0
        )
        for degree in range(7)
    ]
    return bareiss_determinant(
        [
            [derivatives[row + column] for column in range(4)]
            for row in range(4)
        ]
    )


def restricted_family_coefficients(
    roots: list[int | Fraction],
) -> list[int | Fraction]:
    psi: list[int | Fraction] = [1]
    for root in roots:
        psi = convolve(psi, [root, 1])
    return [0, 0] + psi


def verify_published_and_restricted_witnesses() -> None:
    published_coefficients = [0, 0] + [
        comb(11, degree) for degree in range(12)
    ]
    assert polynomial_wronskian_at_zero(published_coefficients) == -1584
    assert centered_product_determinant([162, 162] + [1] * 11) == (
        EXPECTED_STRICT_WITNESS
    )
    assert centered_product_determinant([161, 161] + [1] * 11) > 0

    for roots in (
        [1] * 11,
        [2] * 12,
        [Fraction(1, 2), Fraction(2, 3), Fraction(3, 5)],
        [1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
    ):
        psi_zero = 1
        for root in roots:
            psi_zero *= root
        reciprocals = [Fraction(1, root) for root in roots]
        b_value = sum(value**2 for value in reciprocals)
        d_value = sum(value**4 for value in reciprocals)
        source_formula = 144 * psi_zero**4 * (10 * d_value - b_value**2)
        determinant = polynomial_wronskian_at_zero(
            restricted_family_coefficients(roots)
        )
        assert determinant == source_formula

    assert 10 * 11 - 11**2 < 0
    assert all(10 * count - count**2 >= 0 for count in range(11))
    assert polynomial_wronskian_at_zero(
        restricted_family_coefficients([2] * 12)
    ) == -3456 * 2**44


def verify_transfer_identity() -> None:
    roots = [0, 1, 2, 5]
    point = 3
    coefficients: list[int | Fraction] = [1]
    for root in roots:
        coefficients = convolve(coefficients, [root, 1])
    derivatives = []
    for order in range(7):
        value = 0
        for degree in range(order, len(coefficients)):
            value += (
                coefficients[degree]
                * factorial(degree)
                // factorial(degree - order)
                * point ** (degree - order)
            )
        derivatives.append(value)
    hankel = bareiss_determinant(
        [
            [derivatives[row + column] for column in range(4)]
            for row in range(4)
        ]
    )
    polynomial_value = derivatives[0]
    reciprocal_roots = [
        Fraction(1, point + root)
        for root in roots
    ]
    assert hankel == (
        polynomial_value**4
        * centered_product_determinant(reciprocal_roots)
    )
    assert (-1) ** (4 * 3 // 2) == 1


def verify_scope_and_equality_boundary() -> None:
    assert centered_product_determinant([]) == 0
    assert centered_product_determinant([1]) == 0
    assert centered_product_determinant([1, 2]) == 0
    assert centered_product_determinant([1, 2, 3]) > 0
    assert centered_product_determinant([1, 2, 3, 4]) > 0

    structural_boundary = [
        [16, 128, 0, 0],
        [1, 16, 128, 0],
        [0, 1, 16, 128],
        [0, 0, 1, 16],
    ]
    assert bareiss_determinant(structural_boundary) == -16_384


def sha256(path: Path) -> str:
    digest = hashlib.sha256()
    with path.open("rb") as handle:
        for block in iter(lambda: handle.read(1 << 20), b""):
            digest.update(block)
    return digest.hexdigest()


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

    source_2002 = (
        source_directory / "craven-csordas-2002-iterated-laguerre-turan.txt"
    ).read_text()
    source_2005 = (
        source_directory / "craven-csordas-2005-karlin-conjecture.txt"
    ).read_text()
    for fragment in (
        "Remarks 5.6. (a)",
        "K(x) = 3(10 D − B 2 )",
        "such that 10 D < B 2 holds at x = 0",
        "x2 (x + a)12",
    ):
        assert fragment in source_2002
    for fragment in (
        "Example 3.3.",
        "The lowest",
        "degree example we have found",
        "Wronskian at x = 0 is −1584.",
    ):
        assert fragment in source_2005


def main() -> None:
    verify_sources()
    power_sum_determinant = normalized_power_sum_determinant()
    root_determinant = substitute_power_sums(power_sum_determinant)
    assert orbit_coefficients(root_determinant) == EXPECTED_ORBITS
    verify_am_gm()
    verify_johnson_certificate()
    verify_published_and_restricted_witnesses()
    verify_transfer_identity()
    verify_scope_and_equality_boundary()
    print("primary-source hashes and passages: PASS")
    print("dependency-free determinant orbit expansion: PASS")
    print("AM-GM orbit: PASS")
    print("direct 220x220 Johnson identity A^2=360A: PASS")
    print("published -1584 and strict-root normalization: PASS")
    print("restricted-family degree-13 threshold: PASS")
    print("polynomial transfer and endpoint setup: PASS")
    print("equality and full-PF4 scope boundary: PASS")


if __name__ == "__main__":
    main()
