#!/usr/bin/env python3
"""Check the hand part of the classical genus-five certificate exactly."""

from fractions import Fraction


def dot(x, y):
    return sum(a * b for a, b in zip(x, y))


def matvec(A, x):
    return tuple(dot(row, x) for row in A)


def add(x, y):
    return tuple(a + b for a, b in zip(x, y))


def determinant(A):
    """Exact Gaussian-elimination determinant; matrices here are at most 5x5."""
    A = [list(map(Fraction, row)) for row in A]
    answer = Fraction(1)
    for j in range(len(A)):
        pivot = next((i for i in range(j, len(A)) if A[i][j]), None)
        if pivot is None:
            return Fraction(0)
        if pivot != j:
            A[j], A[pivot] = A[pivot], A[j]
            answer = -answer
        q = A[j][j]
        answer *= q
        A[j] = [a / q for a in A[j]]
        for i in range(j + 1, len(A)):
            q = A[i][j]
            A[i] = [a - q * b for a, b in zip(A[i], A[j])]
    return answer


A_cl = (
    (4, 0, 0, 0, -2),
    (0, 4, 0, 0, -2),
    (0, 0, 4, 0, -2),
    (0, 0, 0, 4, -2),
    (-2, -2, -2, -2, 5),
)

C = (
    (-1, -1, -1, 0, 2),
    (-1, -1, 0, -1, 2),
    (-1, 0, -1, -1, 2),
    (0, -1, -1, -1, 2),
    (-1, -1, -1, -1, 2),
    (-1, -1, -1, -1, 3),
    (-2, -1, -1, -1, 3),
    (-1, -2, -1, -1, 3),
    (-1, -1, -2, -1, 3),
    (-1, -1, -1, -2, 3),
)
b = (1, 1, 1, 1, 1, 2, 2, 2, 2, 2)

y_a = (0, 0, -1, -1, 0)
y_b = (0, 1, 1, 1, 2)
p_a = (0, 0, 1, 1, 1)
p_b = (0, -1, -1, -1, -1)


def slack(y):
    return tuple(a - d for a, d in zip(matvec(C, y), b))


# Sylvester gives an independent positive-definiteness check.
assert tuple(determinant([row[:n] for row in A_cl[:n]]) for n in range(1, 6)) \
    == (4, 16, 64, 256, 256)

assert slack(y_a) == (0, 0, 1, 1, 1, 0, 0, 0, 1, 1)
assert slack(y_b) == (1, 1, 1, 0, 0, 1, 1, 0, 0, 0)

# Tight rows span the whole quotient lattice, so nothing nonzero in P can be
# subtracted from either floor.  Unimodularity is stronger than needed here.
tight_a = (0, 1, 5, 6, 7)
tight_b = (3, 4, 7, 8, 9)
assert determinant([C[i] for i in tight_a]) == -1
assert determinant([C[i] for i in tight_b]) == -1

assert matvec(C, p_a) == (1, 1, 0, 0, 0, 1, 1, 1, 0, 0)
assert matvec(C, p_b) == (0, 0, 0, 1, 1, 0, 0, 1, 1, 1)
assert p_a != (0,) * 5 and p_b != (0,) * 5
assert add(y_a, p_a) == add(y_b, p_b) == (0, 0, 0, 0, 1)

print("classical arithmetic certificate: verified")
print("  two tight-row determinants: -1, -1")
print("  common orbit point:", add(y_a, p_a))