#!/usr/bin/env python3
"""Check the hand part of the rho=0 theta certificate exactly."""

from fractions import Fraction
from math import gcd


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):
    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


def rank(A):
    A = [list(map(Fraction, row)) for row in A]
    row = 0
    for col in range(len(A[0])):
        pivot = next((i for i in range(row, len(A)) if A[i][col]), None)
        if pivot is None:
            continue
        A[row], A[pivot] = A[pivot], A[row]
        q = A[row][col]
        A[row] = [a / q for a in A[row]]
        for i in range(len(A)):
            if i != row:
                q = A[i][col]
                A[i] = [a - q * b for a, b in zip(A[i], A[row])]
        row += 1
        if row == len(A):
            break
    return row


pairs = tuple((i, j) for i in range(5) for j in range(i, 5))
K_diag = tuple(2 if i == j else 1 for i, j in pairs)

B_theta = (
    (14, 9, 7, 7, 14),
    (9, 16, 7, 7, 14),
    (7, 7, 16, 7, 14),
    (7, 7, 7, 16, 14),
    (14, 14, 14, 14, 24),
)

vertices = (
    (1, 0, 0, 0, 0),
    (0, 1, 0, 0, 0),
    (0, 0, 1, 0, 0),
    (0, 0, 0, 1, 0),
    (1, 1, 1, 1, -2),
)
c = (1, 1, 1, 1, -1)

rays = (
    (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0),
    (0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0),
    (0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 1, 2, 1),
    (0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
    (1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 2, 1, 2, 2),
    (1, 2, 0, 2, 2, 1, 0, 2, 2, 0, 0, 0, 1, 2, 1),
    (1, 2, 2, 0, 2, 1, 2, 0, 2, 1, 0, 2, 0, 0, 1),
    (4, 4, 4, 4, 8, 4, 4, 4, 8, 4, 4, 8, 4, 8, 7),
)
beta = (1, 1, 1, 1, 3, 2, 2, 10)

f_a = (-4, 1, 0, 0, 0, 1, 0, 0, 0, 1, -1, 0, 1, 0, 2)
f_b = (-1, 0, 0, 0, 0, 1, 0, 0, 0, 1, -1, 0, 1, 0, 1)
p_a = (4, -1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -2)
p_b = (1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1)


def epsilon(v):
    """Coefficients of B(v,v)/2 when vector(B)=K_diag*z."""
    return tuple(v[i] * v[j] for i, j in pairs)


# The form lies in Q_0 and is positive definite by Sylvester's criterion.
assert all(B_theta[i][i] % 2 == 0 for i in range(5))
assert tuple(determinant([row[:n] for row in B_theta[:n]])
             for n in range(1, 6)) == (14, 143, 1700, 19017, 57744)

# The Delaunay simplex has index two and c is the half-sum of its generators.
assert determinant(list(zip(*vertices))) == -2
assert tuple(sum(v[i] for v in vertices) for i in range(5)) \
    == tuple(2 * x for x in c)

# The rays are primitive of rank eight and their positive sum is B_theta.
assert rank(rays) == 8
assert all(gcd(*map(abs, ray)) == 1 for ray in rays)
ray_sum = tuple(sum(ray[j] for ray in rays) for j in range(15))
assert tuple(k * x for k, x in zip(K_diag, ray_sum)) \
    == tuple(B_theta[i][j] for i, j in pairs)

# Compute the interpolation character and its eight ray values from scratch.
u = tuple(sum(Fraction(epsilon(v)[j], 2) for v in vertices)
          for j in range(15))
assert u == (1, Fraction(1, 2), Fraction(1, 2), Fraction(1, 2), -1,
             1, Fraction(1, 2), Fraction(1, 2), -1, 1,
             Fraction(1, 2), -1, 1, -1, 2)
assert matvec(rays, u) == beta

# The parity identity is the quantified minimality argument in one line.
parity_indices = (0, 1, 2, 7)
even_sum = tuple(sum(rays[i][j] for i in parity_indices) for j in range(15))
assert all(x % 2 == 0 for x in even_sum)
assert sum(beta[i] for i in parity_indices) == 13

values_a, values_b = matvec(rays, f_a), matvec(rays, f_b)
slack_a = tuple(x - y for x, y in zip(values_a, beta))
slack_b = tuple(x - y for x, y in zip(values_b, beta))
assert slack_a == (0, 0, 1, 0, 0, 0, 0, 0)
assert slack_b == (0, 0, 0, 0, 0, 0, 0, 1)
assert tuple(min(a, b) for a, b in zip(values_a, values_b)) == beta

assert matvec(rays, p_a) == (0, 0, 0, 0, 0, 0, 0, 2)
assert matvec(rays, p_b) == (0, 0, 1, 0, 0, 0, 0, 1)
overlap = (0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0)
assert add(f_a, p_a) == add(f_b, p_b) == overlap

print("theta arithmetic certificate: verified")
print("  even ray sum / odd bound sum:", even_sum, "/ 13")
print("  common orbit point:", overlap)