#!/usr/bin/env python3
"""verify_P306.py -- Verifier for Addendum 306 (P18-T2 sub-result (ii),
additive-composition uniqueness).

Asserts, recomputing from scratch, that ADDITIVE composition (P18 S4.1,
O_hat = sum of sector operators) is the unique way to compose the named
sectors GIVEN the independent-decoupling / no-forced-cross-coupling
axiom. Two sectors of P18's radial operator are used (declared): the
density coupling alpha*rho(r) [R5] and the renormalization generator
zeta*(r d/dr) [R7], over the fixed (kinetic + V_self) operator K, at the
canonical couplings alpha = 1/137.036, zeta = alpha^{1.25}.

The discriminator is the OPERATOR mixed second partial d^2 O / da dz, the
forced cross-coupling. Additive O is linear in (a, z) so this vanishes as
an exact operator identity; the multiplicative alternative
O = K (I + a Rho_n)(I + z Ren_n) carries the cross-term K Rho_n Ren_n.
The decoupling limit confirms it at the spectral level: setting one
coupling to zero recovers the bare reduced operator for additive, but
the multiplicative form distorts the surviving sector.

Tolerances: tol_zero = 1e-4 (an operator/spectrum residual below this is
treated as vanishing), tol_nonzero = 1e-2 (a cross-coupling above this is
treated as forced). The additive operator mixed partial is 0 to machine
precision; the multiplicative one matches ||K Rho_n Ren_n||_F to 1e-6
relative. This is a partial on P18-T2 sub-result (ii) ONLY: it does not
rule out operadic compositions engineered to decouple, does not prove the
decoupling axiom from deeper principles, and A294's re-typing of the
whole P18-T2 stands.

  S1  Sector pieces and the two compositions   - checks 1-2
  S2  T1: forced cross-coupling at the operator - checks 3-5
  S3  T2: clean decoupling limit                - checks 6-8
"""
import sys
import math

import numpy as np

PASS = FAIL = 0


def check(n, desc, cond):
    global PASS, FAIL
    ok = bool(cond)
    PASS += ok
    FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")


PI = math.pi
N = 500
NMODES = 6
A0 = 1.0 / 137.036
Z0 = A0 ** 1.25
ESELF = 13.177
M0 = 137.036
R0 = A0 ** 1.25
TOL_ZERO = 1e-4
TOL_NONZERO = 1e-2


def rho(r):
    return 16 * PI**3 * r**3 + 3 * PI**2 * r**2 + 2 * PI * r


def build_pieces(l=0):
    dr = 1.0 / N
    r = np.linspace(dr, 1 - dr, N - 1)
    n = len(r)
    d2 = (np.diag(-2.0 * np.ones(n)) + np.diag(np.ones(n - 1), 1)
          + np.diag(np.ones(n - 1), -1)) / dr**2
    d1 = (np.diag(np.ones(n - 1), 1) - np.diag(np.ones(n - 1), -1)) / (2 * dr)
    Vself = (ESELF / M0**2) * (1 - np.exp(-r / R0))
    K = (-d2 - np.diag(3.0 / r) @ d1 + np.diag(l * (l + 2) / r**2)
         + np.diag(Vself))
    return K, np.diag(rho(r)), np.diag(r) @ d1


def lowest_eigs(M, m=NMODES):
    w = np.linalg.eigvals(M)
    return np.sort(w.real)[:m]


K, Rho, Ren = build_pieces(l=0)
Rho_n = Rho / float(np.linalg.norm(Rho, 2))
Ren_n = Ren / float(np.linalg.norm(Ren, 2))
EYE = np.eye(K.shape[0])


def O_add(a, z):
    return K + a * Rho_n + z * Ren_n


def O_mult(a, z):
    return K @ (EYE + a * Rho_n) @ (EYE + z * Ren_n)


def lam0(builder, a, z):
    return float(lowest_eigs(builder(a, z), 1)[0])


h = 1e-3


def mixed_d2_operator(builder):
    return (builder(+h, +h) - builder(+h, -h) - builder(-h, +h)
            + builder(-h, -h)) / (4 * h * h)


d2O_add = float(np.linalg.norm(mixed_d2_operator(O_add), "fro"))
d2O_mult = float(np.linalg.norm(mixed_d2_operator(O_mult), "fro"))
cross_analytic = float(np.linalg.norm(K @ Rho_n @ Ren_n, "fro"))

# Hellmann-Feynman: additive per-sector response is invariant under z.
def d_da(builder, z):
    return (lam0(builder, +h, z) - lam0(builder, -h, z)) / (2 * h)
add_inv = abs(d_da(O_add, Z0) - d_da(O_add, 0.0))
mult_drift = abs(d_da(O_mult, Z0) - d_da(O_mult, 0.0))

# Decoupling limit at canonical alpha.
spec_add_z0 = lowest_eigs(O_add(A0, 0.0))
spec_reduced = lowest_eigs(K + A0 * Rho_n)
spec_mult_z0 = lowest_eigs(O_mult(A0, 0.0))
add_resid = float(np.max(np.abs(spec_add_z0 - spec_reduced)))
mult_vs_add = float(np.max(np.abs(spec_mult_z0 - spec_add_z0)))

print("S1  Sector pieces and the two compositions")
check(1, "two distinct sectors built: density Rho [R5], renorm Ren [R7] "
      "(||Rho||,||Ren|| > 0)",
      np.linalg.norm(Rho) > 0 and np.linalg.norm(Ren) > 0
      and not np.allclose(Rho, Ren))
check(2, "additive and multiplicative compositions agree to first order "
      "but differ (||O_add - O_mult|| > 0 at canonical)",
      not np.allclose(O_add(A0, Z0), O_mult(A0, Z0)))

print("S2  T1: forced cross-coupling at the operator level")
check(3, "additive operator mixed partial ||d2O/dadz||_F = %.3e vanishes "
      "(< %.0e, exact identity: additive O is linear in couplings)"
      % (d2O_add, TOL_ZERO),
      d2O_add < TOL_ZERO)
check(4, "multiplicative ||d2O/dadz||_F = %.3e is a forced cross-coupling "
      "(> %.0e), matching analytic ||K Rho_n Ren_n||_F = %.3e to 1e-6 rel"
      % (d2O_mult, TOL_NONZERO, cross_analytic),
      d2O_mult > TOL_NONZERO
      and abs(d2O_mult - cross_analytic) / cross_analytic < 1e-6)
check(5, "Hellmann-Feynman: additive per-sector response dlambda/da is "
      "invariant under z (drift %.3e < %.0e), multiplicative drifts more "
      "(%.3e)" % (add_inv, TOL_ZERO, mult_drift),
      add_inv < TOL_ZERO and mult_drift > add_inv)

print("S3  T2: clean decoupling limit (z -> 0 at canonical alpha)")
check(6, "additive z->0 recovers the bare reduced operator exactly "
      "(spectral residual %.3e < %.0e)" % (add_resid, TOL_ZERO),
      add_resid < TOL_ZERO)
check(7, "multiplicative z->0 distorts the surviving sector "
      "(||spec[mult] - spec[add]|| = %.3e > %.0e)"
      % (mult_vs_add, TOL_NONZERO),
      mult_vs_add > TOL_NONZERO)
check(8, "additivity is the unique composition meeting independent "
      "decoupling AND no forced cross-coupling (both T1 and T2 hold); "
      "honest boundary: operadic decoupling-respecting forms not excluded",
      d2O_add < TOL_ZERO and d2O_mult > TOL_NONZERO
      and add_resid < TOL_ZERO and mult_vs_add > TOL_NONZERO)

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)
