#!/usr/bin/env python3
"""verify_P310.py -- Verifier for Addendum 310 (P18-T2 sub-result (ii),
discharging the linear-coupling assumption of A308).

A308 reduced sub-result (ii) to one remaining assumption: each coupling c_j
enters as a LINEAR STRENGTH, c_j * S_j, rather than f_j(c_j) * S_j for a
nonlinear f_j. A308 retained it. A310 tests whether the LINEAR FORM is
FORCED by P18's own definition of the couplings.

P18 S4.2 ("Coefficient Values") treats alpha, gamma, zeta, beta as coupling
STRENGTHS, each fixed by geometric/dimensional consistency (alpha =
1/137.036 from int rho = alpha^{-1}; zeta = alpha^{5/4} from dimensional
analysis; etc.), and S4.1 writes the density sector as alpha * rho(r). A
coupling strength multiplying a fixed sector operator is, by definition, a
degree-1 homogeneous linear prefactor. The operator analog of "strength =
linear prefactor" is two-fold:

  (L1) the per-sector operator response is constant in the coupling:
       d O_j / d c_j = S_j, so the operator second derivative
       d^2 O_j / d c_j^2 = 0 (the zero operator). Linear is the unique
       member of the declared family whose own operator curvature vanishes.

  (L2) doubling the coupling doubles the sector's contribution:
       f_j(2c) = 2 f_j(c) (degree-1 homogeneity). Only the linear form is
       degree-1 homogeneous.

  (L3) combined with A308 (IC2 / no category bleed forces additivity within
       the linear class), additive composition is forced with both hinges
       grounded in P18 requirements: R1-R8 independent categories => IC2,
       and couplings are dimensional strengths => linear.

Recomputed from scratch on the density sector S1 = rho(r) [R5] over the
fixed (kinetic + V_self) operator K, with the renormalization sector
S2 = r d/dr [R7] cross-checked, at canonical alpha = 1/137.036,
zeta = alpha^{5/4}. Declared family of coupling forms (per-sector operator
O_j(c) = K + f_j(c) * S_j):

  LINEAR      f(c) = c            QUADRATIC  f(c) = c^2
  SQRT        f(c) = sqrt(c)      EXP-1      f(c) = exp(c) - 1

Tolerances: tol_zero = 1e-2 (vanishing operator curvature), tol_nonzero =
1e-1 (present curvature), tol_hom = 1e-6 (homogeneity residual). The linear
form's d^2 O/dc^2 is the exact zero operator; its measured value (~9e-4) is
finite-difference roundoff, floored by the radial d2 stencil's ~1/dr^2
entries; the smallest nonlinear curvature (quadratic, ~17) is four orders
above, so tol_zero = 1e-2 cleanly separates exact zero from genuine
curvature.

Honest boundary (Case C): "couplings are strengths" is P18's S4.2
definition made precise, not derived from deeper principles; a finite
declared family of f_j is tested, not every nonlinear form; sub-results
(i) = A304 and (iii) = A302 are unchanged; A294's re-typing stands;
P18-T2 is NOT closed.

  S1  Sector pieces and the declared coupling family   - checks 1-2
  S2  Per-sector operator curvature d^2 O/dc^2 (L1)      - checks 3-5
  S3  Degree-1 homogeneity (L2) and the assembled chain  - 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
A0 = 1.0 / 137.036
Z0 = A0 ** 1.25
ESELF = 13.177
M0 = 137.036
R0 = A0 ** 1.25
TOL_ZERO = 1e-2
TOL_NONZERO = 1e-1
TOL_HOM = 1e-6
C0 = A0
H = 1e-3


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


K, Rho, Ren = build_pieces(l=0)
S1 = Rho / float(np.linalg.norm(Rho, 2))      # density [R5]
S2 = Ren / float(np.linalg.norm(Ren, 2))      # renorm [R7]
nS1 = float(np.linalg.norm(S1, "fro"))

ffuncs = {
    "linear":     lambda c: c,
    "quadratic":  lambda c: c * c,
    "sqrt":       lambda c: math.sqrt(c),
    "exp_minus1": lambda c: math.exp(c) - 1.0,
}


def O1(fname, c):
    return K + ffuncs[fname](c) * S1


def d1_op(fname, c):
    return (O1(fname, c + H) - O1(fname, c - H)) / (2 * H)


def d2_op(fname, c):
    return (O1(fname, c + H) - 2 * O1(fname, c) + O1(fname, c - H)) / (H * H)


d1n = {f: float(np.linalg.norm(d1_op(f, C0), "fro")) for f in ffuncs}
d2n = {f: float(np.linalg.norm(d2_op(f, C0), "fro")) for f in ffuncs}

# analytic operator second-derivative references
ref_d2 = {
    "linear":     0.0,
    "quadratic":  float(np.linalg.norm(2.0 * S1, "fro")),
    "sqrt":       float(np.linalg.norm((0.25 * C0 ** -1.5) * S1, "fro")),
    "exp_minus1": float(np.linalg.norm(math.exp(C0) * S1, "fro")),
}

# degree-1 homogeneity residual max_c |f(2c) - 2 f(c)| / |2 f(c)|
cs = np.linspace(0.2 * C0, 5.0 * C0, 25)
hom = {}
for f, fn in ffuncs.items():
    rs = []
    for c in cs:
        rhs = 2.0 * fn(c)
        den = abs(rhs) if abs(rhs) > 1e-300 else 1.0
        rs.append(abs(fn(2.0 * c) - rhs) / den)
    hom[f] = float(max(rs))

is_hom = {f: hom[f] < TOL_HOM for f in ffuncs}
is_d2zero = {f: d2n[f] < TOL_ZERO for f in ffuncs}
nonlinear = [f for f in ffuncs if f != "linear"]

# renorm sector S2 under linear coupling: same vanishing
d2_S2 = float(np.linalg.norm(
    (K + (Z0 + H) * S2 - 2 * (K + Z0 * S2) + K + (Z0 - H) * S2) / (H * H),
    "fro"))

both = {f: (is_d2zero[f] and is_hom[f]) for f in ffuncs}
linear_unique = all(both[f] == (f == "linear") for f in ffuncs)

print("S1  Sector pieces and the declared coupling family")
check(1, "two distinct sectors built: density S1 [R5], renorm S2 [R7] "
      "(norms > 0, S1 != S2); canonical (alpha, zeta) = (%.6g, %.6g)"
      % (A0, Z0),
      np.linalg.norm(S1) > 0 and np.linalg.norm(S2) > 0
      and not np.allclose(S1, S2))
check(2, "declared family of 4 coupling forms built (linear, quadratic, "
      "sqrt, exp-1); the linear dO/dc = %.4e matches the constant operator "
      "||S1||_F = %.4e (linear response is the bare sector, the strength "
      "reading)" % (d1n["linear"], nS1),
      len(ffuncs) == 4 and abs(d1n["linear"] - nS1) / nS1 < 1e-3)

print("S2  Per-sector operator curvature d^2 O/dc^2 = ||d2 O/dc2||_F (L1)")
check(3, "linear d^2 O/dc^2 = %.3e is the zero operator (< %.0e: the linear "
      "per-sector response is constant = S1, exact identity; measured value "
      "is finite-difference roundoff)" % (d2n["linear"], TOL_ZERO),
      d2n["linear"] < TOL_ZERO)
check(4, "quadratic d^2 O/dc^2 = %.3e and exp-1 d^2 O/dc^2 = %.3e are nonzero "
      "curvature (> %.0e), matching analytic ||2 S1||_F = %.3e and "
      "e^{c0}||S1||_F = %.3e to 1e-4 rel"
      % (d2n["quadratic"], d2n["exp_minus1"], TOL_NONZERO,
         ref_d2["quadratic"], ref_d2["exp_minus1"]),
      d2n["quadratic"] > TOL_NONZERO and d2n["exp_minus1"] > TOL_NONZERO
      and abs(d2n["quadratic"] - ref_d2["quadratic"]) / ref_d2["quadratic"] < 1e-4
      and abs(d2n["exp_minus1"] - ref_d2["exp_minus1"]) / ref_d2["exp_minus1"] < 1e-4)
check(5, "sqrt d^2 O/dc^2 = %.3e is nonzero curvature (> %.0e), matching "
      "analytic (1/4)c0^{-3/2}||S1||_F = %.3e to 1e-2 rel (finite-difference "
      "O(H^2) residual on the strongly-curved sqrt)"
      % (d2n["sqrt"], TOL_NONZERO, ref_d2["sqrt"]),
      d2n["sqrt"] > TOL_NONZERO
      and abs(d2n["sqrt"] - ref_d2["sqrt"]) / ref_d2["sqrt"] < 1e-2)

print("S3  Degree-1 homogeneity (L2) and the assembled chain (L3)")
check(6, "linear is degree-1 homogeneous: max|f(2c)-2f(c)|/|2f(c)| = %.3e "
      "(< %.0e); every nonlinear fails (quadratic %.3e, sqrt %.3e, exp-1 %.3e "
      "all >= %.0e)"
      % (hom["linear"], TOL_HOM, hom["quadratic"], hom["sqrt"],
         hom["exp_minus1"], TOL_HOM),
      is_hom["linear"] and all(not is_hom[f] for f in nonlinear))
check(7, "every nonlinear form is disqualified on BOTH criteria (nonzero "
      "operator curvature AND fails degree-1 homogeneity); linear is the "
      "UNIQUE form in the family that is both d^2 O/dc^2 = 0 and degree-1 "
      "homogeneous",
      all((d2n[f] > TOL_NONZERO) and (not is_hom[f]) for f in nonlinear)
      and linear_unique)
check(8, "renorm sector S2 [R7] under linear coupling has d^2 O/dc^2 = %.3e "
      "(zero, < %.0e), confirming the vanishing is sector-independent; "
      "assembled chain: a coupling strength (degree-1 homogeneous, no "
      "operator curvature) is the LINEAR form, so with A308 (IC2 forces "
      "additivity within the linear class) additive composition is forced "
      "with BOTH hinges grounded in P18 requirements; honest boundary: "
      "'couplings are strengths' is P18 S4.2 made precise not derived, "
      "finite declared family, (i)=A304/(iii)=A302 unchanged, A294 stands, "
      "P18-T2 NOT closed" % (d2_S2, TOL_ZERO),
      d2_S2 < TOL_ZERO and linear_unique
      and d2n["linear"] < TOL_ZERO and is_hom["linear"]
      and all((d2n[f] > TOL_NONZERO) and (not is_hom[f]) for f in nonlinear))

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