from fractions import Fraction
from math import factorial


class PhiField:
    """Exact a + b*phi arithmetic, where phi**2 = phi + 1."""

    def __init__(self, a=0, b=0):
        self.a = Fraction(a)
        self.b = Fraction(b)

    @staticmethod
    def coerce(value):
        return value if isinstance(value, PhiField) else PhiField(value)

    def __add__(self, other):
        other = self.coerce(other)
        return PhiField(self.a + other.a, self.b + other.b)

    __radd__ = __add__

    def __neg__(self):
        return PhiField(-self.a, -self.b)

    def __sub__(self, other):
        return self + (-self.coerce(other))

    def __mul__(self, other):
        other = self.coerce(other)
        return PhiField(
            self.a * other.a + self.b * other.b,
            self.a * other.b + self.b * other.a + self.b * other.b,
        )

    __rmul__ = __mul__

    def inverse(self):
        norm = self.a**2 + self.a * self.b - self.b**2
        if norm == 0:
            raise ZeroDivisionError("Zero field element")
        return PhiField((self.a + self.b) / norm, -self.b / norm)

    def __truediv__(self, other):
        return self * self.coerce(other).inverse()

    def __pow__(self, exponent):
        if not isinstance(exponent, int) or exponent < 0:
            raise ValueError("Expected a nonnegative integer exponent")
        result = PhiField(1)
        for _ in range(exponent):
            result = result * self
        return result

    def __eq__(self, other):
        other = self.coerce(other)
        return (self.a, self.b) == (other.a, other.b)


def require(condition, label):
    if not condition:
        raise RuntimeError(label)


def inverse_and_determinant(matrix):
    n = len(matrix)
    augmented = [
        row[:] + [PhiField(int(i == j)) for j in range(n)]
        for i, row in enumerate(matrix)
    ]
    determinant = PhiField(1)
    for j in range(n):
        pivot_row = next(
            (i for i in range(j, n) if augmented[i][j] != 0), None
        )
        if pivot_row is None:
            return PhiField(0), None
        if pivot_row != j:
            augmented[j], augmented[pivot_row] = (
                augmented[pivot_row], augmented[j]
            )
            determinant = -determinant
        pivot = augmented[j][j]
        determinant = determinant * pivot
        augmented[j] = [value / pivot for value in augmented[j]]
        for i in range(n):
            if i != j:
                multiplier = augmented[i][j]
                augmented[i] = [
                    value - multiplier * basis
                    for value, basis in zip(augmented[i], augmented[j])
                ]
    return determinant, [row[n:] for row in augmented]


def transpose(matrix):
    return [list(column) for column in zip(*matrix)]


def matrix_product(left, right):
    return [
        [sum((a * b for a, b in zip(row, column)), PhiField())
         for column in transpose(right)]
        for row in left
    ]


def matrix_vector(matrix, vector):
    return [
        sum((a * b for a, b in zip(row, vector)), PhiField())
        for row in matrix
    ]


def pairing(left, right):
    return sum((a * b for a, b in zip(left, right)), PhiField())


def run_checks():
    phi = PhiField(0, 1)
    identity = [[PhiField(int(i == j)) for j in range(6)]
                for i in range(6)]
    for seed in (Fraction(1, 7), Fraction(1, 2), Fraction(1, 648000)):
        s = PhiField(seed)
        B = s + 1
        A = s + phi**2
        nodes = [-A, -B, -s, s, B, A]
        V = [[node**j for j in range(6)] for node in nodes]
        determinant, inverse = inverse_and_determinant(V)
        expected = (
            8 * A * B * s * (A*A - B*B)**2
            * (A*A - s*s)**2 * (B*B - s*s)**2
        )
        require(determinant == expected, "Determinant factorization")
        require(matrix_product(V, inverse) == identity, "Right inverse")
        require(matrix_product(inverse, V) == identity, "Left inverse")
        require(
            matrix_product(transpose(V), transpose(inverse)) == identity,
            "Symplectic block identity",
        )
        beta = [PhiField(j - 2, Fraction(j, 11)) for j in range(6)]
        h = [PhiField(1 - j, Fraction(2*j, 9)) for j in range(6)]
        y = matrix_vector(V, beta)
        eta = matrix_vector(transpose(inverse), h)
        require(matrix_vector(inverse, y) == beta, "Reconstruction")
        require(pairing(eta, y) == pairing(h, beta), "Dual pairing")
        volume_product = (
            PhiField(64) / determinant
            * (64 * determinant / PhiField(factorial(6)))
        )
        require(
            volume_product == PhiField(Fraction(256, 45)),
            "Volume product",
        )
        print(f"s={seed}: exact algebra checks PASS")

    A, B, s = phi**2, PhiField(1), PhiField()
    nodes = [-A, -B, -s, s, B, A]
    determinant, inverse = inverse_and_determinant(
        [[node**j for j in range(6)] for node in nodes]
    )
    require(determinant == 0 and inverse is None, "Zero-seed boundary")
    print("s=0: excluded singular boundary PASS")


if __name__ == "__main__":
    run_checks()