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

from __future__ import annotations

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


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
)
EXPECTED_SCOPE_COUNTERFEIT = -222_236_584_443_224_648_570_880
EXPECTED_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,
}
EXPECTED_SOURCE_SHA256 = {
    "cc2002-page16.png":
        "6c7b50158ff09e43285c05ad6e3ecead08c8b708d2bd624d2e4d33884713438f",
    "cc2002-page17.png":
        "bf1e3c33ab0130949222ab1ad6b445d9c07f4ffa1a864926c68c5977dcc4a993",
    "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":
        "81705c43d0db60f5d5dbed189ae8220bb4943cf307407cf74f3bee82c8ec4623",
}

Exponent = tuple[int, ...]
Polynomial = dict[Exponent, int]
ZERO_EXPONENT = (0,) * VARIABLES
MODULUS = 2**61 - 1


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


def add_polynomials(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 {
        exponents: coefficient
        for exponents, coefficient in result.items()
        if coefficient
    }


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


def multiply_polynomials(left: Polynomial, right: Polynomial) -> Polynomial:
    result: defaultdict[Exponent, 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 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:
    result: Polynomial = {}
    for permutation in permutations(range(len(matrix))):
        term = {ZERO_EXPONENT: 1}
        for row, column in enumerate(permutation):
            term = multiply_polynomials(term, matrix[row][column])
            if not term:
                break
        result = add_polynomials(
            result,
            scale_polynomial(term, permutation_sign(permutation)),
        )
    return result


def elementary_determinant() -> Polynomial:
    elementary = [{ZERO_EXPONENT: 1}]
    for index in range(VARIABLES):
        exponents = [0] * VARIABLES
        exponents[index] = 1
        elementary.append({tuple(exponents): 1})

    matrix: list[list[Polynomial]] = []
    for row in range(5):
        matrix_row = []
        for column in range(5):
            degree = 7 - row - column
            if degree < 0:
                matrix_row.append({})
            else:
                matrix_row.append(
                    scale_polynomial(
                        elementary[degree],
                        factorial(row + column),
                    )
                )
        matrix.append(matrix_row)

    determinant = determinant_polynomial(matrix)
    content = 0
    for coefficient in determinant.values():
        content = gcd(content, abs(coefficient))
    assert content == DETERMINANT_CONTENT
    assert all(
        coefficient % DETERMINANT_CONTENT == 0
        for coefficient in determinant.values()
    )
    return {
        exponents: coefficient // DETERMINANT_CONTENT
        for exponents, coefficient in determinant.items()
    }


def bounded_partitions(
    total: int,
    maximum_part: int,
    maximum_length: int,
    prefix: Exponent = (),
) -> list[Exponent]:
    if total == 0:
        return [prefix]
    if maximum_length == 0:
        return []
    result = []
    for part in range(min(total, maximum_part), 0, -1):
        result.extend(
            bounded_partitions(
                total - part,
                part,
                maximum_length - 1,
                prefix + (part,),
            )
        )
    return result


def orbit_permutations(partition: Exponent) -> tuple[Exponent, ...]:
    padded = partition + (0,) * (VARIABLES - len(partition))
    return tuple(set(permutations(padded)))


def orbit_descriptor(
    partition: Exponent,
) -> tuple[tuple[int, ...], tuple[int, ...]]:
    padded = partition + (0,) * (VARIABLES - len(partition))
    exponents = tuple(sorted(set(padded)))
    counts = tuple(padded.count(exponent) for exponent in exponents)
    return exponents, counts


def monomial_symmetric_value(
    point: tuple[int, ...],
    descriptor: tuple[tuple[int, ...], tuple[int, ...]],
) -> int:
    exponents, initial_counts = descriptor
    states = {initial_counts: 1}
    for value in point:
        updated: defaultdict[tuple[int, ...], int] = defaultdict(int)
        powers = {
            exponent: value**exponent
            for exponent in exponents
        }
        for counts, coefficient in states.items():
            for index, count in enumerate(counts):
                if count:
                    next_counts = list(counts)
                    next_counts[index] -= 1
                    updated[tuple(next_counts)] += (
                        coefficient * powers[exponents[index]]
                    )
        states = dict(updated)
    return states[(0,) * len(exponents)]


def monomial_symmetric_row(
    point: tuple[int, ...],
    descriptors: list[tuple[tuple[int, ...], tuple[int, ...]]],
) -> list[int]:
    return [
        monomial_symmetric_value(point, descriptor)
        for descriptor in descriptors
    ]


def elementary_values(point: tuple[int, ...]) -> list[int]:
    values = [0] * (VARIABLES + 1)
    values[0] = 1
    for root in point:
        for degree in range(VARIABLES, 0, -1):
            values[degree] += root * values[degree - 1]
    return values


def evaluate_elementary_polynomial(
    polynomial: Polynomial,
    point: tuple[int, ...],
) -> int:
    values = elementary_values(point)
    return sum(
        coefficient
        * product(
            values[index + 1] ** exponent
            for index, exponent in enumerate(exponents)
        )
        for exponents, coefficient in polynomial.items()
    )


def add_modular_rank_row(
    row: list[int],
    echelon: dict[int, list[int]],
) -> bool:
    reduced = [value % MODULUS for value in row]
    for pivot in sorted(echelon):
        factor = reduced[pivot]
        if factor:
            basis = echelon[pivot]
            reduced = [
                (value - factor * basis_value) % MODULUS
                for value, basis_value in zip(
                    reduced,
                    basis,
                    strict=True,
                )
            ]
    pivot = next(
        (index for index, value in enumerate(reduced) if value),
        None,
    )
    if pivot is None:
        return False
    inverse = pow(reduced[pivot], -1, MODULUS)
    echelon[pivot] = [
        value * inverse % MODULUS
        for value in reduced
    ]
    return True


def solve_modular_system(
    matrix: list[list[int]],
    right_hand_side: list[int],
) -> list[int]:
    size = len(matrix)
    augmented = [
        [value % MODULUS for value in row]
        + [right_hand_side[index] % MODULUS]
        for index, row in enumerate(matrix)
    ]
    for column in range(size):
        pivot = next(
            row
            for row in range(column, size)
            if augmented[row][column]
        )
        augmented[column], augmented[pivot] = (
            augmented[pivot],
            augmented[column],
        )
        inverse = pow(augmented[column][column], -1, MODULUS)
        augmented[column] = [
            value * inverse % MODULUS
            for value in augmented[column]
        ]
        for row in range(size):
            if row == column:
                continue
            factor = augmented[row][column]
            if factor:
                augmented[row] = [
                    (value - factor * pivot_value) % MODULUS
                    for value, pivot_value in zip(
                        augmented[row],
                        augmented[column],
                        strict=True,
                    )
                ]
    return [augmented[index][-1] for index in range(size)]


def derive_source_orbits(
    elementary: Polynomial,
) -> tuple[dict[Exponent, int], Polynomial]:
    assert all(
        sum((index + 1) * exponent for index, exponent in enumerate(exponents))
        == 15
        for exponents in elementary
    )
    assert all(sum(exponents) <= 5 for exponents in elementary)

    partitions = bounded_partitions(15, 5, VARIABLES)
    assert len(partitions) == 46
    descriptors = [orbit_descriptor(partition) for partition in partitions]

    selected_rows: list[list[int]] = []
    selected_values: list[int] = []
    echelon: dict[int, list[int]] = {}
    for point in combinations_with_replacement(range(8), VARIABLES):
        row = monomial_symmetric_row(point, descriptors)
        if add_modular_rank_row(row, echelon):
            selected_rows.append(row)
            selected_values.append(
                evaluate_elementary_polynomial(elementary, point)
            )
            if len(selected_rows) == len(partitions):
                break
    assert len(selected_rows) == len(partitions)

    residues = solve_modular_system(selected_rows, selected_values)
    coefficients = [
        residue if residue <= MODULUS // 2 else residue - MODULUS
        for residue in residues
    ]
    assert all(
        sum(
            row[index] * coefficients[index]
            for index in range(len(partitions))
        ) == value
        for row, value in zip(
            selected_rows,
            selected_values,
            strict=True,
        )
    )

    orbit_coefficients_result = {
        partition: coefficient
        for partition, coefficient in zip(
            partitions,
            coefficients,
            strict=True,
        )
        if coefficient
    }
    source: Polynomial = {}
    for partition, coefficient in orbit_coefficients_result.items():
        for exponents in orbit_permutations(partition):
            source[exponents] = coefficient
    return orbit_coefficients_result, source


def orbit_coefficients(polynomial: Polynomial) -> dict[Exponent, int]:
    result: dict[Exponent, int] = {}
    counts: defaultdict[Exponent, 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,) * (VARIABLES - len(partition))
        assert count == len(set(permutations(padded)))
    return result


def ordered_difference_transform(
    coefficients: Polynomial,
) -> Polynomial:
    transformed = coefficients
    for index in range(VARIABLES - 1):
        updated: defaultdict[Exponent, 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: Polynomial) -> 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: Polynomial,
    values: tuple[int | Fraction, ...],
) -> int | Fraction:
    return sum(
        coefficient
        * product(
            value**exponent
            for value, exponent in zip(values, exponents, strict=True)
        )
        for exponents, coefficient in coefficients.items()
    )


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


def verify_transform(source: Polynomial, chamber: Polynomial) -> None:
    assert len(source) == EXPECTED_SOURCE_TERMS
    assert len(chamber) == EXPECTED_CHAMBER_TERMS
    assert all(sum(exponents) == 15 for exponents in source)
    assert all(sum(exponents) == 15 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, 0, 15)]
        == EXPECTED_CONSTANT_RAY_COEFFICIENT
    )
    assert dictionary_digest(chamber) == EXPECTED_CHAMBER_SHA256

    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 bareiss_determinant(
    matrix: list[list[int | Fraction]],
) -> int | Fraction:
    work = [
        [Fraction(entry) for entry in row]
        for row in matrix
    ]
    sign = 1
    previous = Fraction(1)
    for pivot_index in range(len(work) - 1):
        if work[pivot_index][pivot_index] == 0:
            swap_index = next(
                (
                    index
                    for index in range(pivot_index + 1, len(work))
                    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, len(work)):
            for column in range(pivot_index + 1, len(work)):
                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: tuple[int | Fraction, ...],
) -> list[int | Fraction]:
    coefficients: list[int | Fraction] = [1]
    for root in reciprocal_roots:
        coefficients = convolve(coefficients, [1, root])
    return coefficients


def centered_determinant(
    reciprocal_roots: tuple[int | Fraction, ...],
) -> int | Fraction:
    coefficients = product_coefficients(reciprocal_roots)
    derivatives = [
        factorial(degree) * (
            coefficients[degree] if degree < len(coefficients) else 0
        )
        for degree in range(9)
    ]
    matrix = [
        [
            derivatives[4 + column - row]
            if 4 + column - row >= 0
            else 0
            for column in range(5)
        ]
        for row in range(5)
    ]
    return bareiss_determinant(matrix)


def polynomial_hankel_at_zero(
    coefficients: tuple[int | Fraction, ...],
) -> int | Fraction:
    derivatives = [
        factorial(degree) * (
            coefficients[degree] if degree < len(coefficients) else 0
        )
        for degree in range(9)
    ]
    return bareiss_determinant(
        [
            [derivatives[row + column] for column in range(5)]
            for row in range(5)
        ]
    )


def verify_reciprocal_relation(source: Polynomial) -> None:
    alpha = (1, 2, 3, 4, 5, 6, 7)
    beta = tuple(Fraction(1, value) for value in alpha)
    assert centered_determinant(alpha) == (
        DETERMINANT_CONTENT
        * product(alpha) ** 5
        * evaluate_dictionary(source, beta)
    )


def verify_published_witnesses() -> None:
    boundary_coefficients = (0, 0, 0) + tuple(
        comb(5, degree) for degree in range(6)
    )
    assert polynomial_hankel_at_zero(boundary_coefficients) == (
        EXPECTED_BOUNDARY_WITNESS
    )
    assert centered_determinant(
        (50, 50, 50, 1, 1, 1, 1, 1)
    ) == EXPECTED_DEGREE_EIGHT_WITNESS


def polynomial_derivatives_at(
    coefficients: list[int | Fraction],
    point: int,
) -> list[int | Fraction]:
    derivatives = []
    for order in range(9):
        value: int | Fraction = 0
        for degree in range(order, len(coefficients)):
            value += (
                coefficients[degree]
                * factorial(degree)
                // factorial(degree - order)
                * point ** (degree - order)
            )
        derivatives.append(value)
    return derivatives


def verify_transfer_and_sign() -> None:
    roots = (0, 1, 2, 4, 7, 8, 11)
    coefficients: list[int | Fraction] = [1]
    for root in roots:
        coefficients = convolve(coefficients, [root, 1])
    derivatives = polynomial_derivatives_at(coefficients, 3)
    hankel = bareiss_determinant(
        [
            [derivatives[row + column] for column in range(5)]
            for row in range(5)
        ]
    )
    reciprocal_roots = tuple(
        Fraction(1, 3 + root) for root in roots
    )
    assert hankel == derivatives[0] ** 5 * centered_determinant(
        reciprocal_roots
    )
    negative_hankel = bareiss_determinant(
        [
            [-derivatives[row + column] for column in range(5)]
            for row in range(5)
        ]
    )
    assert negative_hankel == -hankel
    assert (-1) ** (5 * 4 // 2) == 1


def verify_continuity_and_scope() -> None:
    for support in range(8):
        roots = tuple(range(1, support + 1))
        assert centered_determinant(roots) >= 0
    assert centered_determinant(
        (1,) * 10 + (13, 14, 14)
    ) == EXPECTED_SCOPE_COUNTERFEIT


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 )",
    ):
        assert fragment in source_2002
    for fragment in (
        "Example 3.1. Let f (x) := x3 (x + 1)5",
        "= 414720(x+1)5",
        "+ 101015x3 + 12540x2 + 600x − 24),",
        "One can avoid the root at 0; (50x + 1)3 (x + 1)5",
        "perturbation of the zeros will give examples with simple zeros.",
    ):
        assert fragment in source_2005


def main() -> None:
    verify_sources()
    elementary = elementary_determinant()
    recovered_orbits, source = derive_source_orbits(elementary)
    assert recovered_orbits == EXPECTED_ORBITS
    assert orbit_coefficients(source) == EXPECTED_ORBITS
    chamber = ordered_difference_transform(source)
    verify_transform(source, chamber)
    verify_reciprocal_relation(source)
    verify_published_witnesses()
    verify_transfer_and_sign()
    verify_continuity_and_scope()
    print("primary-source hashes and passages: PASS")
    print("dependency-free reciprocal determinant and 20 orbits: 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("published boundary and strict-root witnesses: PASS")
    print("positive-leading transfer and sign reversal: PASS")
    print("continuity support and PF5 scope counterfeit: PASS")


if __name__ == "__main__":
    main()
