#!/usr/bin/env python3
# SPDX-License-Identifier: MIT
r"""
AN EXPLICIT CUBIC-HOMOGENEOUS JACOBIAN COUNTEREXAMPLE IN 19 VARIABLES
=====================================================================

Result
------
This script constructs an explicit map

    Ghat = X + H : C^19 -> C^19

and verifies, in exact rational polynomial arithmetic, that

    * every nonzero component of H is homogeneous of degree 3;
    * (JH)^18 = 0 but (JH)^17 != 0;
    * three distinct rational points have the same image.

It follows that det J(Ghat) = 1 identically and Ghat is not injective.

Attribution and relationship to public constructions
----------------------------------------------------
1. Levent Alpoge publicly posted the three-variable Keller map from which the
   immediate reductions originate (July 20, 2026 UTC; the announcement credits
   Fable 5 for finding the example):
       https://x.com/__alpoge__/status/2079028340955197566

2. Ali Demirci (depyronick) publicly posted a direct, generic implementation of
   the classical degree-reduction and homogenization pipeline, producing a
   Yagzhev-form example in dimension 79:
       https://x.com/alidemirci/status/2079304070297477493
       https://gist.github.com/depyronick/d20a019d915c8b3a612a45e0846dede8

3. Reddit user u/Distinct-Soft-3991 posted an explicit degree-three reduction
   in dimension 17 in the immediate discussion:
       https://www.reddit.com/r/math/comments/1v1aix1/

4. Joe Atkins-Turkish (GitHub/Reddit: Spacerat) posted the factor-aware
   degree-three map in dimension 11 used as the input below. The Gist states
   that the explicit construction, simplification, and verification code were
   generated by ChatGPT:
       https://gist.github.com/Spacerat/08b4a43f6b6ca57178efabc220170ce8

5. The dimension-19 lift and its verification were developed through
   iterative prompting and symbolic computation. Harris Chan curated the
   resulting construction, checked the generated artifacts, and prepared the
   exposition.

AI assistance
-------------
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.

This file is intended to make the final object independently checkable: its
proof does not require trusting the derivation history once the displayed
polynomial map has been constructed.

The classical reductions are due to Yagzhev, Bass-Connell-Wright,
Druzkowski, and the formulation used here follows Campbell:
    https://arxiv.org/abs/1303.3853

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

Tested environment
------------------
    Python 3.13.5
    SymPy 1.14.0

Typical runtime is dominated by exact multiplication of the 19 x 19
polynomial Jacobian matrix through the eighteenth power.
"""

from __future__ import annotations

from collections import defaultdict
import hashlib
import json
from pathlib import Path
import sys
from typing import Dict, Iterable, Mapping, Sequence, Tuple

import sympy as sp

EXPECTED_JSON_SHA256 = "f88a0883ce9509339c53a6a721390053b424039a8dd6c41ebc7770aeb955cced"

PolyMatrix = Dict[Tuple[int, int], sp.Expr]


def homogeneous_part(
    expr: sp.Expr,
    variables: Sequence[sp.Symbol],
    degree: int,
) -> sp.Expr:
    """Return the homogeneous part of ``expr`` of the requested total degree."""
    poly = sp.Poly(sp.expand(expr), *variables)
    result = sp.Integer(0)
    for monomial, coefficient in poly.terms():
        if sum(monomial) != degree:
            continue
        term = coefficient
        for variable, exponent in zip(variables, monomial):
            if exponent:
                term *= variable**exponent
        result += term
    return sp.expand(result)


def evaluate_map(
    components: Sequence[sp.Expr],
    variables: Sequence[sp.Symbol],
    point: Sequence[sp.Expr],
) -> Tuple[sp.Expr, ...]:
    if len(variables) != len(point):
        raise ValueError("Point dimension does not match the map dimension")
    substitution = dict(zip(variables, point))
    return tuple(sp.expand(component.subs(substitution)) for component in components)


def sparse_jacobian(
    components: Sequence[sp.Expr],
    variables: Sequence[sp.Symbol],
) -> PolyMatrix:
    result: PolyMatrix = {}
    for row, component in enumerate(components):
        if component == 0:
            continue
        for column, variable in enumerate(variables):
            derivative = sp.diff(component, variable)
            if derivative != 0:
                result[(row, column)] = sp.expand(derivative)
    return result


def sparse_multiply(left: PolyMatrix, right: PolyMatrix) -> PolyMatrix:
    """Multiply sparse polynomial matrices exactly, cancelling zero entries."""
    right_by_row: Dict[int, list[Tuple[int, sp.Expr]]] = defaultdict(list)
    for (row, column), value in right.items():
        right_by_row[row].append((column, value))

    accumulated: Dict[Tuple[int, int], sp.Expr] = defaultdict(lambda: sp.Integer(0))
    for (row, middle), left_value in left.items():
        for column, right_value in right_by_row.get(middle, []):
            accumulated[(row, column)] += left_value * right_value

    product: PolyMatrix = {}
    for index, value in accumulated.items():
        expanded = sp.expand(value)
        if expanded != 0:
            product[index] = expanded
    return product


# ---------------------------------------------------------------------------
# 1. Explicit degree-three map Phi : C^11 -> C^11 posted by Joe Atkins-Turkish.
# ---------------------------------------------------------------------------

x, y, z, a, b, c, d, q, s, h, k = sp.symbols("x y z a b c d q s h k")
variables_11 = [x, y, z, a, b, c, d, q, s, h, k]

phi = [
    -a*c - a*d*z - 3*a*y**2 - 2*a*z
    - c*d**2 + d**2*z - d*s + 7*d*y**2
    + s*x*y + 3*x*y*z + 4*y**2 + z,

    -b*c - b*d*z - 3*b*y**2 - 2*b*z
    - 3*c*d*x - d*q + q*x*y
    + 12*x*y**2 + 3*x*z + y,

    -h*k - h*x*z + k*x**2 - 3*x**2*y + 2*x,

    a - d**2 + 2*d*x*y,
    b + 3*x**2*y,
    c + x*y*z + 3*y**2 + 2*z,
    d - x*y,
    b*z + 3*c*x + q,
    s + a*z + c*x*y - x*y*z - 7*y**2 + c*d - d*z,
    h - x**2,
    k + x*z,
]
phi = [sp.expand(component) for component in phi]

points_11 = [
    (
        sp.Rational(0), sp.Rational(0), sp.Rational(-1, 4),
        sp.Rational(0), sp.Rational(0), sp.Rational(1, 2),
        sp.Rational(0), sp.Rational(0), sp.Rational(0),
        sp.Rational(0), sp.Rational(0),
    ),
    (
        sp.Rational(1), sp.Rational(-3, 2), sp.Rational(13, 2),
        sp.Rational(-9, 4), sp.Rational(9, 2), sp.Rational(-10),
        sp.Rational(-3, 2), sp.Rational(3, 4), sp.Rational(-153, 8),
        sp.Rational(1), sp.Rational(-13, 2),
    ),
    (
        sp.Rational(-1), sp.Rational(3, 2), sp.Rational(13, 2),
        sp.Rational(-9, 4), sp.Rational(-9, 2), sp.Rational(-10),
        sp.Rational(-3, 2), sp.Rational(-3, 4), sp.Rational(-153, 8),
        sp.Rational(1), sp.Rational(13, 2),
    ),
]

images_11 = [evaluate_map(phi, variables_11, point) for point in points_11]
assert len(set(points_11)) == 3
assert images_11[0] == images_11[1] == images_11[2]
assert images_11[0] == (sp.Rational(-1, 4),) + (sp.Rational(0),) * 10
assert max(sp.Poly(component, *variables_11).total_degree() for component in phi) == 3


# ---------------------------------------------------------------------------
# 2. Normalize the linear part: F = L^{-1} o Phi = X + Q + C.
# ---------------------------------------------------------------------------

origin_11 = {variable: 0 for variable in variables_11}
linear_part = sp.Matrix(phi).jacobian(variables_11).subs(origin_11)
assert linear_part.det() == -2

normalized = [
    sp.expand(component)
    for component in list(linear_part.inv() * sp.Matrix(phi))
]
assert sp.Matrix(normalized).jacobian(variables_11).subs(origin_11) == sp.eye(11)

quadratic_parts: list[sp.Expr] = []
cubic_parts: list[sp.Expr] = []
for index, component in enumerate(normalized):
    nonlinear = sp.expand(component - variables_11[index])
    quadratic = homogeneous_part(nonlinear, variables_11, 2)
    cubic = homogeneous_part(nonlinear, variables_11, 3)
    assert sp.expand(component - variables_11[index] - quadratic - cubic) == 0
    quadratic_parts.append(quadratic)
    cubic_parts.append(cubic)

cubic_indices = [index for index, component in enumerate(cubic_parts) if component != 0]
assert cubic_indices == [0, 1, 2, 3, 4, 5, 8]

# The seven nonzero cubic components are linearly independent over Q.
cubic_monomials = sorted(
    set().union(*(set(sp.Poly(cubic_parts[i], *variables_11).monoms()) for i in cubic_indices))
)
cubic_coefficient_matrix = sp.Matrix([
    [sp.Poly(cubic_parts[i], *variables_11).coeff_monomial(m) for m in cubic_monomials]
    for i in cubic_indices
])
assert cubic_coefficient_matrix.rank() == 7


# ---------------------------------------------------------------------------
# 3. Sparse Segre/homogenization step.
#
# For S = {1,2,3,4,5,6,9}, introduce Y_i only when C_i is nonzero:
#
#   Ghat(X,Y,t) = (X + t Q(X) - t^2 iota(Y), Y + C_S(X), t).
#
# Dimension = 11 + 7 + 1 = 19.
# ---------------------------------------------------------------------------

Y1, Y2, Y3, Y4, Y5, Y6, Y9 = sp.symbols("Y1 Y2 Y3 Y4 Y5 Y6 Y9")
y_variables = [Y1, Y2, Y3, Y4, Y5, Y6, Y9]
t = sp.Symbol("t")
y_for_index = dict(zip(cubic_indices, y_variables))

final_map: list[sp.Expr] = []
for index, variable in enumerate(variables_11):
    component = variable + t * quadratic_parts[index]
    if index in y_for_index:
        component -= t**2 * y_for_index[index]
    final_map.append(sp.expand(component))

for index in cubic_indices:
    final_map.append(sp.expand(y_for_index[index] + cubic_parts[index]))
final_map.append(t)

variables_19 = variables_11 + y_variables + [t]
assert len(final_map) == len(variables_19) == 19

H = [sp.expand(component - variable) for component, variable in zip(final_map, variables_19)]
for index, component in enumerate(H):
    if component == 0:
        continue
    poly = sp.Poly(component, *variables_19)
    assert poly.is_homogeneous, index
    assert poly.total_degree() == 3, index


# ---------------------------------------------------------------------------
# 4. Lift and verify the three-point collision exactly.
# ---------------------------------------------------------------------------

points_19: list[Tuple[sp.Expr, ...]] = []
for point in points_11:
    substitution = dict(zip(variables_11, point))
    lifted_y = tuple(-sp.expand(cubic_parts[index].subs(substitution)) for index in cubic_indices)
    points_19.append(tuple(point) + lifted_y + (sp.Integer(1),))

images_19 = [evaluate_map(final_map, variables_19, point) for point in points_19]
assert len(set(points_19)) == 3
assert images_19[0] == images_19[1] == images_19[2]
common_image_19 = images_19[0]
assert common_image_19 == (
    sp.Integer(0), sp.Integer(0), sp.Rational(-1, 4),
    sp.Integer(0), sp.Integer(0), sp.Rational(1, 2),
    sp.Integer(0), sp.Integer(0), sp.Integer(0), sp.Integer(0),
    sp.Integer(0), sp.Integer(0), sp.Integer(0), sp.Integer(0),
    sp.Integer(0), sp.Integer(0), sp.Integer(0), sp.Integer(0),
    sp.Integer(1),
)


# ---------------------------------------------------------------------------
# 5. Compare the derived object with the canonical JSON statement.
# ---------------------------------------------------------------------------

json_path = Path(__file__).with_name("yagzhev19_map.json")
json_raw = json_path.read_bytes()
assert hashlib.sha256(json_raw).hexdigest() == EXPECTED_JSON_SHA256
canonical = json.loads(json_raw)
assert canonical["dimension"] == 19
assert canonical["variable_order"] == [str(variable) for variable in variables_19]

def json_component(terms: Sequence[Mapping[str, object]]) -> sp.Expr:
    expression = sp.Integer(0)
    for term in terms:
        coefficient = sp.Rational(str(term["coeff"]))
        monomial = coefficient
        for variable, exponent in zip(variables_19, term["exponents"]):
            monomial *= variable ** int(exponent)
        expression += monomial
    return sp.expand(expression)

assert [json_component(terms) for terms in canonical["H"]] == H
assert [tuple(sp.Rational(value) for value in row) for row in canonical["collision_points"]] == points_19
assert tuple(sp.Rational(value) for value in canonical["common_image"]) == common_image_19


# ---------------------------------------------------------------------------
# 6. Global exact nilpotency certificate for JH.
# ---------------------------------------------------------------------------

JH = sparse_jacobian(H, variables_19)
power = JH
power_17: PolyMatrix | None = None
for exponent in range(1, 19):
    if exponent == 17:
        power_17 = dict(power)
    if exponent < 18:
        power = sparse_multiply(power, JH)

assert power_17 is not None
assert power_17 != {}
assert power == {}  # (JH)^18 = 0 as a polynomial matrix.

# Paper matrix-entry indices are one-based; Python uses the zero-based key (3, 18).
explicit_nonzero_entry = sp.expand(power_17[(3, 18)])
assert explicit_nonzero_entry == -18 * t**16 * x**10 * y**6 * z**2

# Since JH is nilpotent, I + JH has the finite inverse
# I - JH + ... + (-1)^17 (JH)^17. Therefore det(I + JH) is a unit.
# Evaluating at the origin gives det(I + JH) = 1.


# ---------------------------------------------------------------------------
# 7. Generate and verify the LaTeX appendices from the checked expressions.
# ---------------------------------------------------------------------------

def latex_rational(value: sp.Expr) -> str:
    return sp.latex(sp.Rational(value))


def render_generated_appendices() -> str:
    lines: list[str] = []
    lines.extend([
        r"% GENERATED FILE -- DO NOT EDIT BY HAND.",
        r"% Generated and byte-for-byte checked by verify_yagzhev19_sympy.py.",
        r"\section{The 11-variable map}\label{app:phi}",
        r"With coordinate order $(x,y,z,a,b,c,d,q,s,h,k)$, the components are",
        r"\begin{align*}",
    ])
    for index, component in enumerate(phi, start=1):
        suffix = r",\\" if index < len(phi) else r"."
        lines.append(rf"\Phi_{{{index}}}={{}}&{sp.latex(component)}{suffix}")
    lines.extend([
        r"\end{align*}",
        r"",
        r"\section{Quadratic and cubic parts after normalization}\label{app:qc}",
        r"The nonzero components of $Q$ are",
        r"\begin{align*}",
    ])
    q_indices = [i for i, component in enumerate(quadratic_parts) if component != 0]
    for position, index in enumerate(q_indices):
        suffix = r",\\" if position < len(q_indices) - 1 else r"."
        lines.append(rf"Q_{{{index + 1}}}={{}}&{sp.latex(quadratic_parts[index])}{suffix}")
    lines.extend([
        r"\end{align*}",
        r"All omitted $Q_i$ are zero. The nonzero components of $C$ are",
        r"\begin{align*}",
    ])
    for position, index in enumerate(cubic_indices):
        suffix = r",\\" if position < len(cubic_indices) - 1 else r"."
        lines.append(rf"C_{{{index + 1}}}={{}}&{sp.latex(cubic_parts[index])}{suffix}")
    lines.extend([
        r"\end{align*}",
        r"",
        r"\section{The lifted collision points}\label{app:points}",
        r"The lifted points are $\widehat p_i=(p_i,-C_S(p_i),1)$. In the coordinate order",
        r"\[",
        r"(x,y,z,a,b,c,d,q,s,h,k,Y_1,Y_2,Y_3,Y_4,Y_5,Y_6,Y_9,t),",
        r"\]",
        r"they are",
        r"\begin{align*}",
    ])
    for point_index, point in enumerate(points_19, start=1):
        values = [latex_rational(value) for value in point]
        if point_index == 1:
            body = ",".join(values)
            suffix = r",\\"
            lines.append(rf"\widehat p_{{{point_index}}}={{}}&({body}){suffix}")
        else:
            first = ",".join(values[:11])
            second = ",".join(values[11:])
            suffix = r",\\" if point_index < len(points_19) else r"."
            lines.append(rf"\widehat p_{{{point_index}}}={{}}&({first},\\")
            lines.append(rf"&\hspace{{1cm}}{second}){suffix}")
    lines.extend([
        r"\end{align*}",
        r"Their common image is",
        r"\[",
        "(" + ",".join(latex_rational(value) for value in common_image_19) + ").",
        r"\]",
        r"",
    ])
    return "\n".join(lines)


generated_appendices = render_generated_appendices()
generated_path = Path(__file__).with_name("yagzhev19_appendices.tex")
if "--write-generated" in sys.argv:
    generated_path.write_text(generated_appendices, encoding="utf-8")
else:
    assert generated_path.read_text(encoding="utf-8") == generated_appendices


print("Verified explicit cubic-homogeneous Jacobian counterexample:")
print("  canonical JSON SHA-256:", EXPECTED_JSON_SHA256)
print("  dimension: 19")
print("  cubic indices after normalization: 1, 2, 3, 4, 5, 6, 9")
print("  every nonzero component of H is homogeneous of degree 3")
print("  (JH)^18 = 0")
print("  (JH)^17 != 0; for example")
print("    ((JH)^17)_[4,19] =", explicit_nonzero_entry)
print("  generated LaTeX appendices match the checked expressions byte-for-byte")
print("  three distinct rational inputs share the image:")
print("   ", common_image_19)
print("  hence det J(Ghat) = 1 identically and Ghat is non-injective")