#!/usr/bin/env python3
"""
verify_P154.py — P154 Laplace-type spectral transform candidates for G2 root measure.
Investigates whether any natural smearing of the G2 discrete root spectral measure
produces a second moment matching BREATH_PERIOD/9, closing gap G1.
"""

import math
import numpy as np

# ─── constants ────────────────────────────────────────────────────────────────
ALPHA_INV   = 4*math.pi**3 + math.pi**2 + math.pi   # ≈ 137.036
BREATH_PERIOD = math.pi * ALPHA_INV                   # ≈ 430.806
ALPHA       = 1.0 / ALPHA_INV                         # fine-structure constant ≈ 0.007297

# G2 root projections onto H_s (coroot of short simple root alpha_1)
PROJ = [-3, -3, -2, -1, -1, 0, 0, 1, 1, 2, 3, 3]
N    = len(PROJ)                     # 12 roots
PROJ_SQ = [p**2 for p in PROJ]
J_discrete = sum(PROJ_SQ)           # discrete Killing form = 48
TARGET = BREATH_PERIOD / 9          # the transcendental target ≈ 47.867

G1 = (432 - BREATH_PERIOD) / BREATH_PERIOD

print("=" * 65)
print("P154 SPECTRAL LAPLACE VERIFICATION")
print("=" * 65)
print(f"ALPHA_INV        = {ALPHA_INV:.10f}")
print(f"BREATH_PERIOD    = {BREATH_PERIOD:.10f}")
print(f"BREATH_PERIOD/9  = {TARGET:.10f}   ← transcendental target")
print(f"J_discrete       = {J_discrete}          ← Killing form (integer)")
print(f"G1 = (432-BP)/BP = {G1:.8e}  ({G1*100:.6f}%)")
print()


# ─── utility ──────────────────────────────────────────────────────────────────
def bisect(f, a, b, tol=1e-14, maxiter=200):
    """Bisection root-finder; assumes f(a)*f(b) <= 0."""
    fa, fb = f(a), f(b)
    if fa == 0: return a
    if fb == 0: return b
    if fa * fb > 0:
        return None          # no sign change
    for _ in range(maxiter):
        mid = 0.5*(a + b)
        fm  = f(mid)
        if abs(fm) < tol or (b - a) < tol:
            return mid
        if fa * fm < 0:
            b, fb = mid, fm
        else:
            a, fa = mid, fm
    return 0.5*(a + b)


PASS = FAIL = 0
_N = 0


def check(label, result, tgt, note=""):
    global PASS, FAIL, _N
    ratio  = result / tgt if tgt != 0 else float('nan')
    ok     = abs(ratio - 1) < 0.01
    PASS  += ok
    FAIL  += (not ok)
    _N    += 1
    tag    = f"  [{note}]" if note else ""
    print(f"  [{'PASS' if ok else 'FAIL'}] {_N:>2}. Candidate {label}: "
          f"result = {result:.8f}, target = {tgt:.8f}, "
          f"ratio = {ratio:.8f}{tag}")


# ══════════════════════════════════════════════════════════════════════════════
# CANDIDATE A — Gaussian smearing at σ² = π
# M2_A = Σ_α [⟨α,H_s⟩² + π] = 48 + 12π
# ══════════════════════════════════════════════════════════════════════════════
print("─" * 65)
print("CANDIDATE A — Gaussian smearing σ² = π")
M2_A = J_discrete + N * math.pi
print(f"  M2_A = 48 + 12π = {M2_A:.8f}")
check("A", M2_A, TARGET, "smearing variance = π")

# ══════════════════════════════════════════════════════════════════════════════
# CANDIDATE B — MGF-Laplace:  L(s) = Σ_α exp(-s·⟨α,H_s⟩²)
# Range: L(s) ∈ (2, 12] for s ≥ 0  (2 zero-projection roots floor it at 2)
# Seek s* where the weighted mean-square -L'(s*)/L(s*) equals TARGET/N
# ══════════════════════════════════════════════════════════════════════════════
print("\n─" * 2 + "─" * 63)
print("CANDIDATE B — Laplace MGF  L(s) = Σ_α exp(-s·⟨α,H_s⟩²)")

def L(s):
    return sum(math.exp(-s * p2) for p2 in PROJ_SQ)

def L_prime(s):
    return sum(-p2 * math.exp(-s * p2) for p2 in PROJ_SQ)

def L_double(s):
    return sum(p2**2 * math.exp(-s * p2) for p2 in PROJ_SQ)

print(f"  L(0)   = {L(0):.4f}   (= N = {N})")
print(f"  lim_∞ L(s) → {L(1000):.4f}  (2 zero-projection roots)")
print(f"  Range of L(s) on [0,∞): [{L(1000):.4f}, {L(0):.4f}]")
print(f"  -L'(0)/L(0) = {-L_prime(0)/L(0):.6f} = J_discrete/N = 48/12 = 4")
print()

# Target for weighted mean-square:  TARGET/N = BREATH_PERIOD/108
wms_target = TARGET / N
print(f"  Seek s* where -L'(s*)/L(s*) = BREATH_PERIOD/108 = {wms_target:.10f}")

f_B = lambda s: -L_prime(s)/L(s) - wms_target
# f_B(0)   = 4 - 3.9888... > 0
# f_B(∞)   → 0 - 3.9888... < 0  (non-zero-projection terms vanish, leaving zero)
# Actually at large s, L→2 (zero-proj), L'→0, so -L'/L → 0. Sign change exists.
s_B = bisect(f_B, 1e-12, 50.0)
if s_B is not None:
    wms_at_sB = -L_prime(s_B)/L(s_B)
    print(f"  s* = {s_B:.10f}")
    print(f"  s* / (1/BREATH_PERIOD) = {s_B * BREATH_PERIOD:.6f}")
    print(f"  s* / (π/432)           = {s_B / (math.pi/432):.6f}")
    print(f"  s* / (G1/4π)           = {s_B / (G1/(4*math.pi)):.6f}")
    print(f"  Analytical approx (linear): s* ≈ (J/N - TARGET/N)/0.5")
    # Second-order correction for -L'/L near s=0:
    # d/ds[-L'/L] = [-L''·L + (L')²]/L² = [-(Σp⁴)·N + (Σp²)²]/N² at s=0
    sum_p4 = sum(p**4 for p in PROJ)
    deriv_wms_0 = (-L_double(0)*L(0) + L_prime(0)**2) / L(0)**2
    s_B_linear = (4.0 - wms_target) / (-deriv_wms_0)
    print(f"  d(-L'/L)/ds|_{{s=0}} = {deriv_wms_0:.6f}")
    print(f"  Linear estimate s* ≈ {s_B_linear:.8f}")
    print(f"  s* × BREATH_PERIOD = {s_B * BREATH_PERIOD:.8f}")
    print(f"  s* × 432           = {s_B * 432:.8f}")
    check("B", wms_at_sB * N, TARGET, f"s*={s_B:.6e}")

# ══════════════════════════════════════════════════════════════════════════════
# CANDIDATE C — Bilateral Laplace / MGF
# M(s) = 4cosh(3s) + 2cosh(2s) + 4cosh(s) + 2;  M''(0) = 48, M(0) = 12
# Seek s* where M''(s*)/M(s*) = BREATH_PERIOD/432 × (J/N) = BREATH_PERIOD/108
# But note: d/ds[M''/M] has the same monotone-increase problem.
# ══════════════════════════════════════════════════════════════════════════════
print("\n─" * 2 + "─" * 63)
print("CANDIDATE C — Bilateral MGF  M(s) = 4cosh(3s)+2cosh(2s)+4cosh(s)+2")

def M(s):
    return 4*math.cosh(3*s) + 2*math.cosh(2*s) + 4*math.cosh(s) + 2

def M_pp(s):
    return 36*math.cosh(3*s) + 8*math.cosh(2*s) + 4*math.cosh(s)

print(f"  M(0)   = {M(0):.4f}  (= N = 12)")
print(f"  M''(0) = {M_pp(0):.4f}  (= J_discrete = 48)")
print(f"  M''(0)/M(0) = {M_pp(0)/M(0):.6f}")
print()

# As |s|→∞, M and M'' both grow as 4cosh(3s) → M''/M → 36/4 = 9
# So M''/M increases monotonically from 4 (s=0) to 9 (s→∞).
# The ratio BREATH_PERIOD/432 ≈ 0.998 < 4: no solution for s > 0.
ratio_C_target = BREATH_PERIOD / 432
print(f"  Ratio target BREATH_PERIOD/432 = {ratio_C_target:.8f}")
print(f"  M''(s)/M(s) range: [{M_pp(0)/M(0):.4f}, 9.000]  (monotone increasing)")
print(f"  Target {ratio_C_target:.4f} < minimum 4.000: no real s* satisfies M''/M = BP/432")

# Corrected interpretation: seek s* where M''(s*)/M(s*) = (TARGET/N) × (M(0)/M(0))
# i.e., weighted second moment = TARGET/N
# Since M'' is even and increases with |s|, and M''(0)/M(0) = 4 > TARGET/N ≈ 3.989:
# No solution for s ≠ 0 with M''/M < 4 (ratio can only grow).
print(f"  Alt target M''(s*)/M(s*) = TARGET/N = {TARGET/N:.8f} < 4: also no solution")
check("C", M_pp(0)/M(0), TARGET/N, "evaluated at s=0; no s* exists")

# ══════════════════════════════════════════════════════════════════════════════
# CANDIDATE D — S³ heat-kernel smearing
# M2_D(t) = J_discrete + 2t;  solve for t*: 48 + 2t* = BREATH_PERIOD/9
# ══════════════════════════════════════════════════════════════════════════════
print("\n─" * 2 + "─" * 63)
print("CANDIDATE D — 1D heat-kernel smearing  M2_D(t) = 48 + 2t")

t_star = (TARGET - J_discrete) / 2.0
M2_D_tstar = J_discrete + 2 * t_star   # = TARGET by construction
print(f"  t* = (BP/9 - 48)/2 = {t_star:.10f}")
print(f"  t* = -(432 - BP)/18 = {-(432-BREATH_PERIOD)/18:.10f}")
print(f"  t* / π  = {t_star/math.pi:.8f}")
print(f"  |t*| = (432-BP)/18 = G1 · BREATH_PERIOD/18 = {G1*BREATH_PERIOD/18:.10f}")
print(f"  t* < 0: heat-kernel smearing always increases M2; t* < 0 is unphysical.")
print(f"  M2_D(t=π)        = {J_discrete + 2*math.pi:.8f}  (too large)")
print(f"  M2_D(t=BREATH_PERIOD) = {J_discrete + 2*BREATH_PERIOD:.6f}  (far too large)")

# By construction t* gives the target, but it is negative — unphysical.
check("D", M2_D_tstar, TARGET, "t* < 0: unphysical")

# ══════════════════════════════════════════════════════════════════════════════
# CANDIDATE E — Continuous uniform measure on circle of radius R
# Each root on a circle: continuous M2 = R²/2.
# Seek R such that (R²/2) × N = BREATH_PERIOD
# ══════════════════════════════════════════════════════════════════════════════
print("\n─" * 2 + "─" * 63)
print("CANDIDATE E — Uniform circular measure; total M2 = N·R²/2 = BREATH_PERIOD")

R_sq = BREATH_PERIOD / (N / 2.0)          # = 2·BREATH_PERIOD/N = BREATH_PERIOD/6
R    = math.sqrt(R_sq)
print(f"  R² = BREATH_PERIOD/6 = {R_sq:.10f}")
print(f"  R  = √(π·ALPHA_INV/6) = {R:.10f}")
print(f"  R / π  = {R/math.pi:.8f}")
print(f"  R² / π = {R_sq/math.pi:.8f}  (= ALPHA_INV/6 = {ALPHA_INV/6:.8f})")
print(f"  Note: R² = π·α⁻¹/6 — transcendental, no simple integer form")
M2_E = (R_sq / 2.0) * N
# This matches BREATH_PERIOD by construction; target for this candidate is BREATH_PERIOD
check("E", M2_E, BREATH_PERIOD, "matches BP by construction, not BP/9")

# What R would give M2 = TARGET = BREATH_PERIOD/9?
R_sq_9 = TARGET / (N / 2.0)
R_9    = math.sqrt(R_sq_9)
print(f"\n  To match BREATH_PERIOD/9: R² = {R_sq_9:.8f}, R = {R_9:.8f}")
print(f"  R_9 / π = {R_9/math.pi:.8f}")
check("E'", (R_sq_9/2)*N, TARGET, "BP/9 version: R²=BP/54")

# ══════════════════════════════════════════════════════════════════════════════
# CANDIDATE F — Weyl character at imaginary argument
# χ_adj(s) = 4cos(3s) + 2cos(2s) + 4cos(s) + 4
# (12 root contributions + 2 Cartan generators → dim G₂ = 14)
# Seek s* in [0,π] where χ_adj(s*) = α (fine-structure constant ≈ 0.007297)
# ══════════════════════════════════════════════════════════════════════════════
print("\n─" * 2 + "─" * 63)
print("CANDIDATE F — Weyl character  χ_adj(s) = 4cos(3s)+2cos(2s)+4cos(s)+4")

def chi_adj(s):
    return 4*math.cos(3*s) + 2*math.cos(2*s) + 4*math.cos(s) + 4

print(f"  χ_adj(0) = {chi_adj(0):.4f}  (= dim G₂ adjoint = 14)")
print(f"  χ_adj(π) = {chi_adj(math.pi):.4f}")
print(f"  Min on [0,π]: scanning...")

s_grid = np.linspace(0, math.pi, 100000)
chi_grid = np.array([chi_adj(float(s)) for s in s_grid])
print(f"  min(χ_adj) on [0,π] ≈ {chi_grid.min():.6f}  at s ≈ {float(s_grid[chi_grid.argmin()]):.6f}")

# Values at s = π/n for n = 1..12
print()
print("  χ_adj(π/n) for n = 1 … 12:")
for n in range(1, 13):
    s = math.pi / n
    v = chi_adj(s)
    print(f"    n={n:2d}: s=π/{n:<2d} ({s:.5f}),  χ_adj = {v:+.8f}")

# Zeros of χ_adj in [0, π]
print()
print("  Zeros of χ_adj in [0,π]:")
zeros_F = []
for i in range(len(chi_grid)-1):
    if chi_grid[i]*chi_grid[i+1] < 0:
        s_zero = bisect(chi_adj, float(s_grid[i]), float(s_grid[i+1]))
        if s_zero is not None:
            zeros_F.append(s_zero)
            print(f"    s_zero = {s_zero:.10f},  s/π = {s_zero/math.pi:.10f}")

# Crossings of χ_adj(s) = ALPHA (the fine-structure constant α ≈ 0.007297)
print()
target_chi = ALPHA
print(f"  Crossings of χ_adj(s) = α = 1/ALPHA_INV = {target_chi:.8e}:")
chi_crossings = []
for i in range(len(chi_grid)-1):
    if (chi_grid[i] - target_chi)*(chi_grid[i+1] - target_chi) < 0:
        sc = bisect(lambda s: chi_adj(s) - target_chi, float(s_grid[i]), float(s_grid[i+1]))
        if sc is not None:
            chi_crossings.append(sc)
            ratio_to_pi = sc / math.pi
            print(f"    s* = {sc:.10f},  s*/π = {ratio_to_pi:.10f}")
            # Check proximity to simple fractions of π
            for denom in range(2, 24):
                for num in range(1, denom):
                    frac = num / denom
                    if abs(ratio_to_pi - frac) < 5e-4:
                        print(f"      ≈ {num}π/{denom}")

if chi_crossings:
    check("F", chi_adj(chi_crossings[0]), target_chi, "first crossing")
else:
    print("  No crossings found in [0,π].")
    check("F", chi_adj(0.0), target_chi, "no crossing; value at s=0")

# ══════════════════════════════════════════════════════════════════════════════
# SUMMARY TABLE
# ══════════════════════════════════════════════════════════════════════════════
print()
print("=" * 65)
print("SUMMARY")
print("=" * 65)
print(f"G1               = (432 - BREATH_PERIOD)/BREATH_PERIOD = {G1:.8e}")
print(f"BREATH_PERIOD/9  = {TARGET:.10f}")
print(f"J_discrete       = {J_discrete}   (integer; discrete Killing form)")
print(f"48 - BP/9        = {J_discrete - TARGET:.10f}  (size of gap in moment space)")
print()
check("A", M2_A, TARGET, "48+12π")
wms_B = -L_prime(s_B)/L(s_B) * N if s_B else float('nan')
check("B", wms_B, TARGET, f"s*={s_B:.4e}" if s_B else "no s*")
check("C", M_pp(0)/M(0) * N, TARGET, "no finite s*; evaluated at s=0")
check("D", M2_D_tstar, TARGET, "t*<0 unphysical; trivially by construction")
check("E", (R_sq_9/2)*N, TARGET, "R²=BP/54 by construction")
if chi_crossings:
    check("F", chi_adj(chi_crossings[0]), target_chi, "Weyl char crosses α; target is α not BP/9")
else:
    check("F", chi_adj(0), target_chi)

print()
print("CONCLUSION: No candidate naturally closes G1.")
print(f"G1 = (432-π·α⁻¹)/(π·α⁻¹) requires transcendental input beyond G₂ geometry.")
print(f"The gap {G1*100:.6f}% remains open after P154.")

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")