"""omega_solver.py
Solve for the Mersenne-domain field tension Omega for each known Mersenne prime.

Model claim: all constants in a domain share a near-constant Omega.
For Mersenne primes, with base=2, n+beta = PLN omega:

  D_n(omega, r) = sqrt(phi * F_omega * 2^omega * P_omega * Omega) * r^k = target_value

Solving for Omega:
  Omega = target^2 / (phi * F_omega * 2^omega * P_omega * r^(2k))

If Mersenne primes form a coherent domain, the derived log(Omega) values cluster tightly.
Compare against random non-Mersenne-prime exponents in the same scale range.
"""
import math
import random
import statistics

PHI = (1 + math.sqrt(5)) / 2

# ================================================================
# FIRST 50 PRIMES for cyclic P_n indexing
# ================================================================
def _sieve(n):
    s = [True] * (n + 1)
    s[0] = s[1] = False
    for i in range(2, int(n**0.5) + 1):
        if s[i]:
            for j in range(i*i, n+1, i):
                s[j] = False
    return [i for i in range(2, n+1) if s[i]]

PRIMES_50 = _sieve(250)[:50]

# ================================================================
# GENERALIZED FIBONACCI (Binet — non-integer x)
# ================================================================
def gen_fib(x):
    """F_x via Binet's formula for real x."""
    pos = PHI ** x
    try:
        neg = ((-PHI) ** (-x)) * math.cos(math.pi * x)
    except Exception:
        neg = 0.0
    return (pos - neg) / math.sqrt(5)

def cyclic_prime(x):
    idx = (int(math.floor(x)) + 50) % 50
    return PRIMES_50[idx]

# ================================================================
# PLN COORDINATES
# ================================================================
def compute_pln(p):
    omega = math.log(math.log(p * math.log(2)) / math.log(PHI)) / math.log(PHI)
    dw = omega % 1.0
    r_raw = dw / (1.0 + dw)
    r = r_raw
    for _ in range(8):
        r = r_raw * (1.0 + math.exp(-r * r / (2.0 * PHI)))
    D = int(omega % 8) + 1
    return omega, r, D

# ================================================================
# OMEGA SOLVER — four k-mode variants
# ================================================================
def solve_domain_omega(p, base=2, k_mode='phi01', target='p'):
    """
    Solve for the domain field tension Omega given exponent p.

    k_mode options:
      'phi01'    : k = phi^0.10 ~ 1.0494 (universal cluster from MICRO_TUNING_SUMMARY)
      'pln_r'    : k = PLN orbit radius r  (domain-specific, ~ 0.3-0.6)
      'pln_rk'   : k = (omega+1)/8        (PLN D_n convention)
      'unity'    : k = 1.0               (baseline)

    target options:
      'p'        : D_n = p  (the Mersenne exponent)
      'logp'     : D_n = log2(p)  (bit-depth of p)
      'omega'    : D_n = omega   (PLN coord; dimensionless)
    """
    omega, r, D = compute_pln(p)

    F_w = gen_fib(omega)
    P_w = cyclic_prime(omega)
    base_w = base ** omega   # 2^7.56 ~ 188 when omega~7.56

    if target == 'p':
        T = float(p)
    elif target == 'logp':
        T = math.log2(p)
    elif target == 'omega':
        T = omega
    else:
        raise ValueError(f"unknown target {target}")

    if k_mode == 'phi01':
        k = PHI ** 0.10
    elif k_mode == 'pln_r':
        k = r
    elif k_mode == 'pln_rk':
        k = (omega + 1) / 8.0
    elif k_mode == 'unity':
        k = 1.0
    else:
        raise ValueError(f"unknown k_mode {k_mode}")

    r_pow = r ** (2.0 * k)

    denom = PHI * abs(F_w) * base_w * P_w * r_pow
    if denom < 1e-300:
        return None

    Omega = (T * T) / denom
    return Omega, omega, r, k

# ================================================================
# SIMPLE PRIMALITY CHECK
# ================================================================
def is_probably_prime(n):
    if n < 2: return False
    if n == 2: return True
    if n % 2 == 0: return False
    for d in range(3, min(int(n**0.5) + 1, 1 << 18), 2):
        if n % d == 0:
            return False
    return True

# ================================================================
# ALL 52 KNOWN MERSENNE EXPONENTS
# ================================================================
ALL_P = [
    2,3,5,7,13,17,19,31,61,89,107,127,521,607,1279,2203,2281,3217,4253,
    4423,9689,9941,11213,19937,21701,23209,44497,86243,110503,132049,
    216091,756839,859433,1257787,1398269,2976221,3021377,6972593,13466917,
    20996011,24036583,25964951,30402457,32582657,37156667,42643801,43112609,
    57885161,74207281,77232917,82589933,136279841
]

# ================================================================
# MAIN ANALYSIS
# ================================================================
def analyze(base=2, k_mode='phi01', target='p', n_random=1000):
    print(f"\n{'='*65}")
    print(f"  base={base}, k_mode={k_mode}, target={target}")
    print(f"{'='*65}")

    # --- Known Mersenne primes ---
    mersenne_omegas = []
    mersenne_log_omegas = []

    print(f"\n  {'p':>12}  {'omega':>7}  {'r':>6}  {'k':>6}  {'log10(Omega)':>14}")
    print("  " + "-"*55)

    for p in ALL_P:
        res = solve_domain_omega(p, base=base, k_mode=k_mode, target=target)
        if res is None:
            print(f"  {p:>12d}  FAILED")
            continue
        Omega, omega, r, k = res
        if Omega <= 0:
            continue
        log_Om = math.log10(Omega)
        mersenne_omegas.append(Omega)
        mersenne_log_omegas.append(log_Om)
        print(f"  {p:>12d}  {omega:7.4f}  {r:6.4f}  {k:6.4f}  {log_Om:14.6f}")

    if len(mersenne_log_omegas) < 3:
        print("  Not enough data")
        return

    mean_log = statistics.mean(mersenne_log_omegas)
    stdev_log = statistics.stdev(mersenne_log_omegas)
    cv_log = abs(stdev_log / mean_log) if mean_log != 0 else float('inf')

    print(f"\n  Mersenne known primes (n={len(mersenne_log_omegas)}):")
    print(f"    mean(log10 Omega)  = {mean_log:.6f}")
    print(f"    stdev(log10 Omega) = {stdev_log:.6f}")
    print(f"    c.v.(log10)        = {cv_log:.6f}")
    print(f"    range              = {min(mersenne_log_omegas):.4f}  to  {max(mersenne_log_omegas):.4f}")

    # --- Random prime exponents ---
    random.seed(42)
    random_log_omegas = []
    attempts = 0
    p_set = set(ALL_P)

    # Sample proportionally across the exponent range
    lo, hi = min(ALL_P), max(ALL_P)
    while len(random_log_omegas) < n_random and attempts < n_random * 40:
        attempts += 1
        cand_raw = random.randint(lo, hi)
        cand = cand_raw | 1  # make odd
        if cand in p_set:
            continue
        if not is_probably_prime(cand):
            continue
        res = solve_domain_omega(cand, base=base, k_mode=k_mode, target=target)
        if res is None:
            continue
        Omega, _, _, _ = res
        if Omega > 0:
            random_log_omegas.append(math.log10(Omega))

    if len(random_log_omegas) < 10:
        print("  (not enough random samples)")
        return

    rand_mean = statistics.mean(random_log_omegas)
    rand_stdev = statistics.stdev(random_log_omegas)
    rand_cv = abs(rand_stdev / rand_mean) if rand_mean != 0 else float('inf')

    print(f"\n  Random prime exponents (n={len(random_log_omegas)}):")
    print(f"    mean(log10 Omega)  = {rand_mean:.6f}")
    print(f"    stdev(log10 Omega) = {rand_stdev:.6f}")
    print(f"    c.v.(log10)        = {rand_cv:.6f}")
    print(f"    range              = {min(random_log_omegas):.4f}  to  {max(random_log_omegas):.4f}")

    # --- Discrimination metric ---
    # Cohen's d between Mersenne log-Omega and random log-Omega
    pooled_std = math.sqrt((stdev_log**2 + rand_stdev**2) / 2)
    cohens_d = abs(mean_log - rand_mean) / pooled_std if pooled_std > 0 else 0.0

    # Ratio of variances (Levene-like)
    var_ratio = stdev_log / rand_stdev if rand_stdev > 0 else float('inf')

    # F-test: Mersenne sigma / random sigma  (< 1 means Mersenne tighter)
    print(f"\n  DISCRIMINATION METRICS:")
    print(f"    stdev ratio (Mersenne/Random) = {var_ratio:.4f}  (< 1.0 = tighter cluster)")
    print(f"    Cohen's d (mean separation)   = {cohens_d:.4f}  (> 1.0 = strong separation)")

    # What fraction of random primes fall within 2 sigma of Mersenne mean?
    lo_bound = mean_log - 2 * stdev_log
    hi_bound = mean_log + 2 * stdev_log
    inside = sum(1 for x in random_log_omegas if lo_bound <= x <= hi_bound) / len(random_log_omegas)
    print(f"    Random primes in Mersenne 2-sigma band: {100*inside:.1f}%  (< 5% = well separated)")

    return {
        'mersenne_mean': mean_log, 'mersenne_stdev': stdev_log,
        'random_mean': rand_mean, 'random_stdev': rand_stdev,
        'var_ratio': var_ratio, 'cohens_d': cohens_d, 'inside_pct': inside
    }


if __name__ == '__main__':
    import time
    t0 = time.time()

    print("MERSENNE-DOMAIN OMEGA CLUSTERING ANALYSIS")
    print("Based on BigG + Fudge10 model: phi is the only universal constant")
    print("All others (including Omega = field tension) are domain-emergent")

    results = []

    # Sweep all combinations to find the most discriminating configuration
    for base in [2]:
        for k_mode in ['phi01', 'pln_r', 'pln_rk', 'unity']:
            for target in ['p', 'logp', 'omega']:
                res = analyze(base=base, k_mode=k_mode, target=target, n_random=800)
                if res:
                    results.append((base, k_mode, target, res))

    # Summary table
    print("\n\n" + "="*75)
    print("SUMMARY: All configurations (ranked by stdev ratio — lower = tighter)")
    print("="*75)
    print(f"  {'base':>4}  {'k_mode':>10}  {'target':>6}  {'Mers.std':>9}  {'Rand.std':>9}  {'std_ratio':>9}  {'d':>6}  {'inside%':>8}")
    print("  " + "-"*72)

    results.sort(key=lambda x: x[3]['var_ratio'])
    for base, k_mode, target, r in results:
        print(f"  {base:>4}  {k_mode:>10}  {target:>6}  "
              f"{r['mersenne_stdev']:9.4f}  {r['random_stdev']:9.4f}  "
              f"{r['var_ratio']:9.4f}  {r['cohens_d']:6.3f}  "
              f"{100*r['inside_pct']:8.1f}%")

    print(f"\nTotal time: {time.time()-t0:.1f}s")

    # Best config deep-dive
    best = results[0]
    base, k_mode, target, r = best
    print(f"\nBest configuration: base={base}, k_mode={k_mode}, target={target}")
    print(f"  This means: Mersenne primes cluster {1/r['var_ratio']:.1f}x tighter in Omega-space")
    if r['var_ratio'] < 0.5:
        print("  => STRONG DOMAIN SIGNAL: Mersenne primes occupy a coherent Omega band")
    elif r['var_ratio'] < 0.8:
        print("  => MODERATE SIGNAL: some Omega clustering")
    else:
        print("  => WEAK/NONE: Omega does not cleanly separate Mersenne from random primes")
