#!/usr/bin/env python3
"""This-archive checks for the 21 September 2026 boundary-reconstruction note.

The nine-page PDF is the written proof. This script corroborates selected
identities and the later exact remarks on Ψ², nine-node ranks, and the
surviving kernel after q_*. It is not a proof-assistant certificate.

Revision 2026-09-22.r2: the J-squared spherical moment is computed exactly
with fractions.Fraction; the other checks remain numerical NumPy
corroboration. This is not an exact certificate for the whole theorem.

Requires: Python 3 with numpy.
"""

from __future__ import annotations

from fractions import Fraction
from itertools import combinations, product
from math import log, prod, sqrt

import numpy as np

PASSES: list[str] = []


def require(cond: bool, name: str) -> None:
    if not cond:
        raise SystemExit(f"FAIL\t{name}")
    PASSES.append(name)
    print(f"PASS\t{name}")


phi = (1.0 + sqrt(5.0)) / 2.0
phip = 1.0 - phi
tau = 1.0 / (phi + np.pi)
r = sqrt(1.0 + phi * phi)
rp = sqrt(1.0 + phi * phi)
rm = sqrt(1.0 + phip * phip)

M = np.array([[0.0, 1.0], [1.0, 1.0]])
K = (2.0 * M - np.eye(2)) / sqrt(5.0)
A = np.array([[0.0, tau, 0.0], [tau, tau, 0.0], [0.0, 0.0, -(tau**-2)]])
Ainv = np.linalg.inv(A)
D = A @ A.T
Dinv = np.linalg.inv(D)
R = np.array([[1.0 / rp, phi / rp, 0.0], [1.0 / rm, phip / rm, 0.0], [0.0, 0.0, 1.0]])


def Nxy(x: float, y: float) -> float:
    return x * x + x * y - y * y


def J(P) -> float:
    return float(P[2] * Nxy(P[0], P[1]))


def basis9(P) -> np.ndarray:
    x, y, z = P
    return np.array(
        [1.0, x, y, z, x * x - y * y, 2 * x * y, x * z, y * z, 2 * z * z - x * x - y * y]
    )


def basis10(P) -> np.ndarray:
    return np.concatenate([basis9(P), [J(P)]])


def KA(y, q) -> float:
    s = float(y @ Dinv @ q)
    yy = float(y @ Dinv @ y)
    return 1.0 + 3.0 * s + 2.5 * (3.0 * s * s - yy)


raw = []
for s1, s2 in product((-1.0, 1.0), repeat=2):
    raw.append(np.array([0.0, s1, s2 * phi]))
    raw.append(np.array([s1, s2 * phi, 0.0]))
    raw.append(np.array([s2 * phi, 0.0, s1]))
v_nodes = np.array([u / r for u in raw])
n_nodes = np.array([R.T @ v for v in v_nodes])
q_nodes = np.array([A @ n for n in n_nodes])
pstar = R.T @ (np.array([1.0, 1.0, 1.0]) / sqrt(3.0))
qstar = A @ pstar

pairs = [(0, 9), (1, 10), (2, 11), (3, 6), (4, 7), (5, 8)]


def pair_id(i: int) -> int:
    for k, pair in enumerate(pairs):
        if i in pair:
            return k
    raise KeyError(i)


def Ej(j: int, x: np.ndarray) -> float:
    aj = n_nodes[pairs[j][0]]
    return 1.0 / 6.0 + 5.0 / 4.0 * ((float(x @ aj) ** 2) - float(x @ x) / 3.0)


rng = np.random.default_rng(0)

# 1. Volume and J
require(abs(np.linalg.det(A) - 1.0) < 1e-12, "det_A_is_one")
maxJ = 0.0
for _ in range(40):
    xi = rng.normal(size=3)
    maxJ = max(maxJ, abs(J(A @ xi) - J(xi)))
require(maxJ < 1e-12, "J_invariant_under_A")
Ainv_claimed = np.array([[-1.0 / tau, 1.0 / tau, 0.0], [1.0 / tau, 0.0, 0.0], [0.0, 0.0, -(tau**2)]])
require(np.max(np.abs(Ainv - Ainv_claimed)) < 1e-12, "A_inverse")

# 2. Product form
require(np.max(np.abs(R @ R.T - np.eye(3))) < 1e-12, "R_orthogonal")
max_prod = 0.0
for _ in range(40):
    x, y, z = rng.normal(size=3)
    XYZ = R @ np.array([x, y, z])
    max_prod = max(max_prod, abs(J([x, y, z]) - sqrt(5.0) * XYZ[0] * XYZ[1] * XYZ[2]))
require(max_prod < 1e-12, "J_equals_sqrt5_XYZ")

# 3. Two-step identity
M2 = M @ M
combo = 1.5 * np.eye(2) + (sqrt(5.0) / 2.0) * K
require(np.max(np.abs(M2 - combo)) < 1e-15, "M2_equals_polynomial_in_K")
wK, VK = np.linalg.eig(K)
expK = VK @ np.diag(np.exp(2.0 * log(phi) * wK)) @ np.linalg.inv(VK)
require(np.max(np.abs(M2 - expK.real)) < 1e-12, "M2_equals_exp_2logphi_K")
max_psi = 0.0
c = tau**2
for _ in range(40):
    w = rng.normal(size=2)
    z = float(rng.normal())
    w2 = tau * (M @ (tau * (M @ w)))
    z2 = (-(tau**-2)) * (-(tau**-2) * z)
    ew = (c * (expK.real @ w), c ** (-2) * z)
    max_psi = max(max_psi, float(np.linalg.norm(w2 - ew[0])) + abs(z2 - ew[1]))
require(max_psi < 1e-12, "Psi_squared_equals_Phi")
require(abs(np.linalg.det(M) + 1.0) < 1e-15, "Psi_not_in_Phi_family")

# 4. Aligned vanishing
require(max(abs(np.dot(n, n) - 1.0) for n in n_nodes) < 1e-12, "aligned_nodes_on_sphere")
require(max(abs(J(n)) for n in n_nodes) < 1e-12, "J_vanishes_on_aligned_nodes")
require(max(abs(J(q)) for q in q_nodes) < 1e-12, "J_vanishes_on_transported_nodes")
claimed_J = sqrt(5.0) / (3.0 * sqrt(3.0))
require(abs(J(pstar) - claimed_J) < 1e-12, "J_pstar")
require(abs(J(qstar) - claimed_J) < 1e-12, "J_qstar")

# 5. E_j interpolation
interp_ok = True
for j in range(6):
    for i in range(6):
        for sign in (0, 1):
            val = Ej(j, n_nodes[pairs[i][sign]])
            expect = 1.0 if i == j else 0.0
            if abs(val - expect) > 1e-10:
                interp_ok = False
require(interp_ok, "Ej_interpolates_antipodal_pairs")

# 6. Gram determinants of three axes
s5 = sqrt(5.0)
claimed = {2.0 / 5.0 + 2.0 / 25.0 * s5, 2.0 / 5.0 - 2.0 / 25.0 * s5}
seen = set()
gram_ok = True
axes = [n_nodes[p[0]] for p in pairs]
for ijk in combinations(range(6), 3):
    for signs in product((-1.0, 1.0), repeat=3):
        vecs = [signs[t] * axes[ax] for t, ax in enumerate(ijk)]
        G = np.array([[float(vecs[p] @ vecs[q]) for q in range(3)] for p in range(3)])
        det = float(np.linalg.det(G))
        if det <= 1e-12:
            gram_ok = False
        matched = any(abs(det - c) < 1e-10 for c in claimed)
        if not matched:
            gram_ok = False
        seen.add(round(det, 12))
require(gram_ok and len(seen) == 2, "three_axes_Gram_positive")

# 7–8. Nine-subset ranks
hist9 = {8: 0, 9: 0}
hist10 = {(True, 9, 10): 0, (False, 8, 9): 0}
iff_ok = True
for comb in combinations(range(12), 9):
    represented = len({pair_id(i) for i in comb})
    M9 = np.stack([basis9(n_nodes[i]) for i in comb])
    M10 = np.vstack([np.stack([basis10(n_nodes[i]) for i in comb]), basis10(pstar)])
    r9 = int(np.linalg.matrix_rank(M9, tol=1e-8))
    r10 = int(np.linalg.matrix_rank(M10, tol=1e-8))
    hist9[r9] = hist9.get(r9, 0) + 1
    key = (represented == 6, r9, r10)
    hist10[key] = hist10.get(key, 0) + 1
    if (represented == 6) != (r9 == 9):
        iff_ok = False
    if represented < 6 and r9 != 8:
        iff_ok = False
require(iff_ok, "nine_nodes_iff_every_antipodal_pair")
require(hist9.get(9, 0) == 160 and hist9.get(8, 0) == 60, "counts_160_and_60")
require(
    hist10.get((True, 9, 10), 0) == 160 and hist10.get((False, 8, 9), 0) == 60,
    "adding_qstar_ranks_10_or_9",
)

# 9–10. Surviving kernel and E_j(p*)
plus = (1.0 + s5) / 6.0
minus = (1.0 - s5) / 6.0
ej_vals = [Ej(j, pstar) for j in range(6)]
coeff_ok = all(min(abs(v - plus), abs(v - minus)) < 1e-10 for v in ej_vals)
require(coeff_ok and all(abs(v) > 1e-10 for v in ej_vals), "Ej_pstar_is_one_pm_sqrt5_over_6")

# One failing selection: omit pair 0 and one extra
omit = set(pairs[0]) | {pairs[1][0]}
comb = [i for i in range(12) if i not in omit]
c_mix = Ej(0, pstar) / J(pstar)
max_h = 0.0
for i in comb:
    h = Ej(0, n_nodes[i]) - c_mix * J(n_nodes[i])
    max_h = max(max_h, abs(h))
h_star = Ej(0, pstar) - c_mix * J(pstar)
# quadratic part remains E_0, evaluated away from its pair
probe = n_nodes[pairs[0][0]]
h_probe = Ej(0, probe) - c_mix * J(probe)
require(max_h < 1e-10 and abs(h_star) < 1e-10, "surviving_kernel_vanishes_on_I_and_pstar")
require(abs(h_probe - 1.0) < 1e-10, "surviving_kernel_quadratic_part_nonzero")

# Transport: same mixture in y-coordinates
c_q = Ej(0, pstar) / J(qstar)
max_ht = 0.0
for i in comb:
    y = q_nodes[i]
    ht = Ej(0, Ainv @ y) - c_q * J(y)
    max_ht = max(max_ht, abs(ht))
ht_star = Ej(0, Ainv @ qstar) - c_q * J(qstar)
require(max_ht < 1e-9 and abs(ht_star) < 1e-9, "surviving_kernel_after_transport")

# 11. Twelve-sample quadratic decoder on H_A


def u_quad(X) -> float:
    return float(X[0] ** 2 - X[1] ** 2 + 3.0 * X[0] + 2.0)


def v_field(Y) -> float:
    return u_quad(Ainv @ Y)


max12 = 0.0
for _ in range(30):
    y = rng.normal(size=3)
    rec = sum(v_field(q) * KA(y, q) for q in q_nodes) / 12.0
    max12 = max(max12, abs(rec - v_field(y)))
require(max12 < 1e-9, "twelve_sample_quadratic_decoder")

# 12. Thirteen-sample decoder on F_A


def f_field(P) -> float:
    return v_field(P) + 4.0 * J(P)


max13 = 0.0
max12f = 0.0
for _ in range(30):
    y = rng.normal(size=3)
    vhat = sum(f_field(q) * KA(y, q) for q in q_nodes) / 12.0
    vhat_s = sum(f_field(q) * KA(qstar, q) for q in q_nodes) / 12.0
    rec = vhat + J(y) / J(qstar) * (f_field(qstar) - vhat_s)
    max13 = max(max13, abs(rec - f_field(y)))
    max12f = max(max12f, abs(vhat - f_field(y)))
require(max13 < 1e-9, "thirteen_sample_decoder_on_FA")
require(max12f > 0.1, "twelve_samples_do_not_determine_J")

# 13. Degree-six obstruction
require(max(J(n) ** 2 for n in n_nodes) < 1e-24, "J2_zero_on_twelve_nodes")
# Exact normalized spherical moments on S^2. Odd exponents integrate to zero.
# For even exponents e=(2a,2b,2c), the mean is
# (2a-1)!! (2b-1)!! (2c-1)!! / (2a+2b+2c+1)!!, with (-1)!!=1.
def spherical_monomial_mean(exponents: tuple[int, int, int]) -> Fraction:
    if len(exponents) != 3 or any(type(e) is not int or e < 0 for e in exponents):
        raise ValueError("Expected three nonnegative integer exponents")
    if any(e % 2 for e in exponents):
        return Fraction(0)
    numerator = prod(prod(range(1, e, 2)) for e in exponents)
    denominator = prod(range(1, sum(exponents) + 2, 2))
    return Fraction(numerator, denominator)


def multiply_polynomials(left, right):
    result = {}
    for a, ca in left.items():
        for b, cb in right.items():
            powers = tuple(x + y for x, y in zip(a, b))
            result[powers] = result.get(powers, Fraction(0)) + ca * cb
    return {powers: coefficient for powers, coefficient in result.items() if coefficient}


# J(x,y,z) = x^2 z + xyz - y^2 z. Square the actual coefficient record.
J_POLYNOMIAL = {
    (2, 0, 1): Fraction(1),
    (1, 1, 1): Fraction(1),
    (0, 2, 1): Fraction(-1),
}
J2_POLYNOMIAL = multiply_polynomials(J_POLYNOMIAL, J_POLYNOMIAL)
J2_SPHERICAL_MEAN = sum(
    (coefficient * spherical_monomial_mean(powers)
     for powers, coefficient in J2_POLYNOMIAL.items()),
    Fraction(0),
)
require(J2_SPHERICAL_MEAN == Fraction(1, 21), "J2_spherical_mean_one_twenty_first")

# 14. Surface Jacobians
e1 = np.array([1.0, 0.0, 0.0])
e2 = np.array([0.0, 1.0, 0.0])
ATinv = np.linalg.inv(A).T
require(abs(np.linalg.norm(ATinv @ e1) - sqrt(2.0) / tau) < 1e-12, "jacobian_e1")
require(abs(np.linalg.norm(ATinv @ e2) - 1.0 / tau) < 1e-12, "jacobian_e2")

print(f"TOTAL\t{len(PASSES)}\tPASS")