#!/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.

This file is website revision r2 of the original certificate verifier
Omega_Naught_ClassI_f22_exact_verifier.py (16972 bytes, SHA-256
ae615d262ecfc564aed119e5002c4fa9c4fadf35812002426298caf7c6ab89bf).
The original remains the zip and SHA256SUMS inventory member. Revision r2
adds only a startup guard that refuses python3 -O / -OO / PYTHONOPTIMIZE
before any assertion. It does not change the certificate, the theorem, or
the canonical solution digests.
"""

if not __debug__:
    raise SystemExit(
        "Verification requires Python assertions. "
        "Run without -O or -OO and unset PYTHONOPTIMIZE."
    )

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