#!/usr/bin/env python3
"""Exact certificate verifier for the frequency-22 Class-I cubature theorem.

This program uses only the Python standard library.  Every equality and every
sign decision is made in Q(sqrt(5)); floating-point numbers are used only for
the short decimal summaries printed at the end.

It verifies two statements for the radially projected Class-I icosahedral
grid X_22, with weights constant on full H3 orbits:

1. A strictly positive rule exists through algebraic degree 69.
2. No nonnegative rule exists through degree 70 (equivalently 71).

The second statement is certified twice: the unique degree-70/71 moment
solution has a negative mass on orbit (21,1,0), and an exact Farkas/cardinal
dual polynomial has the same negative spherical mean.
"""

from dataclasses import dataclass
from decimal import Decimal, localcontext
from fractions import Fraction
from hashlib import sha256
from itertools import product
from math import factorial
from time import perf_counter


@dataclass(frozen=True)
class K5:
    """An exact number a + b*sqrt(5), with a,b rational."""

    a: Fraction = Fraction(0)
    b: Fraction = Fraction(0)

    def __init__(self, a=0, b=0):
        object.__setattr__(self, "a", a if isinstance(a, Fraction) else Fraction(a))
        object.__setattr__(self, "b", b if isinstance(b, Fraction) else Fraction(b))

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

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

    __radd__ = __add__

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

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

    def __rsub__(self, other):
        return K5.coerce(other) - self

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

    __rmul__ = __mul__

    def inv(self):
        denominator = self.a * self.a - 5 * self.b * self.b
        if denominator == 0:
            raise ZeroDivisionError
        return K5(self.a / denominator, -self.b / denominator)

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

    def __rtruediv__(self, other):
        return K5.coerce(other) / self

    def __pow__(self, exponent):
        if exponent < 0:
            return self.inv() ** (-exponent)
        result, base = K5(1), self
        while exponent:
            if exponent & 1:
                result *= base
            base *= base
            exponent //= 2
        return result

    def sign(self):
        """Exact sign under the real embedding sqrt(5)>0."""
        if self.b == 0:
            return (self.a > 0) - (self.a < 0)
        if self.a >= 0 and self.b >= 0:
            return 1
        if self.a <= 0 and self.b <= 0:
            return -1
        aa, five_bb = self.a * self.a, 5 * self.b * self.b
        if aa == five_bb:
            return 0
        if self.a > 0:
            return 1 if aa > five_bb else -1
        return -1 if aa > five_bb else 1

    def decimal(self, precision=80):
        """High-precision orientation decimal; never used for a decision."""
        with localcontext() as context:
            context.prec = precision
            a = Decimal(self.a.numerator) / Decimal(self.a.denominator)
            b = Decimal(self.b.numerator) / Decimal(self.b.denominator)
            return a + b * Decimal(5).sqrt()


SQRT5 = K5(0, 1)
PHI = (K5(1) + SQRT5) / 2
R2 = PHI + 2

# Standard unnormalized icosahedron vertices; every norm squared is R2.
VERTICES = []
for sign_1, sign_2 in product((1, -1), repeat=2):
    VERTICES.extend(
        [
            (K5(0), K5(sign_1), K5(sign_2) * PHI),
            (K5(sign_1), K5(sign_2) * PHI, K5(0)),
            (K5(sign_1) * PHI, K5(0), K5(sign_2)),
        ]
    )

FACE = (
    (K5(0), K5(1), PHI),
    (K5(-1), PHI, K5(0)),
    (K5(1), PHI, K5(0)),
)


def dot(left, right):
    return sum((x * y for x, y in zip(left, right)), K5(0))


# Sparse homogeneous polynomials: exponent triple -> K5 coefficient.
def poly_add(left, right):
    result = dict(left)
    for exponent, coefficient in right.items():
        result[exponent] = result.get(exponent, K5(0)) + coefficient
        if result[exponent] == K5(0):
            del result[exponent]
    return result


def poly_scale(poly, scalar):
    scalar = K5.coerce(scalar)
    return {
        exponent: scalar * coefficient
        for exponent, coefficient in poly.items()
        if scalar * coefficient != K5(0)
    }


def poly_multiply(left, right):
    result = {}
    for (i, j, k), a in left.items():
        for (u, v, w), b in right.items():
            exponent = (i + u, j + v, k + w)
            result[exponent] = result.get(exponent, K5(0)) + a * b
    return {e: c for e, c in result.items() if c != K5(0)}


def multinomial_3(n, i, j, k):
    return factorial(n) // (factorial(i) * factorial(j) * factorial(k))


def linear_power(vector, n):
    result = {}
    vx, vy, vz = vector
    for i in range(n + 1):
        for j in range(n - i + 1):
            k = n - i - j
            coefficient = (
                K5(multinomial_3(n, i, j, k)) * vx**i * vy**j * vz**k
            )
            if coefficient != K5(0):
                result[(i, j, k)] = coefficient
    return result


def invariant_polynomial(n):
    result = {}
    for vertex in VERTICES:
        result = poly_add(result, linear_power(vertex, n))
    return poly_scale(result, R2 ** (-(n // 2)))


U6_POLY = invariant_polynomial(6)
U10_POLY = invariant_polynomial(10)


def odd_double_factorial(n):
    result = 1
    for value in range(n, 0, -2):
        result *= value
    return result


def sphere_monomial_mean(i, j, k):
    if i % 2 or j % 2 or k % 2:
        return Fraction(0)
    a, b, c = i // 2, j // 2, k // 2
    return Fraction(
        odd_double_factorial(2 * a - 1)
        * odd_double_factorial(2 * b - 1)
        * odd_double_factorial(2 * c - 1),
        odd_double_factorial(2 * (a + b + c) + 1),
    )


def sphere_integral(poly):
    return sum(
        (
            coefficient * K5(sphere_monomial_mean(i, j, k))
            for (i, j, k), coefficient in poly.items()
        ),
        K5(0),
    )


POWERS_6 = {0: {(0, 0, 0): K5(1)}, 1: U6_POLY}
POWERS_10 = {0: {(0, 0, 0): K5(1)}, 1: U10_POLY}


def invariant_moment(r, s):
    while max(POWERS_6) < r:
        n = max(POWERS_6) + 1
        POWERS_6[n] = poly_multiply(POWERS_6[n - 1], U6_POLY)
    while max(POWERS_10) < s:
        n = max(POWERS_10) + 1
        POWERS_10[n] = poly_multiply(POWERS_10[n - 1], U10_POLY)
    return sphere_integral(poly_multiply(POWERS_6[r], POWERS_10[s]))


def basis_pairs(degree):
    pairs = []
    for s in range(degree // 10 + 1):
        for r in range((degree - 10 * s) // 6 + 1):
            pairs.append((r, s))
    return sorted(pairs, key=lambda pair: (6 * pair[0] + 10 * pair[1], pair[1], pair[0]))


def partitions_3(total):
    result = []
    for a in range(total, -1, -1):
        for b in range(min(a, total - a), -1, -1):
            c = total - a - b
            if 0 <= c <= b:
                result.append((a, b, c))
    return result


def orbit_invariants(representative):
    a, b, c = representative
    point = tuple(
        K5(a) * FACE[0][i] + K5(b) * FACE[1][i] + K5(c) * FACE[2][i]
        for i in range(3)
    )
    norm_squared = dot(point, point)
    u6 = sum((dot(point, v) ** 6 for v in VERTICES), K5(0)) / (
        R2**3 * norm_squared**3
    )
    u10 = sum((dot(point, v) ** 10 for v in VERTICES), K5(0)) / (
        R2**5 * norm_squared**5
    )
    return u6, u10


def rref(matrix, right_hand_side):
    augmented = [row[:] + [rhs] for row, rhs in zip(matrix, right_hand_side)]
    row_count = len(augmented)
    column_count = len(augmented[0]) - 1
    pivot_row = 0
    pivots = []
    for column in range(column_count):
        nonzero = next(
            (row for row in range(pivot_row, row_count) if augmented[row][column] != K5(0)),
            None,
        )
        if nonzero is None:
            continue
        augmented[pivot_row], augmented[nonzero] = (
            augmented[nonzero],
            augmented[pivot_row],
        )
        inverse = augmented[pivot_row][column].inv()
        augmented[pivot_row] = [value * inverse for value in augmented[pivot_row]]
        for row in range(row_count):
            if row != pivot_row and augmented[row][column] != K5(0):
                factor = augmented[row][column]
                augmented[row] = [
                    augmented[row][j] - factor * augmented[pivot_row][j]
                    for j in range(column_count + 1)
                ]
        pivots.append(column)
        pivot_row += 1
        if pivot_row == row_count:
            break
    for row in range(pivot_row, row_count):
        if (
            all(augmented[row][column] == K5(0) for column in range(column_count))
            and augmented[row][column_count] != K5(0)
        ):
            raise AssertionError("inconsistent exact moment system")
    return augmented, pivots


def choose_solution(column_count, reduced, pivots, free_values):
    free_columns = [column for column in range(column_count) if column not in pivots]
    if len(free_columns) != len(free_values):
        raise AssertionError(("wrong number of free values", free_columns, free_values))
    solution = [K5(0)] * column_count
    for column, value in zip(free_columns, free_values):
        solution[column] = K5(value)
    for row, column in enumerate(pivots):
        value = reduced[row][-1]
        for free_column in free_columns:
            value -= reduced[row][free_column] * solution[free_column]
        solution[column] = value
    return solution, free_columns


def verify(matrix, right_hand_side, solution):
    return all(
        sum((row[j] * solution[j] for j in range(len(solution))), K5(0)) == rhs
        for row, rhs in zip(matrix, right_hand_side)
    )


def moment_system(frequency, degree):
    representatives = partitions_3(frequency)
    pairs = basis_pairs(degree)
    values = {rep: orbit_invariants(rep) for rep in representatives}
    matrix = [
        [(values[rep][0] ** r) * (values[rep][1] ** s) for rep in representatives]
        for r, s in pairs
    ]
    rhs = [invariant_moment(r, s) for r, s in pairs]
    reduced, pivots = rref(matrix, rhs)
    return representatives, pairs, values, matrix, rhs, reduced, pivots


def canonical_fraction(value):
    return f"{value.numerator}/{value.denominator}"


def canonical_k5(value):
    return f"{canonical_fraction(value.a)};{canonical_fraction(value.b)}"


def solution_digest(representatives, weights):
    payload = "\n".join(
        f"{a},{b},{c}:{canonical_k5(weight)}"
        for (a, b, c), weight in zip(representatives, weights)
    )
    return sha256(payload.encode("ascii")).hexdigest()


# These are orbit masses for the three free columns of the degree-69 system.
POSITIVE_FREE_VALUES = [
    Fraction(11607197, 394254842),
    Fraction(2312954, 161188803),
    Fraction(7156114, 472522843),
]

EXPECTED_DEGREE_69_SHA256 = (
    "0865572f76c7b8c9816ce937a79361c8926abe785ac3f41b06e4f289665bf340"
)
EXPECTED_DEGREE_71_SHA256 = (
    "c1e620fc9cf6d99623423cf82450540a5991679e8ded6632bd959914068949f2"
)


def main():
    started = perf_counter()
    frequency = 22

    reps, pairs69, values, matrix69, rhs69, reduced69, pivots69 = moment_system(
        frequency, 69
    )
    weights69, free_columns = choose_solution(
        len(reps), reduced69, pivots69, POSITIVE_FREE_VALUES
    )

    assert len(reps) == 52
    assert len(pairs69) == len(pivots69) == 49
    assert free_columns == [49, 50, 51]
    assert [reps[index] for index in free_columns] == [
        (9, 7, 6),
        (8, 8, 6),
        (8, 7, 7),
    ]
    assert verify(matrix69, rhs69, weights69)
    assert sum(weights69, K5(0)) == K5(1)
    assert all(weight.sign() > 0 for weight in weights69)

    # Select the minimum using exact arithmetic in Q(sqrt(5)).  Converting
    # a + b*sqrt(5) to binary64 before comparison is unsafe here because a
    # and b*sqrt(5) nearly cancel for some orbit masses.
    minimum_index = 0
    for index in range(1, len(weights69)):
        if (weights69[index] - weights69[minimum_index]).sign() < 0:
            minimum_index = index
    minimum69 = (
        weights69[minimum_index].decimal(),
        reps[minimum_index],
        weights69[minimum_index],
    )
    assert minimum69[1] == (11, 11, 0)
    digest69 = solution_digest(reps, weights69)
    assert digest69 == EXPECTED_DEGREE_69_SHA256

    # Degree 70 contributes exactly three new invariant monomials.  Degree 71
    # contributes none because every H3 invariant has even degree.
    pairs71 = basis_pairs(71)
    new_pairs = [pair for pair in pairs71 if pair not in pairs69]
    assert new_pairs == [(10, 1), (5, 4), (0, 7)]

    extra_matrix = [
        [(values[rep][0] ** r) * (values[rep][1] ** s) for rep in reps]
        for r, s in new_pairs
    ]
    extra_rhs = [invariant_moment(r, s) for r, s in new_pairs]

    # Express every weight as constant + sum_k coefficient[k]*y_k, where the
    # y_k are the three free degree-69 orbit masses.
    constants = [K5(0)] * len(reps)
    coefficients = [[K5(0)] * 3 for _ in reps]
    for k, column in enumerate(free_columns):
        coefficients[column][k] = K5(1)
    for row, column in enumerate(pivots69):
        constants[column] = reduced69[row][-1]
        for k, free_column in enumerate(free_columns):
            coefficients[column][k] = -reduced69[row][free_column]

    three_by_three = []
    three_rhs = []
    for row, rhs in zip(extra_matrix, extra_rhs):
        three_by_three.append(
            [
                sum((row[j] * coefficients[j][k] for j in range(len(reps))), K5(0))
                for k in range(3)
            ]
        )
        three_rhs.append(
            rhs - sum((row[j] * constants[j] for j in range(len(reps))), K5(0))
        )

    reduced3, pivots3 = rref(three_by_three, three_rhs)
    assert pivots3 == [0, 1, 2]
    forced_free, _ = choose_solution(3, reduced3, pivots3, [])

    weights71 = [
        constants[j]
        + sum((coefficients[j][k] * forced_free[k] for k in range(3)), K5(0))
        for j in range(len(reps))
    ]
    matrix71 = matrix69 + extra_matrix
    rhs71 = rhs69 + extra_rhs
    assert len(matrix71) == len(reps) == 52
    assert verify(matrix71, rhs71, weights71)
    assert sum(weights71, K5(0)) == K5(1)

    negative = [
        (rep, weight)
        for rep, weight in zip(reps, weights71)
        if weight.sign() < 0
    ]
    assert len(negative) == 1
    bad_rep, bad_weight = negative[0]
    assert bad_rep == (21, 1, 0)

    # Exact dual certificate: solve A^T c = e_bad.  The invariant polynomial
    # P=sum_i c_i B_i is 1 on the bad orbit and 0 on the other 51 orbits.
    # Its spherical mean c.b equals the forced negative orbit mass.
    bad_index = reps.index(bad_rep)
    transpose = [
        [matrix71[row][column] for row in range(len(matrix71))]
        for column in range(len(reps))
    ]
    selector = [K5(1) if j == bad_index else K5(0) for j in range(len(reps))]
    reduced_dual, pivots_dual = rref(transpose, selector)
    assert len(pivots_dual) == 52
    dual_coefficients, _ = choose_solution(52, reduced_dual, pivots_dual, [])
    assert verify(transpose, selector, dual_coefficients)
    dual_mean = sum(
        (dual_coefficients[i] * rhs71[i] for i in range(len(rhs71))), K5(0)
    )
    assert dual_mean == bad_weight
    assert dual_mean.sign() < 0
    digest71 = solution_digest(reps, weights71)
    assert digest71 == EXPECTED_DEGREE_71_SHA256

    print("OMEGA NAUGHT — EXACT CLASS-I f=22 CERTIFICATE")
    print("Arithmetic: Q(sqrt(5)); equality/sign checks: exact")
    print(f"nodes={10 * frequency * frequency + 2}  H3_orbits={len(reps)}")
    print(
        f"degree_69_conditions={len(pairs69)}  rank={len(pivots69)}  "
        f"affine_dimension={len(free_columns)}"
    )
    print(f"free_orbits={[reps[index] for index in free_columns]}")
    print(f"free_masses={POSITIVE_FREE_VALUES}")
    print(
        "degree_69: EXACT and STRICTLY POSITIVE; "
        f"minimum_orbit_mass={minimum69[0]:.16g} at {minimum69[1]}"
    )
    print(f"degree_69_solution_sha256={digest69}")
    print(f"new_degree_70_invariants={new_pairs}")
    print(
        "degree_70/71: UNIQUE exact solution; "
        f"negative_orbit_mass={bad_weight.decimal():.16g} at {bad_rep}"
    )
    print(f"degree_71_solution_sha256={digest71}")
    print("Farkas/cardinal dual: EXACT PASS; mean equals negative forced mass")
    print("THEOREM: maximum nonnegative H3-orbit-constant degree on X_22 is 69")
    print(f"CERTIFICATE PASS ({perf_counter() - started:.2f} seconds)")


if __name__ == "__main__":
    main()