"""riemann_edge.py
Riemann critical edge analysis for Mersenne prime candidates.

In the PLN framework:
  omega(p)  = log_phi(log_phi(p * ln2))        -- PLN coordinate
  beta(p)   = frac(omega(p))                   -- fractional part in [0,1)

beta = 0.5 maps exactly to Re(s) = 1/2 on the Riemann critical LINE.
The Euler/phi/physical bridge:
  e^{i*pi} = -1  =>  e^{i*pi*beta} at beta=0.5 = e^{i*pi/2} = i
  This is the 90-degree rotation point: maximum phase-entropy in the phasor.

The local calibration Gaussian exp(-0.5*(delta_omega/BW)^2) with BW=0.5
is a spectral filter CENTRED ON THE CRITICAL LINE:
  inflection point at delta_omega = BW = 0.5 = Re(s) = 1/2.

Moving the prism (query_omega) shifts which primes are within BW of the
critical line -- the Riemann edge follows the prism.

Geometric progression of prime structure (Ulam spiral insight):
  dot    (0D): individual candidate p
  line   (1D): omega axis -- the PLN spiral
  delta  (2D): consecutive beta gaps -- quadratic acceleration
  square (2D): (omega, beta) plane -- full Riemann strip map
"""

import math, os

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

KP = [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]

# LOO ranks from best run (BW=0.5, 33.4% mean, contamination-cleaned)
LOO_RANK = {
    2:0.0, 3:2.5, 5:13.6, 7:8.3, 13:84.7, 17:82.5, 19:81.9, 31:89.8,
    61:80.5, 89:78.4, 107:70.4, 127:73.9,
    521:47.1, 607:39.0, 1279:51.9, 2203:21.7, 2281:24.2, 3217:33.8,
    4253:39.8, 4423:46.7,
    9689:63.8, 9941:64.2, 11213:67.2, 19937:69.0, 21701:68.8, 23209:66.4,
    44497:8.0, 86243:35.4, 110503:61.6, 132049:11.2,
    216091:31.4, 756839:32.4, 859433:34.4, 1257787:79.6, 1398269:87.6,
    2976221:66.6, 3021377:66.6, 6972593:69.4, 13466917:69.2,
    20996011:14.2, 24036583:4.8, 25964951:12.2, 30402457:27.0,
    32582657:34.0, 37156667:33.4, 42643801:35.8, 43112609:31.8,
    57885161:34.0, 74207281:34.8, 77232917:30.8, 82589933:32.4,
    136279841:33.8
}


def pln_omega(p):
    inner = math.log(p * LN2) / LN_PHI
    return math.log(inner) / LN_PHI


# =========================================================================
# 1. Known primes in (omega, beta) space with Riemann edge markers
# =========================================================================
print("=" * 78)
print("  KNOWN MERSENNE PRIMES: (omega, beta) MAPPED TO RIEMANN STRIP")
print("  Riemann critical line: beta = 0.5 = Re(s) = 1/2")
print("  BW = 0.5 calibration bandwidth = critical strip half-width")
print("=" * 78)
print()
print(f"  {'p':>14}  {'omega':>8}  {'beta':>8}  {'|b-0.5|':>8}  {'LOO%':>6}  note")
print("  " + "-" * 72)

for p in KP:
    om = pln_omega(p)
    b  = om - math.floor(om)
    d  = abs(b - 0.5)
    loo = LOO_RANK.get(p, 50.0)
    note = "*** RIEMANN LINE  " if d < 0.02 else \
           "**  near edge     " if d < 0.08 else \
           "*   approaching   " if d < 0.15 else ""
    print(f"  {p:>14}  {om:8.4f}  {b:8.4f}  {d:8.4f}  {loo:6.1f}  {note}")

# =========================================================================
# 2. Correlation between |beta-0.5| and LOO rank
# =========================================================================
print()
print("=" * 78)
print("  RIEMANN EDGE CORRELATION: |beta - 0.5| vs LOO rank")
print("  Test: primes near beta=0.5 should have lower LOO rank (better signal)")
print("=" * 78)

items = [(abs((pln_omega(p) - math.floor(pln_omega(p))) - 0.5),
          LOO_RANK.get(p, 50.0), p) for p in KP]
items.sort()

# Quantile analysis
n = len(items)
q1, q2, q3 = n//4, n//2, 3*n//4
d_vals  = [x[0] for x in items]
r_vals  = [x[1] for x in items]

def mean(lst): return sum(lst)/len(lst)

print(f"  Quartile analysis (sorted by |beta-0.5|):")
print(f"  Q1 (closest to edge) n={q1}:  |b-0.5| mean={mean(d_vals[:q1]):.3f}  LOO mean={mean(r_vals[:q1]):.1f}%")
print(f"  Q2                   n={q2-q1}:  |b-0.5| mean={mean(d_vals[q1:q2]):.3f}  LOO mean={mean(r_vals[q1:q2]):.1f}%")
print(f"  Q3                   n={q3-q2}:  |b-0.5| mean={mean(d_vals[q2:q3]):.3f}  LOO mean={mean(r_vals[q2:q3]):.1f}%")
print(f"  Q4 (farthest)        n={n-q3}:  |b-0.5| mean={mean(d_vals[q3:]):.3f}  LOO mean={mean(r_vals[q3:]):.1f}%")

# Correlation coefficient
mx = mean(d_vals); my = mean(r_vals)
num   = sum((d_vals[i]-mx)*(r_vals[i]-my) for i in range(n))
denom = math.sqrt(sum((x-mx)**2 for x in d_vals)*sum((y-my)**2 for y in r_vals))
corr  = num/denom if denom > 0 else 0.0
print(f"\n  Pearson r(|beta-0.5|, LOO_rank) = {corr:+.4f}")
print(f"  {'POSITIVE correlation: closer to Riemann edge -> better LOO rank' if corr > 0.2 else 'weak or negative' if abs(corr) < 0.2 else 'NEGATIVE'}")

# =========================================================================
# 3. Where is the Riemann edge relative to the p=136M prism?
# =========================================================================
print()
print("=" * 78)
print("  PRISM AT p=136279841: LOCAL RIEMANN EDGE MAP")
print("=" * 78)

om_prism = pln_omega(136279841)
b_prism  = om_prism - math.floor(om_prism)
print(f"  Prism omega = {om_prism:.6f}   beta = {b_prism:.6f}")
print(f"  Distance from Riemann line: |beta - 0.5| = {abs(b_prism-0.5):.4f}")
print(f"  BW = 0.5: primes within omega +/- 1.0 are in calibration neighborhood")
print()
print(f"  Known primes within BW=1.0 of prism (|omega_i - omega_prism| < 1.0):")
print(f"  {'p':>14}  {'omega':>8}  {'beta':>8}  {'|b-0.5|':>8}  {'d_omega':>8}")
print("  " + "-" * 56)
neighbors = []
for p in KP:
    om = pln_omega(p)
    d_om = abs(om - om_prism)
    if d_om < 1.0 and p != 136279841:
        b = om - math.floor(om)
        neighbors.append((d_om, p, om, b))
neighbors.sort()
for d_om, p, om, b in neighbors:
    print(f"  {p:>14}  {om:8.4f}  {b:8.4f}  {abs(b-0.5):8.4f}  {d_om:+8.4f}")

# =========================================================================
# 4. GPU candidate range: find the p where beta = 0.5
# =========================================================================
print()
print("=" * 78)
print("  RIEMANN CRITICAL LINE IN GPU CANDIDATE RANGE")
print("  Find p such that frac(omega(p)) = 0.5 near p=136M")
print("=" * 78)

# Solve: frac(omega(p)) = 0.5
# omega(p) = n + 0.5  for n = floor(omega(p)) = 7
# phi^(7.5) = log_phi(p * ln2)
# p * ln2 = phi^(phi^7.5)
om_target = math.floor(om_prism) + 0.5   # = 7.5
phi_val = PHI ** om_target
p_crit = (PHI ** phi_val) / LN2
print(f"  omega = {om_target} (n=7, beta=0.5):  p_critical = {p_crit:.2e}  ({int(p_crit):,})")
print()

# Also find the next integer + 0.5 above prism
for n_test in range(7, 10):
    phi_val = PHI ** (n_test + 0.5)
    p_c = (PHI ** phi_val) / LN2
    om_c = pln_omega(int(p_c))
    b_c  = om_c - math.floor(om_c)
    print(f"  n={n_test}, beta=0.5:  p = {p_c:.6e}  ({int(p_c):,})")

# Now scan GPU candidates from file if it exists
# Try both candidate files; ranked_candidates has consensus scores
for cand_file in ("ranked_candidates.txt", "priority_gpu_final.txt"):
    if os.path.exists(cand_file):
        cands = []
        with open(cand_file) as f:
            for line in f:
                line = line.strip()
                if not line or line.startswith("#"):
                    continue
                parts = line.split()
                # ranked_candidates: rank p PLN PLN% BigG consensus
                # priority_gpu_final: rank p BigG source ...
                # In both cases column 1 (index 1) is the prime p
                try:
                    p = int(parts[1])
                    if p < 100000:          # skip header artefacts
                        continue
                    om = pln_omega(p)
                    b  = om - math.floor(om)
                    cands.append((abs(b - 0.5), b, om, p))
                except (ValueError, IndexError):
                    pass
        if cands:
            cands.sort()
            print(f"\n  GPU candidates from '{cand_file}' ranked by |beta-0.5|:")
            print(f"  {'rank':>4}  {'p':>14}  {'omega':>8}  {'beta':>8}  {'|b-0.5|':>9}")
            print("  " + "-" * 54)
            for rank, (d, b, om, p) in enumerate(cands[:30], 1):
                edge_tag = " ***" if d < 0.02 else " **" if d < 0.08 else " *" if d < 0.15 else ""
                print(f"  {rank:4d}  {p:>14}  {om:8.4f}  {b:8.4f}  {d:9.4f}{edge_tag}")
            break
else:
    print(f"\n  (no candidate file found — run scorer to generate candidates)")

# =========================================================================
# 5. Summary
# =========================================================================
print()
print("=" * 78)
print("  RIEMANN EDGE SUMMARY")
print("=" * 78)
print(f"  Critical line: beta = 0.5 = Re(s) = 1/2 (Riemann Hypothesis)")
print(f"  Prism at p=136M: beta = {b_prism:.4f}, distance = {abs(b_prism-0.5):.4f}")
print(f"  p_critical (exact edge): {int(p_crit):,} (omega=7.5, beta=0.5)")
print(f"  BW=0.5 is NOT a tuning parameter -- it IS the Riemann half-strip width")
print(f"  The Gaussian window exp(-0.5*(d_omega/0.5)**2) is the spectral projection")
print(f"  onto the space spanned by primes near the critical line.")
print(f"  Moving the prism moves the center; the edge width (BW=0.5) is fixed by RH.")
