"""
CCT Automaton with Black Hole Energy Monitor
============================================

This implementation extends the gradient descent factorization with:
1. J-coupling computation (ζ-Γ-W network)
2. Energy eigenvalue tracking from oscillating solutions
3. Black hole mass and Hawking temperature calculation
4. Real-time visualization of energy modes

Based on the CCT-ODE framework connecting:
- ζ(s) → Information encoding (primes)
- Γ(s) → Thermodynamic entropy
- W(z) → Energy-time relationships
"""

import torch
import math
import numpy as np
from scipy.special import jv as bessel_j0
import matplotlib.pyplot as plt
from matplotlib.patches import FancyBboxPatch
import warnings
warnings.filterwarnings('ignore')

# ============================================================================
# PART 1: UTILITY FUNCTIONS
# ============================================================================

def is_prime_torch(n: int) -> bool:
    """
    Deterministic Miller-Rabin primality test.
    Valid for n < 2^64.
    """
    if n < 2:
        return False
    
    small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
    for p in small_primes:
        if n % p == 0:
            return n == p
    
    # Write n-1 = d * 2^s
    d = n - 1
    s = 0
    while d % 2 == 0:
        d //= 2
        s += 1
    
    # Test bases
    for a in [2, 325, 9375, 28178, 450775, 9780504, 1795265022]:
        if a % n == 0:
            continue
        
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        
        for _ in range(s - 1):
            x = (x * x) % n
            if x == n - 1:
                break
        else:
            return False
    
    return True


def isqrt(n: int) -> int:
    """Integer square root."""
    return int(math.isqrt(n))


def compute_J_couplings(s_value: complex, M_normalized: float = 1.0):
    """
    Compute the ζ-Γ-W coupling constants for a given complex s.
    
    J_{ζΓ}(s) = |Γ(1-s)| · |ζ'(s)/ζ(s)| / (|ψ(s)| · √(2π) · |s|^{Re(s)-½})
    J_{ζW}(s) = |2^{1-s} - 1|
    J_{ΓW}(s) = |π · ζ(s) · (2^{1-s} - 1)| / (|Γ(s)| · |sin(π s)|)
    
    Args:
        s_value: Complex value of s (typically s = √(d² + 4c))
        M_normalized: Normalized mass parameter for BH observables
    
    Returns:
        dict with J_zg, J_zw, J_gw, and derived quantities
    """
    s = complex(s_value)
    
    # J_{ζW} - simplest coupling
    J_zw = abs(2**(1 - s) - 1)
    
    # J_{ζΓ} - requires zeta derivatives (approximated)
    # For numerical stability, use approximation based on s position
    try:
        # Approximate ζ(s) and ζ'(s) using series expansion near critical strip
        # ζ(s) ≈ 1 + 2^{-s} + 3^{-s} for Re(s) > 1
        # For general s, use reflection formula relationship
        
        # Simplified J_{ζΓ} based on distance to pole at s=1
        distance_to_pole = abs(s - 1)
        
        if distance_to_pole < 0.1:
            # Near pole: use pole approximation
            J_zg = 1.0 / distance_to_pole
        else:
            # Away from pole: use functional equation approximation
            # J_zg ≈ |Γ(1-s)| scaling with distance from critical line
            Gamma_approx = bessel_j0(abs(s)) if abs(s) > 0.5 else 1.0 / abs(s)
            psi_approx = math.log(abs(s)) if abs(s) > 1 else -0.5772 - 1/s.real
            J_zg = Gamma_approx * abs(s)**0.5 / max(abs(psi_approx), 0.1)
    except:
        J_zg = 1.0
    
    # J_{ΓW} - couples Gamma and W
    try:
        sin_pi_s = math.sin(math.pi * s.real)
        if abs(sin_pi_s) < 1e-10:
            J_gw = 10.0  # Near zero of sin(πs) - enhanced coupling
        else:
            Gamma_s = bessel_j0(abs(s)) if abs(s) > 0.5 else 1.0 / abs(s)
            J_gw = abs(math.pi * J_zw) / (abs(Gamma_s) * abs(sin_pi_s) + 1e-10)
    except:
        J_gw = 1.0
    
    # Normalize and add phase information for oscillation tracking
    J_zg_complex = J_zg * math.exp(1j * (s.imag * 0.1))
    J_zw_complex = J_zw * math.exp(1j * (s.real * 0.05))
    J_gw_complex = J_gw * math.exp(1j * (abs(s) * 0.08))
    
    # Coupling product (geometric mean)
    J_product = (J_zg * J_zw * J_gw) ** (1/3)
    
    return {
        'J_zg': J_zg,
        'J_zw': J_zw,
        'J_gw': J_gw,
        'J_zg_complex': J_zg_complex,
        'J_zw_complex': J_zw_complex,
        'J_gw_complex': J_gw_complex,
        'J_product': J_product,
        's_value': s
    }


# ============================================================================
# PART 2: BLACK HOLE ENERGY MONITOR
# ============================================================================

class BlackHoleEnergyMonitor:
    """
    Monitor for computing energy eigenvalues from J-network oscillations.
    
    The oscillating J-values (J_ζΓ, J_ζW, J_ΓW) are interpreted as
    energy-bearing modes of the mathematical black hole network.
    
    Key physics:
    - E_n = (1/2)[1 - J_0(2πA)cos(2πs_0)]  (energy quantum)
    - M_n = E_n · m_P / √(J_ζΓ·J_ζW·J_ΓW)   (BH mass)
    - T_H = (E_n·ħc³)/(8πGM_n·k_B)         (Hawking temp)
    """
    
    def __init__(self, name: str = "BH-1"):
        self.name = name
        self.energy_history = []
        self.mass_history = []
        self.temp_history = []
        self.J_zg_history = []
        self.J_zw_history = []
        self.J_gw_history = []
        self.s_trajectory = []
        self.loss_history = []
        self.time_step = 0
        
        # Physical constants
        self.m_P = 2.176e-8      # Planck mass (kg)
        self.m_P_solar = 2.18e-8 / 1.989e30  # in solar masses
        self.hbar = 1.055e-34
        self.c = 3e8
        self.G = 6.674e-11
        self.k_B = 1.381e-23
        
        # State
        self.current_E_n = 0.0
        self.current_M_n = 0.0
        self.current_T_H = 0.0
        self.is_oscillating = False
        self.oscillation_amplitude = 0.0
        self.oscillation_frequency = 0.0
        
    def compute_energy_eigenvalue(self, J_values: dict, s_value: complex, loss: float):
        """
        Compute the energy eigenvalue E_n from current J-network state.
        
        E_n = (1/2) · [1 - J_0(2πA) · cos(2πs_0)]
        
        Where:
        - A = oscillation amplitude (from J-variance)
        - s_0 = oscillation center (from J-phase)
        """
        s = complex(s_value)
        
        # Extract oscillation amplitude from J-deviation from equilibrium (1,0)
        A_zg = abs(J_values['J_zg_complex'] - 1.0)
        A_zw = abs(J_values['J_zw_complex'] - 1.0) 
        A_gw = abs(J_values['J_gw_complex'] - 1.0)
        
        A_avg = (A_zg + A_zw + A_gw) / 3.0
        
        # Oscillation center from mean J-value
        s_0_zg = np.angle(J_values['J_zg_complex']) / (2 * np.pi)
        s_0_zw = np.angle(J_values['J_zw_complex']) / (2 * np.pi)
        s_0_gw = np.angle(J_values['J_gw_complex']) / (2 * np.pi)
        
        s_0_avg = (s_0_zg + s_0_zw + s_0_gw) / 3.0
        
        # Energy eigenvalue via Bessel formula
        if A_avg < 0.01:
            # Small oscillation limit
            J_0_approx = 1.0 - (2 * np.pi * A_avg)**2 / 4.0
        else:
            J_0_approx = bessel_j0(2 * np.pi * A_avg)
        
        E_n = 0.5 * (1.0 - J_0_approx * np.cos(2 * np.pi * s_0_avg))
        
        # Alternative: direct loss-based energy
        E_n_loss = loss
        
        # Combine both methods (weighted average)
        E_n_combined = 0.6 * E_n + 0.4 * E_n_loss
        
        return max(0.0, min(1.0, E_n_combined))  # Clamp to [0, 1]
    
    def compute_black_hole_mass(self, E_n: float, J_product: float):
        """
        Convert energy eigenvalue to physical black hole mass.
        
        M_n = E_n · m_P / √(J_ζΓ · J_ζW · J_ΓW)
        """
        if J_product <= 0:
            J_product = 1.0
        
        M_n = E_n * self.m_P / np.sqrt(J_product)
        return M_n
    
    def compute_hawking_temperature(self, E_n: float, M_n: float):
        """
        Compute Hawking temperature for this energy eigenvalue.
        
        T_H = (E_n · ħ · c³) / (8π · G · M_n · k_B)
        """
        if M_n <= 0:
            return float('inf')
        
        T_H = (E_n * self.hbar * self.c**3) / (8 * np.pi * self.G * M_n * self.k_B)
        return T_H
    
    def update(self, J_values: dict, s_value: complex, loss: float, d_value: float):
        """
        Update monitor with current state.
        """
        self.time_step += 1
        
        # Compute energy eigenvalue
        E_n = self.compute_energy_eigenvalue(J_values, s_value, loss)
        
        # Compute BH mass
        M_n = self.compute_black_hole_mass(E_n, J_values['J_product'])
        
        # Compute Hawking temperature
        T_H = self.compute_hawking_temperature(E_n, M_n)
        
        # Store values
        self.current_E_n = E_n
        self.current_M_n = M_n
        self.current_T_H = T_H
        
        self.energy_history.append(E_n)
        self.mass_history.append(M_n)
        self.temp_history.append(T_H)
        self.J_zg_history.append(J_values['J_zg'])
        self.J_zw_history.append(J_values['J_zw'])
        self.J_gw_history.append(J_values['J_gw'])
        self.s_trajectory.append(s_value)
        self.loss_history.append(loss)
        
        # Detect oscillation
        self._detect_oscillation()
        
        return {
            'E_n': E_n,
            'M_n_kg': M_n,
            'M_n_solar': M_n / 1.989e30,
            'T_H_K': T_H,
            'J_zg': J_values['J_zg'],
            'J_zw': J_values['J_zw'],
            'J_gw': J_values['J_gw'],
            'd_value': d_value
        }
    
    def _detect_oscillation(self):
        """Detect if J-values are oscillating using autocorrelation."""
        if len(self.energy_history) < 20:
            self.is_oscillating = False
            return
        
        # Use energy history for oscillation detection
        energies = np.array(self.energy_history[-50:])
        
        # Remove trend
        energies = energies - np.mean(energies)
        
        if np.std(energies) < 1e-6:
            self.is_oscillating = False
            self.oscillation_amplitude = 0.0
            return
        
        # Compute autocorrelation
        autocorr = np.correlate(energies, energies, mode='full')
        autocorr = autocorr[len(autocorr)//2:]
        autocorr = autocorr / autocorr[0]
        
        # Find first zero crossing (period)
        zero_crossings = np.where(np.diff(np.sign(autocorr[1:5])))[0]
        if len(zero_crossings) > 0:
            period = zero_crossings[0] + 1
            self.oscillation_frequency = 2 * np.pi / period if period > 0 else 0.0
            self.is_oscillating = True
        else:
            self.is_oscillating = False
        
        # Amplitude from variance
        self.oscillation_amplitude = np.std(energies)
    
    def get_status_string(self) -> str:
        """Get formatted status string."""
        status = f"""
╔══════════════════════════════════════════════════════════════════════╗
║           BLACK HOLE ENERGY MONITOR: {self.name:<20}           ║
╠══════════════════════════════════════════════════════════════════════╣
║                                                                      ║
║   TIME STEP: {self.time_step:<5}                                             ║
║                                                                      ║
║   ─── ENERGY EIGENVALUE ───────────────────────────────────────────  ║
║   E_n = {self.current_E_n:.6f}                                              ║
║   (0 = collapsed, 0.5 = max oscillation, 1 = impossible)            ║
║                                                                      ║
║   ─── BLACK HOLE OBSERVABLES ───────────────────────────────────────  ║
║   Mass:  M_n = {self.current_M_n:.6e} kg                              ║
║   Mass:  M_n = {self.current_M_n / 1.989e30:.6e} M_sun                    ║
║   Temp:  T_H = {self.current_T_H:.6e} K                              ║
║                                                                      ║
║   ─── J-NETWORK COUPLINGS ──────────────────────────────────────────  ║
║   J_ζΓ = {self.J_zg_history[-1] if self.J_zg_history else 0:.6f}                                      ║
║   J_ζW = {self.J_zw_history[-1] if self.J_zw_history else 0:.6f}                                      ║
║   J_ΓW = {self.J_gw_history[-1] if self.J_gw_history else 0:.6f}                                      ║
║                                                                      ║
║   ─── OSCILLATION STATUS ───────────────────────────────────────────  ║
║   Oscillating: {'YES' if self.is_oscillating else 'NO'}                                           ║
║   Amplitude:  {self.oscillation_amplitude:.6f}                                      ║
║   Frequency:  {self.oscillation_frequency:.6f} rad/step                       ║
║                                                                      ║
╚══════════════════════════════════════════════════════════════════════╝
"""
        return status
    
    def plot_energy_dynamics(self, save_path: str = None):
        """Plot energy eigenvalue dynamics."""
        if len(self.energy_history) < 2:
            return
        
        fig, axes = plt.subplots(3, 2, figsize=(14, 10))
        fig.suptitle(f'Black Hole Energy Monitor: {self.name}', fontsize=14, fontweight='bold')
        
        t = np.arange(len(self.energy_history))
        
        # 1. Energy eigenvalue over time
        ax1 = axes[0, 0]
        ax1.plot(t, self.energy_history, 'b-', linewidth=1.5, label='E_n(t)')
        ax1.axhline(y=0.5, color='r', linestyle='--', alpha=0.5, label='Max oscillation')
        ax1.axhline(y=0.0, color='gray', linestyle=':', alpha=0.5, label='Collapse')
        ax1.fill_between(t, 0, self.energy_history, alpha=0.3)
        ax1.set_xlabel('Time Step')
        ax1.set_ylabel('Energy Eigenvalue E_n')
        ax1.set_title('Energy Eigenvalue Dynamics')
        ax1.legend()
        ax1.grid(True, alpha=0.3)
        
        # 2. J-couplings over time
        ax2 = axes[0, 1]
        ax2.plot(t, self.J_zg_history, 'r-', linewidth=1.5, label='J_ζΓ', alpha=0.8)
        ax2.plot(t, self.J_zw_history, 'g-', linewidth=1.5, label='J_ζW', alpha=0.8)
        ax2.plot(t, self.J_gw_history, 'b-', linewidth=1.5, label='J_ΓW', alpha=0.8)
        ax2.axhline(y=1.0, color='k', linestyle='--', alpha=0.5, label='Equilibrium')
        ax2.set_xlabel('Time Step')
        ax2.set_ylabel('Coupling Strength')
        ax2.set_title('ζ-Γ-W Network Couplings')
        ax2.legend()
        ax2.grid(True, alpha=0.3)
        
        # 3. Black hole mass
        ax3 = axes[1, 0]
        mass_solar = np.array(self.mass_history) / 1.989e30
        ax3.semilogy(t, mass_solar, 'purple', linewidth=1.5)
        ax3.set_xlabel('Time Step')
        ax3.set_ylabel('Mass (M_sun)')
        ax3.set_title('Black Hole Mass from Energy Eigenvalue')
        ax3.grid(True, alpha=0.3)
        
        # 4. Hawking temperature
        ax4 = axes[1, 1]
        temps = np.array(self.temp_history)
        temps = np.clip(temps, 1e-100, 1e100)  # Clip for log plot
        ax4.semilogy(t, temps, 'orange', linewidth=1.5)
        ax4.set_xlabel('Time Step')
        ax4.set_ylabel('Temperature (K)')
        ax4.set_title('Hawking Temperature')
        ax4.grid(True, alpha=0.3)
        
        # 5. Loss function
        ax5 = axes[2, 0]
        ax5.semilogy(t, np.array(self.loss_history) + 1e-100, 'k-', linewidth=1.5)
        ax5.set_xlabel('Time Step')
        ax5.set_ylabel('Loss ℒ')
        ax5.set_title('CCT Loss Function')
        ax5.grid(True, alpha=0.3)
        
        # 6. Phase portrait (J_zg vs J_zw)
        ax6 = axes[2, 1]
        J_zg_complex = np.array([compute_J_couplings(s, 1.0)['J_zg_complex'] for s in self.s_trajectory[-100:]])
        ax6.scatter(J_zg_complex.real, J_zg_complex.imag, c=t[-100:], cmap='viridis', s=20, alpha=0.7)
        ax6.scatter([1.0], [0.0], c='red', s=100, marker='x', label='Equilibrium')
        ax6.set_xlabel('Re(J_ζΓ)')
        ax6.set_ylabel('Im(J_ζΓ)')
        ax6.set_title('J_ζΓ Phase Portrait')
        ax6.legend()
        ax6.grid(True, alpha=0.3)
        ax6.set_aspect('equal')
        
        plt.tight_layout()
        
        if save_path:
            plt.savefig(save_path, dpi=150, bbox_inches='tight')
            print(f"Saved plot to {save_path}")
        
        plt.show()
    
    def print_summary(self):
        """Print a summary of the energy dynamics."""
        if not self.energy_history:
            print("No data collected.")
            return
        
        print(f"\n{'='*70}")
        print(f"  BLACK HOLE ENERGY MONITOR SUMMARY: {self.name}")
        print(f"{'='*70}")
        
        print(f"\n  Configuration:")
        print(f"    Time steps: {len(self.energy_history)}")
        print(f"    Oscillating: {'YES' if self.is_oscillating else 'NO'}")
        
        print(f"\n  Energy Statistics:")
        print(f"    E_n (final):   {self.energy_history[-1]:.6f}")
        print(f"    E_n (mean):    {np.mean(self.energy_history):.6f}")
        print(f"    E_n (std):     {np.std(self.energy_history):.6f}")
        print(f"    E_n (min):     {np.min(self.energy_history):.6f}")
        print(f"    E_n (max):     {np.max(self.energy_history):.6f}")
        
        print(f"\n  Black Hole Properties (final state):")
        print(f"    Mass: {self.mass_history[-1]:.6e} kg = {self.mass_history[-1]/1.989e30:.6e} M_sun")
        print(f"    Temp: {self.temp_history[-1]:.6e} K")
        
        print(f"\n  J-Couplings (final state):")
        print(f"    J_ζΓ: {self.J_zg_history[-1]:.6f}")
        print(f"    J_ζW: {self.J_zw_history[-1]:.6f}")
        print(f"    J_ΓW: {self.J_gw_history[-1]:.6f}")
        
        print(f"\n  Oscillation Analysis:")
        print(f"    Amplitude: {self.oscillation_amplitude:.6f}")
        print(f"    Frequency: {self.oscillation_frequency:.6f} rad/step")
        
        print(f"\n{'='*70}\n")


# ============================================================================
# PART 3: CCT GRADIENT FACTORIZATION WITH ENERGY MONITORING
# ============================================================================

def factor_cct_gradient_with_energy_monitor(
    c: int,
    lr: float = 0.1,
    max_iter: int = 2000,
    tol: float = 1e-8,
    num_restarts: int = 5,
    device: str = 'cpu',
    verbose: bool = True,
    track_energy: bool = True
) -> tuple:
    """
    Factor a semiprime c = p*q using gradient descent on the gap d,
    with black hole energy monitoring.
    
    Args:
        c: integer to factor (product of two primes)
        lr: learning rate for Adam optimizer
        max_iter: maximum gradient steps per restart
        tol: loss tolerance for convergence
        num_restarts: number of random initialisations
        device: 'cpu' or 'cuda'
        verbose: print progress
        track_energy: enable energy monitoring
    
    Returns:
        (p, q, energy_monitor) such that p <= q and p*q == c
    """
    c_t = torch.tensor(float(c), device=device, requires_grad=False)
    
    # Loss function
    def loss_fn(d):
        s = torch.sqrt(d**2 + 4.0 * c_t)
        loss = torch.sin(torch.pi * s) ** 2
        return loss, s
    
    # Initialize energy monitor
    energy_monitor = BlackHoleEnergyMonitor(name=f"BH-c={c}") if track_energy else None
    
    best_d = None
    best_loss = float('inf')
    best_s = None
    
    for restart in range(num_restarts):
        if verbose:
            print(f"  Restart {restart + 1}/{num_restarts}...", end=' ')
        
        # Initialize d randomly
        d = torch.randn(1, device=device) * 10.0
        d = d.detach().requires_grad_()
        optimizer = torch.optim.Adam([d], lr=lr)
        
        for step in range(max_iter):
            optimizer.zero_grad()
            loss, s_val = loss_fn(d)
            loss.backward()
            optimizer.step()
            
            # Clamp d to positive
            with torch.no_grad():
                d.clamp_(min=0.0)
            
            # Update energy monitor every 10 steps
            if track_energy and step % 10 == 0:
                d_np = d.item()
                s_np = s_val.item()
                loss_np = loss.item()
                
                J_values = compute_J_couplings(complex(s_np), M_normalized=1.0)
                energy_monitor.update(J_values, complex(s_np), loss_np, d_np)
            
            # Track best
            if loss.item() < best_loss:
                best_loss = loss.item()
                best_d = d.item()
                best_s = s_val.item()
            
            if loss.item() < tol:
                break
        
        if verbose:
            print(f"loss={loss.item():.2e}")
        
        if best_loss < tol:
            break
    
    # Post-processing: round and verify
    d_candidate = round(best_d)
    s_sq = d_candidate * d_candidate + 4 * c
    s = int(round(math.sqrt(s_sq)))
    
    if s * s != s_sq:
        # Fallback to Fermat search
        s = math.isqrt(4 * c) + 1
        while True:
            diff = s * s - 4 * c
            if diff < 0:
                s += 1
                continue
            d2 = math.isqrt(diff)
            if d2 * d2 == diff:
                d_candidate = d2
                break
            s += 1
    
    # Recover factors
    a = s - d_candidate
    b = s + d_candidate
    p = a // 2
    q = b // 2
    
    if p > q:
        p, q = q, p
    
    if p * q == c:
        return p, q, energy_monitor
    else:
        raise ValueError(f"Gradient descent converged to invalid pair: ({p}, {q})")


# ============================================================================
# PART 4: MAIN DEMONSTRATION
# ============================================================================

if __name__ == "__main__":
    print("\n" + "="*70)
    print("  CCT AUTOMATON WITH BLACK HOLE ENERGY MONITOR")
    print("="*70 + "\n")
    
    print("Testing factorization with energy eigenvalue tracking...")
    print("The oscillating J-network solutions represent black hole energy modes.\n")
    
    # Test cases: vary in size to see different oscillation behaviors
    test_cases = [
        293579,      # Small (from original paper)
        1339573,     # Medium
        53817427,    # Larger
        # Add more for variety
    ]
    
    results = []
    
    for i, c in enumerate(test_cases):
        print(f"\n{'─'*70}")
        print(f"TEST {i+1}: c = {c}")
        print(f"{'─'*70}")
        
        try:
            p, q, energy_monitor = factor_cct_gradient_with_energy_monitor(
                c,
                num_restarts=3,
                max_iter=1000,
                verbose=True,
                track_energy=True
            )
            
            print(f"\n  ✓ Factored: {p} × {q} = {p*q}")
            
            if energy_monitor:
                # Print status
                print(energy_monitor.get_status_string())
                
                # Save plot
                plot_path = f"energy_monitor_test_{i+1}.png"
                energy_monitor.plot_energy_dynamics(save_path=plot_path)
                
                # Print summary
                energy_monitor.print_summary()
                
                results.append({
                    'c': c,
                    'p': p,
                    'q': q,
                    'monitor': energy_monitor
                })
                
        except Exception as e:
            print(f"  ✗ Error: {e}")
    
    # Summary comparison
    if results:
        print("\n" + "="*70)
        print("  SUMMARY: ENERGY EIGENVALUES ACROSS TEST CASES")
        print("="*70)
        
        print(f"\n{'c':<15} {'p':<12} {'q':<12} {'E_n':<12} {'M_n (M_sun)':<15} {'Oscillating':<12}")
        print("-" * 70)
        
        for r in results:
            mon = r['monitor']
            E_n = mon.energy_history[-1] if mon.energy_history else 0
            M_n = mon.mass_history[-1] / 1.989e30 if mon.mass_history else 0
            osc = 'YES' if mon.is_oscillating else 'NO'
            
            print(f"{r['c']:<15} {r['p']:<12} {r['q']:<12} {E_n:<12.6f} {M_n:<15.6e} {osc:<12}")
        
        print("\n" + "="*70)
        print("\n  INTERPRETATION:")
        print("  ─────────────────────────────────────────────────────────────────")
        print("  • E_n = 0: Black hole collapsed to singularity (factors found)")
        print("  • E_n ≈ 0.5: Maximum oscillating energy mode")
        print("  • M_n: Mass of mathematical black hole from energy eigenvalue")
        print("  • Oscillating J-values: The 'breathing' of the ζ-Γ-W network")
        print("="*70 + "\n")