"""Exact finite-field rank certificate for the Tiles determinant cubic.

All input coefficients and matrix entries are integers.  Full rank modulo a
prime proves full rank over Q; hence the degree-five part of the Jacobian ideal
is the whole degree-five polynomial space.  This proves that the projective
cubic is smooth over C.  It is an exact certificate, not a numerical test.
"""

from itertools import permutations, product
from math import isqrt

PRIME = 1000003
assert all(PRIME % d for d in range(2, isqrt(PRIME) + 1))

# Variables are (u,v,w,z); S=[[u,u,v],[w,-2(u+v+w+z),v],[w,z,z]].
F = {
    (2, 0, 0, 1): -2,
    (1, 1, 0, 1): -3,
    (1, 0, 1, 1): -3,
    (1, 0, 0, 2): -2,
    (1, 1, 1, 0): 3,
    (0, 2, 1, 0): 2,
    (0, 1, 2, 0): 2,
    (0, 1, 1, 1): 3,
}


def add(left, right, factor=1):
    result = dict(left)
    for exponent, coefficient in right.items():
        result[exponent] = result.get(exponent, 0) + factor * coefficient
    return {e: c for e, c in result.items() if c}


def multiply(left, right):
    result = {}
    for e, c in left.items():
        for f, d in right.items():
            exponent = tuple(a + b for a, b in zip(e, f))
            result[exponent] = result.get(exponent, 0) + c * d
    return {e: c for e, c in result.items() if c}


variables = [{tuple(int(i == j) for i in range(4)): 1} for j in range(4)]
u, v, w, z = variables
middle = {}
for variable in variables:
    middle = add(middle, variable, -2)
seed = [[u, u, v], [w, middle, v], [w, z, z]]
determinant = {}
for permutation in permutations(range(3)):
    term = {(0, 0, 0, 0): 1}
    for i, j in enumerate(permutation):
        term = multiply(term, seed[i][j])
    inversions = sum(permutation[i] > permutation[j] for i in range(3) for j in range(i + 1, 3))
    determinant = add(determinant, term, (-1) ** inversions)
assert determinant == F


def monomials(degree):
    return [e for e in product(range(degree + 1), repeat=4) if sum(e) == degree]


def derivative(poly, variable):
    answer = {}
    for exponent, coefficient in poly.items():
        if exponent[variable]:
            new_exponent = list(exponent)
            new_exponent[variable] -= 1
            answer[tuple(new_exponent)] = coefficient * exponent[variable]
    return answer


rows = monomials(5)
row_index = {e: i for i, e in enumerate(rows)}
columns = []
for variable in range(4):
    partial = derivative(F, variable)
    for monomial in monomials(3):
        column = [0] * len(rows)
        for exponent, coefficient in partial.items():
            output = tuple(a + b for a, b in zip(exponent, monomial))
            column[row_index[output]] = coefficient
        columns.append(column)

matrix = [[column[i] % PRIME for column in columns] for i in range(len(rows))]
pivots = []
pivot_product = 1
sign = 1
rank = 0
for column in range(len(columns)):
    pivot = next((i for i in range(rank, len(rows)) if matrix[i][column]), None)
    if pivot is None:
        continue
    if pivot != rank:
        matrix[pivot], matrix[rank] = matrix[rank], matrix[pivot]
        sign = -sign
    value = matrix[rank][column]
    pivot_product = pivot_product * value % PRIME
    inverse = pow(value, -1, PRIME)
    matrix[rank] = [x * inverse % PRIME for x in matrix[rank]]
    for i in range(rank + 1, len(rows)):
        factor = matrix[i][column]
        if factor:
            matrix[i] = [(x - factor * y) % PRIME for x, y in zip(matrix[i], matrix[rank])]
    pivots.append(column)
    rank += 1
    if rank == len(rows):
        break

assert rank == 56
assert sign * pivot_product % PRIME == 152546
print("prime =", PRIME)
print("integer polynomial determinant identity = verified")
print("Jacobian degree-five matrix shape =", (len(rows), len(columns)))
print("rank =", rank)
print("pivot columns (zero based) =", pivots)
print("selected 56x56 minor modulo prime =", sign * pivot_product % PRIME)
print("Conclusion: all degree-five monomials lie in the rational Jacobian ideal.")
print("Conclusion: the projective Tiles determinant cubic is smooth over C.")