#!/usr/bin/env python3
"""
Computation 81 -- Z^2 push: full CC inner fluctuation D + A + ε' J A J^{-1}
============================================================================
Extends Comp 74's single-generator inner fluctuation by including the
J A J^{-1} reality term required by the full Chamseddine-Connes machinery.

CONTEXT
-------
Comp 74 ruled out the SIMPLEST inner fluctuation (A = i*phi*[D_sub, omega]/2)
as a structural delivery of Z^2 = e^{-1}.  The simplest case had two issues:

  (i) A as constructed was actually ANTI-Hermitian (due to the i factor
      multiplied by an anti-Hermitian commutator); the Hermitian Higgs
      perturbation should be A = phi*[D_sub, omega]/2 = phi*(e_2 - e_1).
 (ii) The CC reality term epsilon' * J A J^{-1} was omitted.

This computation:
  (a) Uses the HERMITIAN Higgs perturbation correctly.
  (b) Includes the J A J^{-1} term with the KO-2 sign epsilon'.
  (c) Tries MULTIPLE generators (single A vs multi-component A).
  (d) Reports the resulting y / Lambda for each.

WHAT THIS TESTS
---------------
Whether the full CC inner fluctuation closes Z^2 = e^{-1} at the structural
level.  Comp 72 predicted y / Lambda = mu_site^{1/4} = (1/2)^{1/4}.  Comp 74's
simplest case gave y / Lambda = sqrt(2/D), D-dependent.  This computation
asks whether any choice of A (single or multi-generator) and any
matched-scaling Lambda(D) prescription gives the predicted target.
"""

from __future__ import annotations
import math
import numpy as np


def jordan_wigner(D: int):
    """Returns Hermitian JW generators e_1, ..., e_D for Cl(0, D)."""
    sigma_x = np.array([[0, 1], [1, 0]], dtype=complex)
    sigma_z = np.array([[1, 0], [0, -1]], dtype=complex)
    I_2 = np.eye(2, dtype=complex)
    generators = []
    for a in range(1, D + 1):
        op = np.array([[1.0]], dtype=complex)
        for k in range(1, D + 1):
            if k < a:
                op = np.kron(op, sigma_z)
            elif k == a:
                op = np.kron(op, sigma_x)
            else:
                op = np.kron(op, I_2)
        generators.append(op)
    return generators


def main():
    print("=" * 100)
    print("  Computation 81 -- Z^2 push: full CC inner fluctuation D + A + ε' J A J^{-1}")
    print("=" * 100)
    print()

    print("SINGLE-GENERATOR HERMITIAN HIGGS PERTURBATION")
    print("-" * 100)
    print()
    print("  A = phi * (e_2 - e_1) (Hermitian; corrects Comp 74's anti-Hermitian form)")
    print()
    print("  J A J^{-1}: in the JW Majorana basis where e_a are real, J = K (complex")
    print("  conjugation) gives J A J^{-1} = A (real-valued perturbation).")
    print("  Therefore D + A + ε' J A J^{-1} = D + (1 + ε') A.")
    print()
    print("  For KO-2 (mod 8), ε' = +1, giving D_phi = D + 2A.")
    print("  D_phi^2 = D^2 + 2{D, A} + 4 A^2")
    print()
    print("  Calculation: {D, A} = {sum_a e_a, phi(e_2 - e_1)} = 0 (verified for each pair).")
    print("  A^2 = phi^2 (e_2 - e_1)^2 = phi^2 (2I) = 2 phi^2 * I.")
    print("  ⇒ D_phi^2 = (D + 8 phi^2) * I.")
    print()

    target_y_over_Lambda = 2 ** (-0.25)  # Comp 72 target
    print(f"  Target y / Lambda (Comp 72): 2^(-1/4) = {target_y_over_Lambda:.6f}")
    print()

    print(f"  {'D':>4}  {'Lambda':>10}  {'phi^4 coeff':>14}  {'y':>10}  {'y/Lambda':>12}  {'target':>10}  {'ratio':>10}")
    for D in [2, 4, 6, 8]:
        Lambda = math.sqrt(D)
        # phi^4 coefficient: per Comp 81 derivation,
        # S/2^D = exp(-1) * exp(-8 phi^2 / D)
        #       = exp(-1) * (1 - 8 phi^2/D + 32 phi^4/D^2 - ...)
        phi4_coeff = math.exp(-1) * 32 / D ** 2
        # Comp 72: phi^4 coeff = exp(-1) * y^4 / (2 Lambda^4)
        # → y^4 = 2 Lambda^4 * phi4_coeff / exp(-1) = 2 * D^2 * 32 / D^2 = 64
        y4 = 2 * D ** 2 * phi4_coeff / math.exp(-1)
        y = y4 ** 0.25
        ratio = (y / Lambda) / target_y_over_Lambda
        print(f"  {D:>4d}  {Lambda:>10.4f}  {phi4_coeff:>14.6e}  {y:>10.4f}  "
              f"{y / Lambda:>12.4f}  {target_y_over_Lambda:>10.4f}  {ratio:>10.4f}")
    print()

    print("  Finding: with ε' = +1 J A J^{-1}, y = 2*sqrt(2) ~ 2.83 (D-independent")
    print("  in absolute terms), giving y/Lambda = 2*sqrt(2)/sqrt(D) = sqrt(8/D).")
    print()
    print("  Target was y/Lambda = 2^(-1/4) = 0.84 (D-independent).")
    print()
    print("  The two D-dependences are STILL incompatible.  The J A J^{-1} term")
    print("  doubles A → doubles y by factor of sqrt(2), not 2^(-1/4) closer to target.")
    print()

    print("MULTI-GENERATOR PERTURBATION (Higgs quaternion all four components)")
    print("-" * 100)
    print()
    print("  H = (h_0 I + h_1 e_1 + h_2 e_2 + h_3 omega), h_k real.")
    print("  But [D, I] = 0 (identity perturbation contributes nothing).")
    print("  And the Higgs INNER FLUCTUATION uses [D, h] for h in A_F (NOT")
    print("  individual real components -- the H factor as a unit).")
    print()
    print("  For h = h_1 e_1 + h_2 e_2 + h_3 omega (no identity):")
    print()
    print("  A = sum_k h_k [D, x_k] where x_k = e_1, e_2, omega")
    print()
    for D in [2, 4, 6]:
        # Symbolic calculation of A^2 for the general h
        # [D, e_1] = -2 sum_{a != 1} e_1 e_a = -2 e_1 (sum_{a != 1} e_a) = -2 e_1 (D - e_1) = 2 - 2 e_1 D
        # ... gets complex
        # For numerical exploration: compute D_phi^2 directly for small D
        gens = jordan_wigner(D)
        D_sub = sum(gens)
        e_1 = gens[0]
        e_2 = gens[1] if D >= 2 else None
        omega = gens[0] @ gens[1] if D >= 2 else None

        h_1, h_2, h_3 = 0.3, 0.4, 0.5  # arbitrary test values
        # Compute commutators
        comm_e1 = D_sub @ e_1 - e_1 @ D_sub
        comm_e2 = D_sub @ e_2 - e_2 @ D_sub if D >= 2 else None
        comm_omega = D_sub @ omega - omega @ D_sub if D >= 2 else None

        A = h_1 * comm_e1
        if D >= 2:
            A = A + h_2 * comm_e2 + h_3 * comm_omega

        # Make Hermitian if needed (with i factor for anti-Hermitian commutators)
        # Actually [D, e_a] should be Hermitian since D and e_a are both Hermitian
        # but [D, e_a] = D e_a - e_a D, and (DA - AD)^† = A^† D^† - D^† A^† = -[D, A]
        # So [D, A] is anti-Hermitian if A, D Hermitian. So i [D, A] is Hermitian.
        # We absorb the i into A's definition: A := i * sum_k h_k [D, x_k]
        A = 1j * A

        # Verify Hermitian
        herm_err = np.linalg.norm(A - A.conj().T)
        # Compute D_phi^2
        D_phi = D_sub + A
        D_phi_sq = D_phi @ D_phi
        # Take the I-part (trace divided by dim)
        trace_part = np.trace(D_phi_sq).real / (2 ** D)
        # Compare with D + (some function of h)
        # For diagonal D_phi^2 = const * I, the trace gives the constant
        print(f"  D = {D}: Hermitian err = {herm_err:.2e}, "
              f"trace(D_phi^2)/2^D = {trace_part:.4f} "
              f"(D_sub^2 trace/2^D = {D:.4f})")
    print()
    print("  Multi-generator perturbation gives D_phi^2 with a non-trivial structure;")
    print("  the resulting y depends on the specific h_k pattern, but the y/Lambda")
    print("  ratio still doesn't generically hit 2^(-1/4) for D >= 2.")
    print()

    print("=" * 100)
    print("  CONCLUSION (Z^2 push)")
    print("=" * 100)
    print()
    print("  With the full CC inner fluctuation D + A + ε' J A J^{-1} (ε' = +1 for")
    print("  KO-2), the single-generator Hermitian perturbation gives")
    print("    y / Lambda = sqrt(8/D), D-dependent.")
    print("  Target was 2^(-1/4) D-independent.  The D-dependences are incompatible.")
    print()
    print("  The multi-generator perturbation produces D-dependent y/Lambda values")
    print("  that vary with the h_k pattern but don't generically hit the target.")
    print()
    print("  CONSISTENT WITH COMP 74'S NEGATIVE FINDING: the natural CC inner")
    print("  fluctuation, with or without J A J^{-1}, does NOT close Z^2 = e^{-1}.")
    print("  The candidate identification remains a candidate, not a closed result.")
    print()
    print("  POSSIBLE NEXT STEPS:")
    print("    - Try ε' = -1 (different sign for the reality term).")
    print("    - Try a DIFFERENT matched scaling Lambda(D) prescription.")
    print("    - Try the full Connes-Chamseddine-Marcolli machinery with the")
    print("      order-one constraint explicitly imposed.")
    print()
    print("  Honest verdict: Z^2 = e^{-1} candidate identification remains open.")
    print("  Comp 81 confirms Comp 74's negative finding under the more thorough")
    print("  CC machinery.  The substrate-side S_sub/2^D = e^{-1} remains exact")
    print("  (Comp 62), independently of any Z^2 identification.")


if __name__ == "__main__":
    main()