#!/usr/bin/env python3
"""
Computation 48 -- Refined bridge closure verification at higher k
==================================================================
Computations 45-47 demonstrated:
  - 1-monomial bridge: ratio = 2 exactly at k = 1 (Comp 45)
  - 2-monomial superposition: gap < 1e-3 at even k = 2, 4, 6, 8 (Comp 46)
  - 3-monomial superposition: gap < 1e-4 at odd k = 3, 7 (Comp 47)

Two gaps in the verification:
  - k = 5: Comp 47 search converged to a local minimum (alpha = beta = 0)
    and didn't find closure with the 3-monomial family
  - k > 8: not tested

This script:
  (a) re-runs k = 5 with a more careful multi-restart search to close
      the local-minimum artifact
  (b) extends the 3-monomial framework to k = 9, 10, 11, 12 to verify
      that the refined-bridge closure pattern continues at higher
      weight classes
"""
import math
import numpy as np


def F_sup_3mon(m, alpha, beta, n_theta=24, n_phi=14):
    best = 0.0
    for i in range(n_theta + 1):
        theta = (i / n_theta) * (np.pi / 2)
        r = np.cos(theta)
        s = np.sin(theta)
        for j in range(n_phi + 1):
            phi0 = (j / n_phi) * (2 * np.pi)
            for k in range(n_phi + 1):
                phi1 = (k / n_phi) * (2 * np.pi)
                z0 = r * np.exp(1j * phi0)
                z1 = s * np.exp(1j * phi1)
                u = z0 * z1
                val = abs(u**(m - 1) + beta * u**m + alpha * u**(m + 1))
                if val > best:
                    best = val
    return best


def M_sup_3mon(m, alpha, beta, n_theta=20, n_phi=12):
    """
    F = u^{m-1} + beta u^m + alpha u^{m+1}, u = z_0 z_1.
    J_+ F = z_0 d/dz_1 F
          = (m-1) z_0^m z_1^{m-2} + beta * m * z_0^{m+1} z_1^{m-1}
            + alpha * (m+1) * z_0^{m+2} z_1^m
    """
    best = 0.0
    for i in range(n_theta + 1):
        theta = (i / n_theta) * (np.pi / 2)
        r = np.cos(theta)
        s = np.sin(theta)
        for j in range(n_phi + 1):
            phi0 = (j / n_phi) * (2 * np.pi)
            for k in range(n_phi + 1):
                phi1 = (k / n_phi) * (2 * np.pi)
                z0 = r * np.exp(1j * phi0)
                z1 = s * np.exp(1j * phi1)
                if m >= 2:
                    t1 = (m - 1) * z0**m * z1**(m - 2)
                else:
                    t1 = 0
                t2 = beta * m * z0**(m + 1) * z1**(m - 1) if m >= 1 else 0
                t3 = alpha * (m + 1) * z0**(m + 2) * z1**m
                jp = t1 + t2 + t3
                if m >= 2:
                    s1 = (m - 1) * z0**(m - 2) * z1**m
                else:
                    s1 = 0
                s2 = beta * m * z0**(m - 1) * z1**(m + 1) if m >= 1 else 0
                s3 = alpha * (m + 1) * z0**m * z1**(m + 2)
                jm = s1 + s2 + s3
                jz_val = 0.0 + 0.0j
                naive = 0.5 * (abs(jp)**2 + abs(jm)**2) + abs(jz_val)**2
                X = 0.5 * (abs(jp)**2 - abs(jm)**2)
                Y_imag = (jz_val.conjugate() * jm).imag
                Y_sq = 4.0 * Y_imag**2
                op2 = naive + math.sqrt(X**2 + Y_sq)
                val = math.sqrt(op2)
                if val > best:
                    best = val
    return best


def search_with_restarts(k, m_range, restart_seeds):
    """Grid search with multiple starting points and refinement."""
    target = 2 * math.sqrt(k)
    best_overall = None
    for m in m_range:
        for (beta0, alpha0) in restart_seeds:
            # Coarse grid centered on each seed
            alpha_grid = [alpha0 + 0.5 * d for d in range(-6, 7)]
            beta_grid = [beta0 + 0.5 * d for d in range(-6, 7)]
            for alpha in alpha_grid:
                for beta in beta_grid:
                    sF = F_sup_3mon(m, alpha, beta)
                    if sF < 1e-9:
                        continue
                    sM = M_sup_3mon(m, alpha, beta)
                    r = sM / sF
                    gap = abs(r - target)
                    if best_overall is None or gap < best_overall[0]:
                        best_overall = (gap, m, alpha, beta, r, sF)
    if best_overall is None:
        return None
    # Refine the best
    gap, m, alpha_best, beta_best, ratio, sF = best_overall
    for scale in [0.2, 0.05, 0.01]:
        for d_alpha in range(-5, 6):
            for d_beta in range(-5, 6):
                a = alpha_best + scale * d_alpha
                b = beta_best + scale * d_beta
                sF = F_sup_3mon(m, a, b)
                if sF < 1e-9:
                    continue
                sM = M_sup_3mon(m, a, b)
                r = sM / sF
                g = abs(r - target)
                if g < gap:
                    gap = g
                    alpha_best = a
                    beta_best = b
                    ratio = r
    return (gap, m, alpha_best, beta_best, ratio)


def main():
    print("=" * 90)
    print("  Computation 48  --  Refined bridge closure verification at higher k")
    print("=" * 90)
    print()
    print("  Symbol: f(z) = (z_0 z_1)^{m-1} + beta (z_0 z_1)^m + alpha (z_0 z_1)^{m+1}")
    print()
    print(f"  {'k':>3}  {'2 sqrt(k)':>11}  {'m':>3}  {'beta':>10}  {'alpha':>10}"
          f"  {'ratio':>12}  {'gap':>12}")

    # Multiple seed points to escape local minima
    restart_seeds = [
        (0.0, 0.0),
        (-1.5, 1.5),
        (1.5, -1.5),
        (-2.0, 2.0),
        (2.0, -2.0),
        (-3.0, 3.0),
        (3.0, -3.0),
        (0.0, 2.0),
        (2.0, 0.0),
        (-2.0, 0.0),
        (0.0, -2.0),
    ]

    target_ks = [5, 9, 10, 11, 12]
    results = []
    for k in target_ks:
        m_range = list(range(max(1, int(math.sqrt(k)) - 1), int(math.sqrt(k)) + 4))
        result = search_with_restarts(k, m_range, restart_seeds)
        if result is None:
            continue
        gap, m, alpha, beta, ratio = result
        results.append((k, m, beta, alpha, ratio, gap))
        target = 2 * math.sqrt(k)
        print(f"  {k:>3}  {target:>11.4f}  {m:>3}  {beta:>10.4f}  {alpha:>10.4f}"
              f"  {ratio:>12.5f}  {gap:>12.5f}")

    print()
    print("=" * 90)
    print("  Verdict")
    print("=" * 90)
    print()
    max_gap = max(r[5] for r in results)
    print(f"  Max gap across k = 5, 9, 10, 11, 12:  {max_gap:.5f}")
    if max_gap < 1e-3:
        print()
        print("  => Refined bridge closure verified to gap < 1e-3 at all tested k.")
        print("     Combined with Comp 45-47 results at k = 1..8, the refined")
        print("     symmetric-monomial bridge closes L_comm uniformly across the")
        print("     weight classes k = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12.")
    elif max_gap < 1e-2:
        print(f"  => essentially exact closure (gap < 1%)")
    else:
        print(f"  => some k retain residual gap; investigate")


if __name__ == "__main__":
    main()