// ============================================================================
// HDGL φ-Resonance Scanner
// ============================================================================
//
// Scans Mersenne exponent resonance nodes, one best prime candidate per node.
//
// Scores:  S ≈ 0   → prime-like (destructive interference at Λ_φ = odd integer)
//          S ≈ 2   → composite-like
//          Gate:   S < 0.25 passes as a candidate
//
// Calibration: all 51 known Mersenne prime exponents are scored at startup.
// p = 22514117 is evaluated in isolation as an EXPERIMENTAL EXPONENT.
//
// Search: for each odd-integer lam node N (lam=41, 43, 45, ...)
//   - Computes the exact PASS window around the node centre
//   - Runs a segmented sieve to enumerate prime exponents in that window
//   - GPU-scans those prime exponents, returns single closest-to-centre PASS
//
// Compile:
//   nvcc -O3 -arch=sm_75 hdgl_scan.cu -o hdgl_scan
//
// Usage:
//   .\hdgl_scan.exe                    → from lam=41, advance until 10 candidates found
//   .\hdgl_scan.exe 43                 → from lam=43, advance until 10 candidates found
//   .\hdgl_scan.exe 41 51              → fixed range lam=41..51 (all nodes)
//   .\hdgl_scan.exe 41 0 5             → from lam=41, advance until 5 candidates found
//   .\hdgl_scan.exe 41 99 10           → lam=41..99, stop early after 10 candidates
// ============================================================================

#include <cuda_runtime.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

// ============================================================================
// BigG + Fudge10 empirical parameters
// k=1.049342  s0=0.994533  alpha=0.340052  Omega0=1.049675  beta=0.360942
//
// Scale-corrected φ-resonance (double precision)
// ============================================================================
//
// The raw X+1=0 gate S(p) = |e^(iπΛ_φ) + 1_eff| uses
//   δ(i) = |cos(π·β·φ)| · ln(n+2) / φ^{n+β}
// which decays exponentially in n = ⌊Λ_φ⌋.  For large p (p ~ 136M) n ~ 39,
// so φ^n ~ 10^8 and δ collapses to zero — the gate loses all sensitivity.
//
// BigG empirical finding: constants are scale-dependent.
//   Omega(p) = Omega0 / a(p)^alpha   where  a(p) = p_ref / p   (p_ref = 61)
//   s(p)     = s0 · (p / p_ref)^{-beta}
//
// This means at large p (small a), Omega *increases* and s *decreases* —
// the coupling strengthens, exactly compensating the δ suppression.
//
// Scale-corrected δ:
//   φ_eff(p) = φ · (Omega(p)/Omega0)^{1/k} · s(p)^{-1}
//            = φ · (p/p_ref)^{alpha/k} · (p/p_ref)^{beta}
//            = φ · (p/p_ref)^{(alpha+beta·k)/k}
//
//   δ_scaled = |cos(π·β·φ_eff)| · ln(n+2) / φ_eff^{n+β}
//
// This restores sensitivity at large p while remaining exact at p=61.

__device__ __host__ double phi_resonance_d(uint64_t p)
{
    const double PHI       = 1.6180339887498949;
    const double LN_PHI    = 0.4812118250596035;
    const double LN2_LNPHI = 1.4404200904125564;  // ln(2) / ln(φ) — Mersenne bridge
    const double INV_2PHI  = 0.3090169943749474;  // 1 / (2φ)
    const double PI        = 3.141592653589793238;

    // BigG + Fudge10 empirical parameters
    const double BIGG_K      = 1.049342;
    const double BIGG_S0     = 0.994533;   // (unused directly; folded into exponent below)
    const double BIGG_ALPHA  = 0.340052;
    const double BIGG_OMEGA0 = 1.049675;   // (ratio cancels; kept for clarity)
    const double BIGG_BETA   = 0.360942;   // s(p) evolution exponent
    const double P_REF       = 61.0;       // calibration anchor (M61 exact)

    (void)BIGG_S0; (void)BIGG_OMEGA0;      // ratio forms used, raw values cancel

    if (p < 2) return 2.0;

    double r    = (double)p;
    double M_U  = r * LN2_LNPHI;
    double lam  = log(M_U) / LN_PHI - INV_2PHI;  // Λ_φ

    double n_f  = floor(lam);
    double beta = lam - n_f;   // fractional part β

    // BigG scale factor: how many octaves above the p=61 calibration point.
    // Use log-scale ratio so it's well-behaved across the full range.
    double log_scale = log(r / P_REF) / LN_PHI;  // in φ-octaves

    // Scale exponent for φ_eff: (alpha + beta_bigg * k) / k
    // = alpha/k + beta_bigg
    double scale_exp = BIGG_ALPHA / BIGG_K + BIGG_BETA;  // ≈ 0.6854

    // Effective φ at this scale: φ_eff = φ · (p/p_ref)^scale_exp
    // To avoid overflow in φ_eff^{n+β}, work in log space throughout:
    //   log(φ_eff^{n+β}) = (n_f+1+β) * (LN_PHI + scale_exp * log(p/p_ref))
    double ln_phi_eff = LN_PHI + scale_exp * log(r / P_REF);

    // δ_scaled = |cos(π·β·φ_eff_arg)| · ln(n+2) / exp((n_val+β)·ln_phi_eff)
    // φ_eff_arg for cosine: we scale β-modulation by the same coupling strength
    double phi_eff_arg = PHI * exp(BIGG_ALPHA / BIGG_K * log(r / P_REF));
    double n_val       = n_f + 1.0;

    double denom = exp((n_val + beta) * ln_phi_eff);
    // Guard against underflow/overflow
    double delta;
    if (denom > 1e300) {
        delta = 0.0;
    } else if (denom < 1e-300) {
        delta = 1.0;  // saturate — strong coupling
    } else {
        delta = fabs(cos(PI * beta * phi_eff_arg))
              * log(n_val + 2.0)
              / denom;
    }

    // Clamp delta to [0, 1] — the correction can't exceed unity by construction
    if (delta > 1.0) delta = 1.0;

    double one_eff = 1.0 + delta;

    // S(p) = |e^{iπΛ_φ} + 1_eff|
    double re = cos(PI * lam) + one_eff;
    double im = sin(PI * lam);
    return sqrt(re * re + im * im);
}

// ============================================================================
// GPU kernel: batch resonance scan
// ============================================================================

__global__ void scan_kernel(
    const uint64_t * __restrict__ d_exps,
    double         * __restrict__ d_scores,
    int n)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i >= n) return;
    d_scores[i] = phi_resonance_d(d_exps[i]);
}

// ============================================================================
// Host helpers
// ============================================================================

typedef struct { uint64_t p; double score; } Result;

// ============================================================================
// Node geometry: compute PASS window for a given odd-integer lam node N
// ============================================================================

// Exact node centre: solve  lam = log(p * LN2_LNPHI)/LN_PHI - INV_2PHI = N
//   => p_centre = exp((N + INV_2PHI)*LN_PHI) / LN2_LNPHI
//
// PASS condition: S < 0.25.  For n_f = N-1 (even), cos(π*lam) = cos(π*β)
// and the gate passes when β > beta_thresh = (2/π)*acos(0.125) ≈ 0.9205.
// Window half-width in p-space: dβ/dp = 1/(p*LN_PHI)
//   => half = (1 - beta_thresh) * p_c * LN_PHI
static void node_geometry(int lam_n,
                          double  *p_centre_d,
                          uint64_t *win_lo, uint64_t *win_hi)
{
    const double LN_PHI    = 0.4812118250596035;
    const double LN2_LNPHI = 1.4404200904125564;
    const double INV_2PHI  = 0.3090169943749474;
    const double PI        = 3.141592653589793238;
    const double beta_thresh = 2.0 * acos(0.125) / PI;  // ≈0.9205

    double p_c = exp(((double)lam_n + INV_2PHI) * LN_PHI) / LN2_LNPHI;
    double half = (1.0 - beta_thresh) * p_c * LN_PHI;

    *p_centre_d = p_c;
    uint64_t lo = (uint64_t)(p_c - half);
    uint64_t hi = (uint64_t)(p_c + half);
    if (lo < 2) lo = 3;
    if ((lo & 1) == 0) lo++;   // ensure odd
    if ((hi & 1) == 0) hi++;   // ensure odd
    *win_lo = lo;
    *win_hi = hi;
}

// ============================================================================
// Segmented sieve: enumerate all primes in [lo, hi]
// Returns heap-allocated array (caller frees), sets *count.
// small_primes[] must cover primes up to sqrt(hi).
// ============================================================================

static uint64_t* sieve_range(uint64_t lo, uint64_t hi, int *count,
                              const uint32_t *sprimes, int n_sprimes)
{
    if (lo < 2) lo = 2;
    if (hi < lo) { *count = 0; return NULL; }

    const uint64_t SEG = 1ULL << 19;  // 512K bits = 64KB segment
    uint64_t range = hi - lo + 1;

    // Worst-case prime count estimate: range / ln(hi) * 1.5 for safety
    double est = (double)range / log((double)hi) * 1.5 + 100.0;
    if (est > 2e9) est = 2e9;
    uint64_t alloc = (uint64_t)est;
    uint64_t *out = (uint64_t*)malloc(alloc * sizeof(uint64_t));
    if (!out) { *count = 0; return NULL; }
    int n = 0;

    static char seg_buf[1 << 19];   // reuse across calls

    for (uint64_t seg_lo = lo; seg_lo <= hi; seg_lo += SEG) {
        uint64_t seg_hi = seg_lo + SEG - 1;
        if (seg_hi > hi) seg_hi = hi;
        uint64_t seg_len = seg_hi - seg_lo + 1;

        memset(seg_buf, 1, (size_t)seg_len);
        if (seg_lo == 0) seg_buf[0] = 0;  // 0 not prime
        if (seg_lo <= 1 && 1 <= seg_hi) seg_buf[1 - seg_lo] = 0;

        for (int si = 0; si < n_sprimes; si++) {
            uint64_t sp = sprimes[si];
            // First multiple of sp in [seg_lo, seg_hi]
            uint64_t start = ((seg_lo + sp - 1) / sp) * sp;
            if (start == sp) start += sp;  // skip sp itself
            for (uint64_t j = start; j <= seg_hi; j += sp)
                seg_buf[j - seg_lo] = 0;
        }
        for (uint64_t j = 0; j < seg_len; j++) {
            if (seg_buf[j]) {
                if ((uint64_t)n >= alloc) break;  // safety
                out[n++] = seg_lo + j;
            }
        }
    }
    *count = n;
    return out;
}

#define CUDA_CHECK(c) do {                                                    \
    cudaError_t _e = (c);                                                     \
    if (_e != cudaSuccess) {                                                  \
        fprintf(stderr, "[scan] CUDA error %s:%d: %s\n",                     \
                __FILE__, __LINE__, cudaGetErrorString(_e));                  \
        exit(EXIT_FAILURE);                                                   \
    }                                                                         \
} while (0)

int main(int argc, char **argv)
{
    // ---- Parse arguments: [lam_start=41] [lam_end=0] [n_goal=10] -----------
    // lam_end=0 (default): advance indefinitely until n_goal candidates found
    // lam_end>0: fixed range; still stops early when n_goal reached
    int lam_start = 41;
    int lam_end   = 0;   // 0 = open-ended
    int n_goal    = 10;  // stop after this many candidates
    if (argc >= 2) lam_start = atoi(argv[1]);
    if (argc >= 3) lam_end   = atoi(argv[2]);
    if (argc >= 4) n_goal    = atoi(argv[3]);
    if (lam_start < 3) lam_start = 3;
    if (n_goal < 1)    n_goal = 1;
    // Nodes must be odd integers (even lam gives cos(π*lam)=+1 → never cancels)
    if ((lam_start & 1) == 0) lam_start++;
    if (lam_end > 0 && (lam_end & 1) == 0) lam_end++;

    // ---- Build small prime table for segmented sieve -----------------------
    // Need primes up to sqrt(win_hi of largest requested node).
    // For lam up to 99: p_c ~ exp((99+0.309)*0.481)/1.440 ~ 1.3e14
    // sqrt(1.3e14) ~ 11.4M  -> sieve to 11.5M.
    // For safety sieve to 12M.
    const int SIEVE_LIMIT = 12000000;
    char *sv = (char*)calloc(SIEVE_LIMIT + 1, 1);
    if (!sv) { fprintf(stderr, "OOM: sieve\n"); return 1; }
    for (int i = 2; i <= SIEVE_LIMIT; i++) sv[i] = 1;
    for (int i = 2; (long long)i*i <= SIEVE_LIMIT; i++)
        if (sv[i]) for (int j = i*i; j <= SIEVE_LIMIT; j += i) sv[j] = 0;
    // Count primes
    int n_sp = 0;
    for (int i = 2; i <= SIEVE_LIMIT; i++) if (sv[i]) n_sp++;
    uint32_t *small_primes = (uint32_t*)malloc((size_t)n_sp * sizeof(uint32_t));
    if (!small_primes) { fprintf(stderr, "OOM: small_primes\n"); return 1; }
    int spi = 0;
    for (int i = 2; i <= SIEVE_LIMIT; i++) if (sv[i]) small_primes[spi++] = (uint32_t)i;
    free(sv);

    // ---- Calibration: score all known Mersenne prime exponents --------------

    static const uint64_t KNOWN_MP[] = {
        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
    };
    int n_known = (int)(sizeof(KNOWN_MP) / sizeof(KNOWN_MP[0]));

    printf("====================================================================\n");
    printf("  HDGL φ-Resonance Scanner — BigG+Fudge10 scale-corrected X+1=0 gate\n");
    printf("  S(p) → 0 : prime-like   S(p) → 2 : composite\n");
    printf("  Gate threshold S < 0.25  |  δ boosted by (p/61)^0.685\n");
    printf("====================================================================\n\n");

    printf("── Known Mersenne prime calibration (all %d) ────────────────────────\n", n_known);
    printf("  %-16s  %-14s  %-10s  %s\n", "p", "S(p)", "gate", "Λ_φ (n, β)");
    int n_pass_known = 0;
    for (int c = 0; c < n_known; c++) {
        uint64_t cp = KNOWN_MP[c];
        double sc = phi_resonance_d(cp);
        const double LN_PHI_C    = 0.4812118250596035;
        const double LN2_LNPHI_C = 1.4404200904125564;
        const double INV_2PHI_C  = 0.3090169943749474;
        double lam_c = log((double)cp * LN2_LNPHI_C) / LN_PHI_C - INV_2PHI_C;
        int pass = (sc < 0.25);
        if (pass) n_pass_known++;
        printf("  M%-15llu  %-14.6f  %-10s  n=%d \xCE\xB2=%.4f\n",
               (unsigned long long)cp, sc,
               pass ? "PASS \xE2\x9C\x93" : "FAIL",
               (int)floor(lam_c), lam_c - floor(lam_c));
    }
    printf("  ── Calibration: %d PASS / %d FAIL out of %d ──\n\n",
           n_pass_known, n_known - n_pass_known, n_known);

    // ---- Experimental exponent (isolated) -----------------------------------
    printf("── Experimental exponent (isolated, NOT a known prime) ────────────────\n");
    printf("  p = 22514117  [EXPERIMENTAL — not in calibration set]\n");
    double s_exp = phi_resonance_d(22514117ULL);
    printf("  S(22514117) = %.8f  →  %s\n\n", s_exp,
           (s_exp < 0.25) ? "PASS — STRONG candidate!" : "FAIL — gate rejects");

    // ---- Node search loop ---------------------------------------------------
    if (lam_end > 0)
        printf("── Node search: lam = %d to %d (fixed range) ──────────────────────────\n",
               lam_start, lam_end);
    else
        printf("── Node search: from lam = %d, advancing until %d candidate(s) found ──\n",
               lam_start, n_goal);
    printf("  One best prime candidate per node (closest to node centre)\n");
    printf("  Window capped at ~8M primes for nodes where full window is larger.\n");
    printf("  %-6s  %-16s  %-14s  %-12s  %-10s  %s\n",
           "node", "p_centre", "best_p", "S(best_p)", "|offset|", "status");
    printf("  %s\n", "------  ----------------  --------------  ------------  ----------  ------");

    int total_nodes = 0, nodes_with_candidate = 0, last_lam = lam_start;

    // Safety ceiling: lam grows by 2 per step; lam=999 puts p_c ~ 10^63
    for (int lam_n = lam_start; lam_n < lam_start + 2000; lam_n += 2) {
        // Stop conditions
        if (lam_end > 0 && lam_n > lam_end) break;
        if (nodes_with_candidate >= n_goal) break;
        last_lam = lam_n;
        total_nodes++;

        double  p_centre_d;
        uint64_t win_lo, win_hi;
        node_geometry(lam_n, &p_centre_d, &win_lo, &win_hi);
        uint64_t p_centre = (uint64_t)p_centre_d;

        // Cap window: never feed the GPU more than MAX_SCAN primes.
        // For large lam (p_c >> 1B) the full PASS window holds hundreds of
        // millions of primes; we shrink it symmetrically around p_centre so
        // GPU alloc stays ≤ ~256 MB and sieve time stays seconds not minutes.
        const int MAX_SCAN = 8000000;
        {
            double est = (double)(win_hi - win_lo) / log((double)p_centre_d);
            if (est > (double)MAX_SCAN) {
                double new_half = (double)MAX_SCAN * 0.5 * log((double)p_centre_d);
                uint64_t new_lo = (p_centre_d > new_half + 2.0)
                                  ? (uint64_t)(p_centre_d - new_half) : 2;
                uint64_t new_hi = (uint64_t)(p_centre_d + new_half);
                if ((new_lo & 1) == 0) new_lo++;
                if ((new_hi & 1) == 0) new_hi++;
                if (new_lo < win_lo) new_lo = win_lo;
                if (new_hi > win_hi) new_hi = win_hi;
                win_lo = new_lo; win_hi = new_hi;
            }
        }

        // Segmented sieve: enumerate prime exponents in this node's window
        int n_primes = 0;
        // We only need small_primes up to sqrt(win_hi) for this node
        uint64_t sq = (uint64_t)sqrt((double)win_hi) + 2;
        // Find cutoff in small_primes table
        int n_sp_node = 0;
        while (n_sp_node < n_sp && (uint64_t)small_primes[n_sp_node] <= sq)
            n_sp_node++;

        uint64_t *prime_exps = sieve_range(win_lo, win_hi, &n_primes,
                                           small_primes, n_sp_node);
        if (!prime_exps || n_primes == 0) {
            printf("  %-6d  %-16llu  (sieve empty or OOM)\n",
                   lam_n, (unsigned long long)p_centre);
            free(prime_exps);
            continue;
        }

        // GPU scan: score all prime candidates
        uint64_t *d_exps   = NULL;
        double   *d_scores = NULL;
        CUDA_CHECK(cudaMalloc(&d_exps,   (size_t)n_primes * sizeof(uint64_t)));
        CUDA_CHECK(cudaMalloc(&d_scores, (size_t)n_primes * sizeof(double)));
        CUDA_CHECK(cudaMemcpy(d_exps, prime_exps, (size_t)n_primes * sizeof(uint64_t),
                              cudaMemcpyHostToDevice));
        int blk = 256, grd = (n_primes + blk - 1) / blk;
        scan_kernel<<<grd, blk>>>(d_exps, d_scores, n_primes);
        CUDA_CHECK(cudaDeviceSynchronize());
        double *h_scores = (double*)malloc((size_t)n_primes * sizeof(double));
        if (!h_scores) { fprintf(stderr, "OOM: h_scores\n"); return 1; }
        CUDA_CHECK(cudaMemcpy(h_scores, d_scores, (size_t)n_primes * sizeof(double),
                              cudaMemcpyDeviceToHost));
        cudaFree(d_exps); cudaFree(d_scores);

        // Find single best: PASS (S<0.25) closest to node centre
        uint64_t best_p    = 0;
        double   best_s    = 2.0;
        uint64_t best_dist = UINT64_MAX;
        for (int i = 0; i < n_primes; i++) {
            if (h_scores[i] >= 0.25) continue;
            uint64_t cp   = prime_exps[i];
            uint64_t dist = (cp >= p_centre) ? (cp - p_centre) : (p_centre - cp);
            if (dist < best_dist) {
                best_dist = dist; best_p = cp; best_s = h_scores[i];
            }
        }
        free(h_scores); free(prime_exps);

        if (best_p == 0) {
            printf("  %-6d  %-16llu  (no prime PASS in window)\n",
                   lam_n, (unsigned long long)p_centre);
        } else {
            nodes_with_candidate++;
            printf("  %-6d  %-16llu  %-14llu  %-12.8f  %-10llu  PASS \xE2\x9C\x93 submit mersenne.org\n",
                   lam_n,
                   (unsigned long long)p_centre,
                   (unsigned long long)best_p,
                   best_s,
                   (unsigned long long)best_dist);
        }
    }

    printf("\n── Summary ──────────────────────────────────────────────────────────\n");
    printf("  Nodes scanned: lam = %d .. %d  (%d nodes)\n", lam_start, last_lam, total_nodes);
    printf("  Nodes with prime candidate: %d / %d\n", nodes_with_candidate, total_nodes);
    if (nodes_with_candidate >= n_goal)
        printf("  Goal of %d candidate(s) reached.\n", n_goal);
    printf("====================================================================\n");

    free(small_primes);
    return 0;
}
