#!/usr/bin/env python3
# SPDX-License-Identifier: MIT
"""Pure-standard-library exact verifier for ``yagzhev19_map.json``.

The JSON file stores H for G = X + H, three rational points, and their common
image. This script uses only ``fractions.Fraction`` and dictionary polynomials.
It checks:

  * every nonzero H_i is homogeneous of degree 3;
  * the three distinct rational points have the displayed common image;
  * (JH)^18 is the zero polynomial matrix;
  * ((JH)^17)_[4,19] = -18 t^16 x^10 y^6 z^2.

Paper indices are one-based; the corresponding Python matrix key is (3, 18).

No computer-algebra package is imported.

Canonical release: https://harrischan.com/docs/yagzhev-19/

Acknowledgment: Claude Fable 5, GPT-5.6 Sol, and Gemini 3.1 Pro assisted with
mathematical exploration, drafting, code generation, and cross-checking. They
are credited as computational tools rather than authors. Harris Chan takes
responsibility for the mathematical claims, citations, attribution decisions,
and conclusions.
"""

from __future__ import annotations

import hashlib
import json
import sys
from collections import defaultdict
from fractions import Fraction
from pathlib import Path
from typing import Dict, Iterable, Tuple

Monomial = Tuple[int, ...]
Polynomial = Dict[Monomial, Fraction]
Matrix = Dict[Tuple[int, int], Polynomial]

EXPECTED_SHA256 = "f88a0883ce9509339c53a6a721390053b424039a8dd6c41ebc7770aeb955cced"
EXPECTED_VARIABLE_ORDER = [
    "x", "y", "z", "a", "b", "c", "d", "q", "s", "h", "k",
    "Y1", "Y2", "Y3", "Y4", "Y5", "Y6", "Y9", "t",
]


def parse_fraction(value: str) -> Fraction:
    return Fraction(value)


def add_monomials(left: Monomial, right: Monomial) -> Monomial:
    return tuple(a + b for a, b in zip(left, right))


def poly_mul(left: Polynomial, right: Polynomial) -> Polynomial:
    if not left or not right:
        return {}
    output: defaultdict[Monomial, Fraction] = defaultdict(Fraction)
    for mon_left, coeff_left in left.items():
        for mon_right, coeff_right in right.items():
            output[add_monomials(mon_left, mon_right)] += coeff_left * coeff_right
    return {mon: coeff for mon, coeff in output.items() if coeff}


def poly_accumulate(target: defaultdict[Monomial, Fraction], source: Polynomial) -> None:
    for monomial, coefficient in source.items():
        target[monomial] += coefficient


def poly_derivative(poly: Polynomial, variable: int) -> Polynomial:
    output: Polynomial = {}
    for monomial, coefficient in poly.items():
        exponent = monomial[variable]
        if exponent == 0:
            continue
        derived = list(monomial)
        derived[variable] -= 1
        output[tuple(derived)] = coefficient * exponent
    return output


def poly_evaluate(poly: Polynomial, point: Tuple[Fraction, ...]) -> Fraction:
    result = Fraction(0)
    for monomial, coefficient in poly.items():
        term = coefficient
        for value, exponent in zip(point, monomial):
            if exponent:
                term *= value**exponent
        result += term
    return result


def matrix_mul(left: Matrix, right: Matrix) -> Matrix:
    right_by_row: defaultdict[int, list[Tuple[int, Polynomial]]] = defaultdict(list)
    for (row, column), polynomial in right.items():
        right_by_row[row].append((column, polynomial))

    accumulators: dict[Tuple[int, int], defaultdict[Monomial, Fraction]] = {}
    for (row, middle), left_poly in left.items():
        for column, right_poly in right_by_row.get(middle, []):
            product = poly_mul(left_poly, right_poly)
            if not product:
                continue
            key = (row, column)
            accumulator = accumulators.setdefault(key, defaultdict(Fraction))
            poly_accumulate(accumulator, product)

    output: Matrix = {}
    for key, accumulator in accumulators.items():
        polynomial = {mon: coeff for mon, coeff in accumulator.items() if coeff}
        if polynomial:
            output[key] = polynomial
    return output


def load_map(path: Path) -> tuple[dict, bytes]:
    raw = path.read_bytes()
    digest = hashlib.sha256(raw).hexdigest()
    if digest != EXPECTED_SHA256:
        raise AssertionError(
            f"unexpected JSON SHA-256: {digest}; expected {EXPECTED_SHA256}"
        )
    return json.loads(raw), raw


def main() -> None:
    path = Path(sys.argv[1]) if len(sys.argv) > 1 else Path(__file__).with_name("yagzhev19_map.json")
    data, _ = load_map(path)

    dimension = int(data["dimension"])
    assert dimension == 19
    assert data["variable_order"] == EXPECTED_VARIABLE_ORDER
    assert len(data["H"]) == dimension

    H: list[Polynomial] = []
    for component in data["H"]:
        polynomial: Polynomial = {}
        for term in component:
            monomial = tuple(int(exponent) for exponent in term["exponents"])
            coefficient = parse_fraction(term["coeff"])
            assert len(monomial) == dimension
            polynomial[monomial] = polynomial.get(monomial, Fraction(0)) + coefficient
        polynomial = {mon: coeff for mon, coeff in polynomial.items() if coeff}
        if polynomial:
            assert all(sum(monomial) == 3 for monomial in polynomial)
        H.append(polynomial)

    points = [tuple(parse_fraction(value) for value in row) for row in data["collision_points"]]
    common_image = tuple(parse_fraction(value) for value in data["common_image"])
    assert len(points) == 3 and len(set(points)) == 3
    assert all(len(point) == dimension for point in points)
    assert len(common_image) == dimension

    for point in points:
        image = tuple(point[i] + poly_evaluate(H[i], point) for i in range(dimension))
        assert image == common_image

    JH: Matrix = {}
    for row, component in enumerate(H):
        for column in range(dimension):
            derivative = poly_derivative(component, column)
            if derivative:
                JH[(row, column)] = derivative

    power = JH
    power_17: Matrix | None = None
    for exponent in range(1, 19):
        if exponent == 17:
            power_17 = {key: dict(poly) for key, poly in power.items()}
        if exponent < 18:
            power = matrix_mul(power, JH)

    assert power_17 is not None and power_17
    assert not power

    expected_monomial = [0] * dimension
    expected_monomial[0] = 10  # x
    expected_monomial[1] = 6   # y
    expected_monomial[2] = 2   # z
    expected_monomial[18] = 16 # t
    expected_entry = {tuple(expected_monomial): Fraction(-18)}
    assert power_17[(3, 18)] == expected_entry

    print("Verified yagzhev19_map.json using only the Python standard library:")
    print("  SHA-256:", EXPECTED_SHA256)
    print("  dimension: 19")
    print("  every nonzero H_i is cubic homogeneous")
    print("  three distinct rational points share the displayed image")
    print("  (JH)^18 = 0")
    print("  ((JH)^17)_[4,19] = -18*t^16*x^10*y^6*z^2")
    print("  therefore det J(X+H) = 1 and X+H is non-injective")


if __name__ == "__main__":
    main()