/*
 * mersenne_sieve.c
 * ============================================================
 * MERSENNE CANDIDATE SIEVE — correct and fast
 *
 * For each candidate exponent p, this tool:
 *
 *  Stage 1 — Exponent primality (instant)
 *    p must itself be prime for 2^p-1 to potentially be prime.
 *    Uses deterministic Miller-Rabin (64-bit, 7 witnesses, proven correct).
 *
 *  Stage 2 — Trial factor sieve (seconds to minutes)
 *    Any factor q of 2^p-1 must satisfy:
 *      q = 2kp + 1  for some k ≥ 1
 *      q ≡ 1 or 7 (mod 8)
 *    We test all such q up to a configurable limit using fast
 *    64-bit modular exponentiation: if 2^p ≡ 1 (mod q), then q | 2^p-1.
 *    This eliminates most composites instantly.
 *
 *  Stage 3 — Fermat probable prime (minutes on CPU, fast on GPU)
 *    Compute 3^(M_p - 1) mod M_p where M_p = 2^p - 1.
 *    Uses Barrett-style mod 2^p-1 reduction (shift-and-add).
 *    False positive rate: ≤ 4^{-1} for a single base.
 *    Run multiple bases for stronger confidence.
 *    NOTE: This is NOT Lucas-Lehmer. It is a probabilistic test.
 *    For a definitive answer, use gpuowl or mlucas.
 *
 *  Stage 4 — PLN resonance score (instant)
 *    Scores p by PLN proximity to an anchor prime.
 *    This is a HEURISTIC RANKING, not a primality test.
 *
 *  Output: ranked candidate list with sieve results
 *    Eliminated candidates are marked COMPOSITE.
 *    Survivors are ranked by PLN score for GPU submission.
 *
 * Compile:
 *   gcc -O3 -std=c99 -o mersenne_sieve mersenne_sieve.c -lm
 *   (no GMP needed — all arithmetic is 64-bit)
 *
 * Usage:
 *   mersenne_sieve --sieve 136279841 8.0 200000 70        # sweep + sieve to 2^70
 *   mersenne_sieve --sieve-file mersenne_candidates.txt 64 # sieve pln_final output
 *   mersenne_sieve --test-p 272559689 272559697 272559701  # test specific exponents
 *   mersenne_sieve --factor-limit 272559689 68             # trial factor to 2^68
 *   mersenne_sieve --help
 *
 * Pipeline with pln_final:
 *   pln --prime-file-candidates M136279841.txt 20 0.005 12
 *   mersenne_sieve --sieve-file mersenne_candidates.txt 64
 *
 * Honest certainty levels:
 *   COMPOSITE (factor found):    100% certain composite
 *   COMPOSITE (p not prime):     100% certain composite
 *   SURVIVED sieve to 2^B:       No factor below 2^B — not a primality statement
 *   Fermat PASS (k bases):       P(false positive) ≤ 4^{-k} — probabilistic only
 *   Lucas-Lehmer:                100% deterministic — requires GPU (months)
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include <time.h>

/* ================================================================
 * CONSTANTS
 * ================================================================ */
#define PHI      1.6180339887498948482L
#define LN_PHI   0.4812118250596034748L
#define LN2      0.6931471805599453094L
#define LOG10_2  0.30102999566398119521L

/* ================================================================
 * FAST 64-BIT MODULAR ARITHMETIC
 * ================================================================ */
static inline uint64_t mulmod64(uint64_t a, uint64_t b, uint64_t m) {
    return (uint64_t)((__uint128_t)a * b % m);
}

static uint64_t powmod64(uint64_t base, uint64_t exp, uint64_t mod) {
    if (mod == 1) return 0;
    uint64_t result = 1;
    base %= mod;
    while (exp) {
        if (exp & 1) result = mulmod64(result, base, mod);
        base = mulmod64(base, base, mod);
        exp >>= 1;
    }
    return result;
}

/* ================================================================
 * STAGE 1: EXPONENT PRIMALITY TEST
 * Deterministic Miller-Rabin for 64-bit numbers.
 * Proven correct for all n < 3,317,044,064,679,887,385,961,981.
 * ================================================================ */
static bool is_prime_u64(uint64_t n) {
    if (n < 2) return false;
    if (n == 2 || n == 3 || n == 5 || n == 7) return true;
    if (n % 2 == 0 || n % 3 == 0 || n % 5 == 0) return false;

    /* Trial division to 1000 first */
    for (uint64_t d = 7; d <= 1000 && d * d <= n; d += 2) {
        if (n % d == 0) return false;
    }

    uint64_t d = n - 1;
    int s = 0;
    while (!(d & 1)) { d >>= 1; s++; }

    static const uint64_t bases[] = {
        2, 325, 9375, 28178, 450775, 9780504, 1795265022
    };
    for (int i = 0; i < 7; i++) {
        uint64_t a = bases[i] % n;
        if (a == 0) return true;
        uint64_t x = powmod64(a, d, n);
        if (x == 1 || x == n - 1) continue;
        bool composite = true;
        for (int r = 1; r < s; r++) {
            x = mulmod64(x, x, n);
            if (x == n - 1) { composite = false; break; }
        }
        if (composite) return false;
    }
    return true;
}

/* ================================================================
 * STAGE 2: TRIAL FACTOR SIEVE
 *
 * Tests all q = 2kp+1 with q ≡ 1 or 7 (mod 8) and q prime-like,
 * up to trial_limit.
 *
 * Key insight: 2^p ≡ 1 (mod q) iff q | 2^p-1.
 * pow(2, p, q) is O(log p) multiplications mod q — microseconds.
 *
 * Returns: factor if found (> 1), or 0 if none found below limit.
 * ================================================================ */
static uint64_t trial_factor(uint64_t p, uint64_t trial_limit) {
    /* k starts at 1: q = 2p+1, 4p+1, 6p+1, ... */
    for (uint64_t k = 1; ; k++) {
        /* q = 2kp + 1 */
        /* Overflow check: if 2kp+1 > trial_limit or overflows u64 */
        if (k > (trial_limit - 1) / (2 * p)) break;
        uint64_t q = 2 * k * p + 1;
        if (q > trial_limit) break;

        /* Must satisfy q ≡ 1 or 7 (mod 8) */
        uint64_t qmod8 = q % 8;
        if (qmod8 != 1 && qmod8 != 7) continue;

        /* Quick composite rejection: divisible by 3, 5, 7? */
        if (q % 3 == 0 || q % 5 == 0 || q % 7 == 0 ||
            q % 11 == 0 || q % 13 == 0) {
            if (q > 13) continue;
        }

        /* Check 2^p ≡ 1 (mod q) */
        if (powmod64(2, p, q) == 1) {
            return q; /* q divides 2^p - 1 */
        }
    }
    return 0; /* no factor found */
}

/* ================================================================
 * STAGE 3: FERMAT PROBABLE PRIME TEST FOR MERSENNE NUMBERS
 *
 * Tests: a^(M_p - 1) ≡ 1 (mod M_p) where M_p = 2^p - 1.
 * Equivalently: a^(2^p - 2) ≡ 1 (mod 2^p - 1)
 * Which means: a^(2^(p-1) - 1) ≡ 1 (mod 2^p - 1) [since M_p - 1 = 2(2^(p-1)-1)]
 *
 * For Mersenne numbers, x mod (2^p-1) = (x >> p) + (x & (2^p-1))
 * exploiting the structure for fast reduction.
 *
 * PRACTICAL LIMIT: For p > ~25M, this requires GMP or GPU arithmetic.
 * For p ≈ 272M on CPU: each squaring is ~272M bits → not feasible on CPU.
 *
 * This implementation uses GMP if available (compile with -DUSE_GMP -lgmp),
 * otherwise falls back gracefully with a clear warning.
 *
 * Returns:
 *   1  = passed (probable prime for this base)
 *   0  = failed (definitely composite)
 *  -1  = skipped (too large for CPU without GMP)
 * ================================================================ */
#ifdef USE_GMP
#include <gmp.h>

static int fermat_mersenne_gmp(uint64_t p, uint64_t base) {
    mpz_t M, a, result;
    mpz_init(M); mpz_init(a); mpz_init(result);

    /* M = 2^p - 1 */
    mpz_set_ui(M, 1);
    mpz_mul_2exp(M, M, (mp_bitcnt_t)p);
    mpz_sub_ui(M, M, 1);

    /* exponent = M - 1 = 2^p - 2 */
    mpz_t exp;
    mpz_init_set(exp, M);
    mpz_sub_ui(exp, exp, 1);

    mpz_set_ui(a, base);
    mpz_powm(result, a, exp, M);

    int pass = (mpz_cmp_ui(result, 1) == 0) ? 1 : 0;

    mpz_clear(M); mpz_clear(a); mpz_clear(result); mpz_clear(exp);
    return pass;
}
#endif

static int fermat_mersenne(uint64_t p, uint64_t base) {
    /* Practical feasibility: O(p^2) bit operations */
    /* For p > 10M on CPU without GPU: skip */
    if (p > 10000000ULL) {
#ifdef USE_GMP
        return fermat_mersenne_gmp(p, base);
#else
        return -1; /* too large for CPU-only */
#endif
    }
    /* For small p: compute directly */
    /* This path handles only small verification cases */
    return -1; /* defer to GPU for all Mersenne sizes we care about */
}

/* ================================================================
 * PLN COORDINATE (log-space, instant)
 * ================================================================ */
typedef struct {
    uint64_t    p;
    long double omega;
    int         D;
    long double theta_pi;
    long double r;
    long double score;   /* resonance vs anchor */
    uint64_t    digits;
} PLNInfo;

static long double orbit_r(long double omega) {
    long double frac = omega - floorl(omega);
    long double theta = frac * 6.2831853071795864769L;
    long double D_cont = fmodl(omega, 8.0L) + 1.0L;
    long double lat = PHI*PHI*(1.0L - powl(PHI, -D_cont)) + 1.0L/D_cont;
    lat = fmaxl(1.0L, fminl(lat, PHI*PHI));
    long double Rre = cosl(theta) - 1.0L + lat;
    long double Rim = sinl(theta);
    return sqrtl(Rre*Rre + Rim*Rim);
}

static PLNInfo compute_pln(uint64_t p) {
    PLNInfo m;
    m.p = p;
    long double log2x = (long double)p;
    long double lnx = log2x * LN2;
    long double arg = lnx / LN_PHI;
    m.omega     = (arg > 1.0L) ? logl(arg) / LN_PHI : 0.0L;
    m.D         = (int)fmodl(floorl(m.omega), 8.0L) + 1;
    m.theta_pi  = fmodl(m.omega, 1.0L);
    m.r         = orbit_r(m.omega);
    m.digits    = (uint64_t)floorl((long double)p * LOG10_2) + 1;
    m.score     = 0.0L;
    return m;
}

static long double pln_resonance(const PLNInfo *a, const PLNInfo *anchor) {
    /* fractional-omega phase distance (was split 0.5+0.2 on identical quantity) */
    long double dw = fabsl(fmodl(a->omega,1.0L) - fmodl(anchor->omega,1.0L));
    if (dw > 0.5L) dw = 1.0L - dw;
    /* orbit radius distance relative to anchor */
    long double dr = fabsl(a->r - anchor->r);
    /* absolute golden orbit proximity: Gaussian centered at |R|=1/phi (~0.618)
     * This term is INDEPENDENT of the anchor and adds real per-candidate signal.
     * sigma ~ 0.2: exp(-dr^2 / (2*0.2^2)) = exp(-dr^2 * 24) */
    long double dr_gold = a->r - 0.6180339887L;
    long double gold = expl(-dr_gold * dr_gold * 24.0L);
    long double score = 0.50L/(1.0L+dw) + 0.20L/(1.0L+dr) + 0.25L*gold;
    if (a->D == anchor->D) score += 0.05L;
    return score;
}

/* ================================================================
 * RESULT RECORD
 * ================================================================ */
typedef enum {
    STATUS_UNKNOWN,
    STATUS_P_COMPOSITE,   /* p itself is composite — 2^p-1 definitely composite */
    STATUS_FACTOR_FOUND,  /* small factor of 2^p-1 found */
    STATUS_SURVIVED,      /* survived trial factoring to limit */
} SieveStatus;

typedef struct {
    PLNInfo     pln;
    SieveStatus status;
    uint64_t    factor;       /* set if STATUS_FACTOR_FOUND */
    uint64_t    sieve_limit;  /* highest factor tested */
    double      sieve_bits;   /* log2(sieve_limit) */
    double      elapsed_s;
} CandResult;

static int cmp_results(const void *a, const void *b) {
    const CandResult *ra = (const CandResult *)a;
    const CandResult *rb = (const CandResult *)b;
    int sa = (ra->status == STATUS_SURVIVED) ? 0 : 1;
    int sb = (rb->status == STATUS_SURVIVED) ? 0 : 1;
    if (sa != sb) return sa - sb;
    if (rb->pln.score > ra->pln.score) return 1;
    if (rb->pln.score < ra->pln.score) return -1;
    return 0;
}

static void print_result(const CandResult *r, int rank) {
    const char *status_str;
    switch (r->status) {
        case STATUS_P_COMPOSITE:  status_str = "COMPOSITE(p)"; break;
        case STATUS_FACTOR_FOUND: status_str = "COMPOSITE(factor)"; break;
        case STATUS_SURVIVED:     status_str = "SURVIVED"; break;
        default:                  status_str = "UNKNOWN"; break;
    }

    printf("  %3d  p=%-14llu  D%d  ω=%.8Lf  digits=%llu\n",
           rank, (unsigned long long)r->pln.p,
           r->pln.D, r->pln.omega, (unsigned long long)r->pln.digits);
    printf("       score=%.5Lf  |R|=%.5Lf  status=%-20s",
           r->pln.score, r->pln.r, status_str);
    if (r->status == STATUS_FACTOR_FOUND)
        printf("  factor=%llu", (unsigned long long)r->factor);
    else if (r->status == STATUS_SURVIVED)
        printf("  sieved_to=2^%.1f", r->sieve_bits);
    printf("  (%.2fs)\n", r->elapsed_s);
}

/* ================================================================
 * MAIN SIEVE FUNCTION
 * ================================================================ */
static void run_sieve(uint64_t anchor_p, long double omega_radius,
                      int sweep_steps, double factor_bits) {
    PLNInfo anchor = compute_pln(anchor_p);
    uint64_t trial_limit;

    /* Build trial limit from bits */
    if (factor_bits >= 63.0) {
        trial_limit = UINT64_MAX;
    } else {
        trial_limit = (uint64_t)1 << (int)factor_bits;
    }

    printf("=== MERSENNE CANDIDATE SIEVE ===\n");
    printf("Anchor: p=%llu  omega=%.12Lf  D=%d\n",
           (unsigned long long)anchor_p, anchor.omega, anchor.D);
    printf("Sweep:  omega ± %.4Lf  steps=%d\n", omega_radius, sweep_steps);
    printf("Trial factor limit: 2^%.1f = %llu\n\n", factor_bits, (unsigned long long)trial_limit);

    /* Generate candidate p values — dynamic allocation, O(1) dedup
     * (p is non-decreasing with omega, so only the last value can duplicate) */
    uint64_t *cands = (uint64_t *)malloc((size_t)sweep_steps * sizeof(uint64_t));
    if (!cands) { fprintf(stderr, "Out of memory\n"); return; }
    int ncands = 0;
    for (int i = 0; i < sweep_steps; i++) {
        long double t = (long double)i / (sweep_steps - 1);
        long double omega = anchor.omega - omega_radius + 2.0L * omega_radius * t;
        if (omega <= 0.0L) continue;
        long double log2x = expl(omega * LN_PHI) * LN_PHI / LN2;
        uint64_t p = (uint64_t)llroundl(log2x);
        if (p <= anchor_p) continue;
        if (ncands > 0 && cands[ncands-1] == p) continue; /* monotonic dedup */
        cands[ncands++] = p;
    }

    printf("Generated %d candidate exponents above anchor.\n", ncands);
    printf("Running sieve...\n\n");

    /* Process each candidate */
    CandResult *results = calloc(ncands, sizeof(CandResult));
    int n_composite_p = 0, n_factor_found = 0, n_survived = 0;

    for (int i = 0; i < ncands; i++) {
        uint64_t p = cands[i];
        clock_t t0, t1;
        t0 = clock();

        CandResult *r = &results[i];
        r->pln = compute_pln(p);

        /* Compute PLN resonance vs anchor */
        r->pln.score = pln_resonance(&r->pln, &anchor);

        /* Stage 1: is p prime? */
        if (!is_prime_u64(p)) {
            r->status = STATUS_P_COMPOSITE;
            r->factor = 0;
            n_composite_p++;
        } else {
            /* Stage 2: trial factoring */
            uint64_t f = trial_factor(p, trial_limit);
            if (f > 1) {
                r->status = STATUS_FACTOR_FOUND;
                r->factor = f;
                n_factor_found++;
            } else {
                r->status = STATUS_SURVIVED;
                r->sieve_limit = trial_limit;
                r->sieve_bits = factor_bits;
                n_survived++;
            }
        }

        t1 = clock();
        r->elapsed_s = (double)(t1-t0)/CLOCKS_PER_SEC;

        /* Progress */
        if ((i+1) % 50 == 0 || i == ncands-1) {
            printf("\r  Processed %d/%d ...", i+1, ncands);
            fflush(stdout);
        }
    }
    printf("\n\n");

    /* Sort: survived first, then by PLN score descending */
    qsort(results, (size_t)ncands, sizeof(CandResult), cmp_results);

    /* Print summary */
    printf("=== SIEVE RESULTS ===\n");
    printf("Total candidates: %d\n", ncands);
    printf("Eliminated (p composite):     %d\n", n_composite_p);
    printf("Eliminated (factor found):    %d\n", n_factor_found);
    printf("Survived (no factor to 2^%.1f): %d\n\n", factor_bits, n_survived);

    /* Print honest certainty statement */
    printf("=== CERTAINTY STATEMENT ===\n");
    printf("COMPOSITE results are 100%% certain.\n");
    printf("SURVIVED results: no factor q = 2kp+1 (q ≡ 1,7 mod 8) below 2^%.1f.\n", factor_bits);
    printf("This does NOT imply primality. Surviving candidates need:\n");
    printf("  → GPU trial factoring to 2^76 (mfactor/GIMPS)\n");
    printf("  → Full Lucas-Lehmer test (gpuowl/mlucas, months on GPU)\n");
    printf("PLN scores are HEURISTIC RANKINGS only, not primality evidence.\n\n");

    /* Print top survivors */
    int shown = 0;
    printf("=== TOP SURVIVORS (ranked by PLN score) ===\n");
    printf("  These are the best candidates to submit to gpuowl:\n\n");
    for (int i = 0; i < ncands && shown < 30; i++) {
        if (results[i].status != STATUS_SURVIVED) continue;
        print_result(&results[i], shown+1);
        shown++;
    }

    /* Write gpuowl commands */
    FILE *f = fopen("mersenne_sieve_candidates.txt", "w");
    if (f) {
        fprintf(f, "# Mersenne Sieve Candidates\n");
        fprintf(f, "# Anchor: p=%llu  sieved_to=2^%.1f\n",
                (unsigned long long)anchor_p, factor_bits);
        fprintf(f, "# CERTAINTY: composite results are final;\n");
        fprintf(f, "# SURVIVED results need GPU LL test for primality.\n\n");
        int written = 0;
        for (int i = 0; i < ncands && written < 100; i++) {
            if (results[i].status != STATUS_SURVIVED) continue;
            fprintf(f, "p = %-14llu  # omega=%.10Lf  D%d  score=%.5Lf  digits=%llu\n",
                    (unsigned long long)results[i].pln.p,
                    results[i].pln.omega, results[i].pln.D,
                    results[i].pln.score,
                    (unsigned long long)results[i].pln.digits);
            fprintf(f, "   gpuowl -p %llu\n\n",
                    (unsigned long long)results[i].pln.p);
            written++;
        }
        fclose(f);
        printf("\nSaved: mersenne_sieve_candidates.txt\n");
    }

    /* Also print eliminated composites briefly */
    if (n_factor_found > 0) {
        printf("\n=== CONFIRMED COMPOSITES (factor found) ===\n");
        int shown2 = 0;
        for (int i = 0; i < ncands && shown2 < 10; i++) {
            if (results[i].status != STATUS_FACTOR_FOUND) continue;
            printf("  p=%-14llu  COMPOSITE: %llu × (2^p-1)/%-14llu\n",
                   (unsigned long long)results[i].pln.p,
                   (unsigned long long)results[i].factor,
                   (unsigned long long)results[i].factor);
            shown2++;
        }
    }

    free(results);
    free(cands);
}

/* ================================================================
 * TEST SPECIFIC EXPONENTS
 * ================================================================ */
static void test_specific(uint64_t *ps, int count, double factor_bits) {
    uint64_t trial_limit = (factor_bits >= 63.0) ?
        UINT64_MAX : ((uint64_t)1 << (int)factor_bits);

    printf("=== SPECIFIC EXPONENT TESTS ===\n");
    printf("Trial factor limit: 2^%.1f\n\n", factor_bits);
    printf("  Rank  p               p_prime  sieve_result      factor\n");
    printf("  ----  --------------  -------  ----------------  ------\n");

    for (int i = 0; i < count; i++) {
        uint64_t p = ps[i];
        bool p_prime = is_prime_u64(p);
        PLNInfo pln = compute_pln(p);

        clock_t t0, t1;
        t0 = clock();

        if (!p_prime) {
            t1 = clock();
            double el = (double)(t1-t0)/CLOCKS_PER_SEC;
            printf("  %4d  %-14llu  NO       COMPOSITE(p is not prime)  (%.3fs)\n",
                   i+1, (unsigned long long)p, el);
            printf("        omega=%.10Lf  D%d  digits=%llu\n",
                   pln.omega, pln.D, (unsigned long long)pln.digits);
            continue;
        }

        uint64_t f = trial_factor(p, trial_limit);
        t1 = clock();
        double el = (double)(t1-t0)/CLOCKS_PER_SEC;

        if (f > 1) {
            printf("  %4d  %-14llu  YES      COMPOSITE(factor=%llu)  (%.3fs)\n",
                   i+1, (unsigned long long)p, (unsigned long long)f, el);
        } else {
            printf("  %4d  %-14llu  YES      SURVIVED to 2^%.1f          (%.3fs)\n",
                   i+1, (unsigned long long)p, factor_bits, el);
        }
        printf("        omega=%.10Lf  D%d  digits=%llu\n",
               pln.omega, pln.D, (unsigned long long)pln.digits);
    }

    printf("\nCERTAINTY: COMPOSITE results are final.\n");
    printf("SURVIVED: needs GPU LL test (gpuowl/mlucas). No CPU shortcut exists.\n");
}

/* ================================================================
 * SIEVE FROM PLN CANDIDATES FILE
 *
 * Reads p values from mersenne_candidates.txt (the output of
 * pln_final --prime-file-candidates) and runs the trial factor
 * sieve on each, then re-ranks survivors by PLN score.
 *
 * File format expected (same as pln_final output):
 *   # Anchor: ...  known_p=136279841  omega=...
 *   p = 136285088      # omega=...  score=...  ...
 *   gpuowl -p 136285088
 *
 * Pipeline:
 *   pln --prime-file-candidates M136279841.txt 20 0.005 12
 *   mersenne_sieve --sieve-file mersenne_candidates.txt 64
 * ================================================================ */
static void sieve_from_file(const char *path, double factor_bits) {
    FILE *f = fopen(path, "r");
    if (!f) {
        fprintf(stderr, "Error: cannot open '%s'\n", path);
        return;
    }

    uint64_t anchor_p = 0;
    uint64_t cands[200];
    int ncands = 0;

    char line[512];
    while (fgets(line, sizeof(line), f) && ncands < 200) {
        /* Parse anchor from pln_final comment: "# Anchor: ...  known_p=136279841  ..." */
        if (line[0] == '#') {
            char *kp = strstr(line, "known_p=");
            if (kp) anchor_p = strtoull(kp + 8, NULL, 10);
            continue;
        }
        /* Parse candidate line: "p = 136285088   # ..." */
        if (strncmp(line, "p = ", 4) == 0 || strncmp(line, "p =", 3) == 0) {
            uint64_t p = strtoull(line + 3, NULL, 10);
            if (p > 1) {
                bool dup = false;
                for (int j = 0; j < ncands; j++) if (cands[j] == p) { dup = true; break; }
                if (!dup) cands[ncands++] = p;
            }
        }
    }
    fclose(f);

    if (ncands == 0) {
        printf("No candidates found in '%s'\n", path);
        printf("Expected lines of the form: p = <number>  # ...\n");
        return;
    }

    if (anchor_p == 0) {
        /* Fall back to the smallest candidate as anchor */
        anchor_p = cands[0];
        for (int i = 1; i < ncands; i++)
            if (cands[i] < anchor_p) anchor_p = cands[i];
    }

    PLNInfo anchor = compute_pln(anchor_p);
    uint64_t trial_limit = (factor_bits >= 63.0) ?
        UINT64_MAX : ((uint64_t)1 << (int)factor_bits);

    printf("=== MERSENNE SIEVE (from file: %s) ===\n", path);
    printf("Anchor: p=%llu  omega=%.12Lf  D=%d\n",
           (unsigned long long)anchor_p, anchor.omega, anchor.D);
    printf("Candidates: %d   Trial factor limit: 2^%.1f\n\n", ncands, factor_bits);

    CandResult *results = calloc(ncands, sizeof(CandResult));
    int n_composite_p = 0, n_factor_found = 0, n_survived = 0;

    for (int i = 0; i < ncands; i++) {
        uint64_t p = cands[i];
        clock_t t0 = clock();

        CandResult *r = &results[i];
        r->pln = compute_pln(p);
        r->pln.score = pln_resonance(&r->pln, &anchor);

        if (!is_prime_u64(p)) {
            r->status = STATUS_P_COMPOSITE;
            n_composite_p++;
        } else {
            uint64_t fac = trial_factor(p, trial_limit);
            if (fac > 1) {
                r->status = STATUS_FACTOR_FOUND;
                r->factor = fac;
                n_factor_found++;
            } else {
                r->status = STATUS_SURVIVED;
                r->sieve_limit = trial_limit;
                r->sieve_bits = factor_bits;
                n_survived++;
            }
        }
        r->elapsed_s = (double)(clock() - t0) / CLOCKS_PER_SEC;

        printf("\r  Processed %d/%d ...", i + 1, ncands);
        fflush(stdout);
    }
    printf("\n\n");

    /* Sort survivors by PLN score descending */
    for (int i = 1; i < ncands; i++) {
        CandResult tmp = results[i];
        int j = i - 1;
        while (j >= 0 && results[j].pln.score < tmp.pln.score &&
               results[j].status == STATUS_SURVIVED && tmp.status == STATUS_SURVIVED) {
            results[j+1] = results[j]; j--;
        }
        results[j+1] = tmp;
    }

    printf("=== SIEVE RESULTS ===\n");
    printf("Total: %d   Eliminated(p composite): %d   Eliminated(factor): %d   Survived: %d\n\n",
           ncands, n_composite_p, n_factor_found, n_survived);

    printf("=== CERTAINTY STATEMENT ===\n");
    printf("COMPOSITE results are 100%% certain.\n");
    printf("SURVIVED: no factor q = 2kp+1 (q ≡ 1,7 mod 8) below 2^%.1f.\n", factor_bits);
    printf("This does NOT imply primality — needs GPU LL test (gpuowl/mlucas).\n");
    printf("PLN scores are HEURISTIC RANKINGS only, not primality evidence.\n\n");

    printf("=== SURVIVORS (ranked by PLN score) ===\n");
    int shown = 0;
    for (int i = 0; i < ncands && shown < 20; i++) {
        if (results[i].status != STATUS_SURVIVED) continue;
        print_result(&results[i], shown + 1);
        shown++;
    }

    if (n_factor_found > 0) {
        printf("\n=== CONFIRMED COMPOSITES (factor found) ===\n");
        int shown2 = 0;
        for (int i = 0; i < ncands && shown2 < 5; i++) {
            if (results[i].status != STATUS_FACTOR_FOUND) continue;
            printf("  p=%-14llu  COMPOSITE: factor=%llu\n",
                   (unsigned long long)results[i].pln.p,
                   (unsigned long long)results[i].factor);
            shown2++;
        }
    }

    /* Write sieve-filtered candidate list */
    FILE *out = fopen("mersenne_sieve_candidates.txt", "w");
    if (out) {
        fprintf(out, "# Mersenne Sieve Results (input: %s)\n", path);
        fprintf(out, "# Anchor: p=%llu  sieved_to=2^%.1f\n",
                (unsigned long long)anchor_p, factor_bits);
        fprintf(out, "# CERTAINTY: composite results final; SURVIVED needs GPU LL test.\n\n");
        int written = 0;
        for (int i = 0; i < ncands; i++) {
            if (results[i].status != STATUS_SURVIVED) continue;
            fprintf(out, "p = %-14llu  # omega=%.10Lf  D%d  score=%.5Lf  sieved_to=2^%.1f  digits=%llu\n",
                    (unsigned long long)results[i].pln.p,
                    results[i].pln.omega, results[i].pln.D,
                    results[i].pln.score, factor_bits,
                    (unsigned long long)results[i].pln.digits);
            fprintf(out, "   gpuowl -p %llu\n\n",
                    (unsigned long long)results[i].pln.p);
            written++;
        }
        fclose(out);
        printf("\nSaved: mersenne_sieve_candidates.txt  (%d survivors)\n", written);
    }

    free(results);
}

/* ================================================================
 * CROSS-VALIDATION  (all 51 known Mersenne prime exponents)
 * See function below for full documentation.
 * ================================================================ */
#define CV_SAMPLE_N  50000   /* prime exponents sampled per anchor */
#define CV_WIN       0.30    /* ±30% window to capture more pairs */
#define CV_MAX_K     64      /* max known primes */

/* All 51 known Mersenne prime exponents (OEIS A000043) */
static const uint64_t ALL_MERSENNE_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
};
#define N_ALL_MERSENNE (int)(sizeof(ALL_MERSENNE_P)/sizeof(ALL_MERSENNE_P[0]))

/* Per-prime cross-validation accumulator */
typedef struct {
    uint64_t p;
    double   pctiles[CV_MAX_K]; /* percentile each time this p was a target */
    int      n;                  /* number of times tested as target */
} PrimeCV;

static int cmp_double_desc(const void *a, const void *b) {
    double da = *(const double*)a, db = *(const double*)b;
    return (db > da) - (db < da);
}

static double median_d(double *arr, int n) {
    if (n <= 0) return -1.0;
    double tmp[CV_MAX_K]; memcpy(tmp, arr, n * sizeof(double));
    qsort(tmp, n, sizeof(double), cmp_double_desc); /* desc */
    /* invert for ascending sort of raw percentiles */
    for (int i = 0; i < n/2; i++) { double t=tmp[i]; tmp[i]=tmp[n-1-i]; tmp[n-1-i]=t; }
    return n % 2 == 1 ? tmp[n/2] : (tmp[n/2-1]+tmp[n/2])/2.0;
}

static void cross_validate(uint64_t *extra, int nextra) {
    /* Build working set: all-51 plus any user-supplied extras (dedup) */
    uint64_t kp[CV_MAX_K * 2]; int nk = 0;
    for (int i = 0; i < N_ALL_MERSENNE && nk < CV_MAX_K*2-1; i++) kp[nk++] = ALL_MERSENNE_P[i];
    for (int i = 0; i < nextra && nk < CV_MAX_K*2-1; i++) {
        bool dup = false;
        for (int j = 0; j < nk; j++) if (kp[j] == extra[i]) { dup = true; break; }
        if (!dup) kp[nk++] = extra[i];
    }
    /* Sort ascending */
    for (int i = 1; i < nk; i++) {
        uint64_t t = kp[i]; int j = i-1;
        while (j >= 0 && kp[j] > t) { kp[j+1] = kp[j]; j--; }
        kp[j+1] = t;
    }

    printf("MERSENNE CANDIDATE SIEVE — honest certainty edition\n\n");
    printf("=== PLN CROSS-VALIDATION ===\n");
    printf("Known primes used : %d (all 51 OEIS A000043", nk);
    if (nextra > 0) printf(" + %d user", nextra);
    printf(")\n");
    printf("Sample per anchor : %d prime exponents\n", CV_SAMPLE_N);
    printf("Window            : ±%.0f%% of anchor p\n", CV_WIN*100);
    printf("Small primes      : computed on-the-fly (no digit files needed)\n\n");

    /* Per-prime accumulators */
    PrimeCV cv[CV_MAX_K*2]; memset(cv, 0, sizeof(cv));
    for (int i = 0; i < nk; i++) cv[i].p = kp[i];

    /* Stride small-prime coprimality table */
    static const int SP[] = {3,5,7,11,13,17,19,23,29,31,37,41,43,47};

    /* Header for pair results */
    printf("  %-12s  %-12s  %-10s  %6s  %6s  %-7s  %s\n",
           "anchor_p", "target_p", "score", "rank", "of_N", "pctile", "signal");
    printf("  %-12s  %-12s  %-10s  %6s  %6s  %-7s  %s\n",
           "------------","------------","----------",
           "------","------","-------","-------");

    int total_pairs = 0, top10 = 0, top25 = 0, top50 = 0, below50 = 0;

    for (int ai = 0; ai < nk; ai++) {
        uint64_t pa = kp[ai];
        PLNInfo anchor = compute_pln(pa);

        uint64_t p_lo = (uint64_t)(pa * (1.0 - CV_WIN));
        uint64_t p_hi = (uint64_t)(pa * (1.0 + CV_WIN));
        if (p_lo < 3) p_lo = 3;

        /* Check if any known target falls in this window */
        int n_targets = 0;
        for (int ti = 0; ti < nk; ti++) {
            if (ti == ai) continue;
            if (kp[ti] >= p_lo && kp[ti] <= p_hi) n_targets++;
        }
        if (n_targets == 0) continue; /* skip anchor — no pairs */

        /* Build odd stride coprime to small primes */
        if (p_lo % 2 == 0) p_lo++;
        uint64_t odd_count = (p_hi > p_lo) ? (p_hi - p_lo) / 2 + 1 : 1;
        uint64_t stride = (odd_count / (uint64_t)CV_SAMPLE_N) * 2;
        if (stride < 2) stride = 2;
        bool adj; do { adj = false;
            for (int si = 0; si < 14; si++)
                if (stride % SP[si] == 0) { stride += 2; adj = true; break; }
        } while (adj);

        /* Sample score distribution */
        long double *scores = (long double *)malloc((CV_SAMPLE_N + 1024) * sizeof(long double));
        if (!scores) { fprintf(stderr, "  OOM — skipping anchor %llu\n", (unsigned long long)pa); continue; }
        int ns = 0;
        for (uint64_t p = p_lo; p <= p_hi && ns < CV_SAMPLE_N + 1024 - 1; p += stride) {
            if (!is_prime_u64(p)) continue;
            PLNInfo cand = compute_pln(p);
            scores[ns++] = pln_resonance(&cand, &anchor);
        }

        if (ns == 0) {
            fprintf(stderr, "  WARNING: zero samples for anchor %llu (stride=%llu) — skipping\n",
                    (unsigned long long)pa, (unsigned long long)stride);
            free(scores); continue;
        }

        /* Report each known target in window */
        for (int ti = 0; ti < nk; ti++) {
            if (ti == ai) continue;
            uint64_t pt = kp[ti];
            if (pt < p_lo || pt > p_hi) continue;

            PLNInfo tgt = compute_pln(pt);
            long double tscore = pln_resonance(&tgt, &anchor);

            int n_better = 0;
            for (int k = 0; k < ns; k++)
                if (scores[k] > tscore) n_better++;

            int rank = n_better + 1;
            double pctile = 100.0 * (ns - n_better) / (double)ns;

            const char *tag =
                pctile >= 95.0 ? "TOP 5%  **" :
                pctile >= 90.0 ? "top 10% *" :
                pctile >= 75.0 ? "top 25%" :
                pctile >= 50.0 ? "top 50%" : "below med";

            printf("  %-12llu  %-12llu  %.8Lf  %6d  %6d  %5.1f%%  %s\n",
                   (unsigned long long)pa, (unsigned long long)pt,
                   tscore, rank, ns, pctile, tag);

            /* Accumulate for per-prime stats */
            cv[ti].pctiles[cv[ti].n < CV_MAX_K ? cv[ti].n : CV_MAX_K-1] = pctile;
            if (cv[ti].n < CV_MAX_K) cv[ti].n++;

            total_pairs++;
            if (pctile >= 90.0) top10++;
            if (pctile >= 75.0) top25++;
            if (pctile >= 50.0) top50++;
            else                below50++;
        }
        printf("  [anchor %llu: sampled %d primes from [%llu,%llu] stride=%llu]\n\n",
               (unsigned long long)pa, ns,
               (unsigned long long)p_lo, (unsigned long long)p_hi,
               (unsigned long long)stride);

        free(scores);
    }

    /* ---- Aggregate summary ---- */
    printf("=== AGGREGATE SIGNAL SUMMARY (%d pairs) ===\n", total_pairs);
    printf("  TOP  5%%: %d/%d = %.1f%%  (expected %.1f%% if random)\n",
           top10,  total_pairs, 100.0*top10/total_pairs,  5.0);
    printf("  TOP 10%%: %d/%d = %.1f%%  (expected %.1f%% if random)\n",
           top10,  total_pairs, 100.0*top10/total_pairs, 10.0);
    printf("  TOP 25%%: %d/%d = %.1f%%  (expected %.1f%% if random)\n",
           top25,  total_pairs, 100.0*top25/total_pairs, 25.0);
    printf("  TOP 50%%: %d/%d = %.1f%%  (expected %.1f%% if random)\n",
           top50,  total_pairs, 100.0*top50/total_pairs, 50.0);
    printf("  BELOW 50%%: %d/%d = %.1f%%\n\n", below50, total_pairs, 100.0*below50/total_pairs);

    /* ---- Per-prime stats (sorted by median percentile desc) ---- */
    printf("=== PER-PRIME TARGET STATS (median percentile as target) ===\n");
    printf("  Outliers = primes that PLN consistently mis-ranks.\n");
    printf("  %-12s  %5s  %8s  %8s  %s\n",
           "p", "n_cv", "med_pct", "min_pct", "flag");
    printf("  %-12s  %5s  %8s  %8s  %s\n",
           "------------","-----","--------","--------","----");

    /* collect medians for IQR outlier detection */
    double medians[CV_MAX_K*2]; int nm = 0;
    for (int i = 0; i < nk; i++) {
        if (cv[i].n == 0) continue;
        medians[nm++] = median_d(cv[i].pctiles, cv[i].n);
    }
    /* sort medians ascending for Q1/Q3 */
    for (int i = 1; i < nm; i++) {
        double t = medians[i]; int j = i-1;
        while (j >= 0 && medians[j] > t) { medians[j+1] = medians[j]; j--; }
        medians[j+1] = t;
    }
    double q1 = nm >= 4 ? medians[nm/4] : 0.0;
    double q3 = nm >= 4 ? medians[3*nm/4] : 100.0;
    double iqr = q3 - q1;
    double lo_fence = q1 - 1.5 * iqr;
    double hi_fence = q3 + 1.5 * iqr;

    /* sort primes by median desc for display */
    int idx[CV_MAX_K*2]; for (int i = 0; i < nk; i++) idx[i] = i;
    for (int i = 1; i < nk; i++) {
        int ti = idx[i]; int j = i-1;
        double mi = cv[ti].n > 0 ? median_d(cv[ti].pctiles, cv[ti].n) : -1.0;
        while (j >= 0) {
            int tj = idx[j];
            double mj = cv[tj].n > 0 ? median_d(cv[tj].pctiles, cv[tj].n) : -1.0;
            if (mj >= mi) break;
            idx[j+1] = idx[j]; j--;
        }
        idx[j+1] = ti;
    }

    for (int ii = 0; ii < nk; ii++) {
        int i = idx[ii];
        if (cv[i].n == 0) continue;
        double med = median_d(cv[i].pctiles, cv[i].n);
        double minp = cv[i].pctiles[0];
        for (int j = 1; j < cv[i].n; j++) if (cv[i].pctiles[j] < minp) minp = cv[i].pctiles[j];
        const char *flag =
            med > hi_fence  ? "HIGH OUTLIER (suspiciously good?)" :
            med < lo_fence  ? "LOW OUTLIER  (PLN mismatch)" :
            med >= 75.0     ? "good" :
            med >= 50.0     ? "ok" : "weak";
        printf("  %-12llu  %5d  %7.1f%%  %7.1f%%  %s\n",
               (unsigned long long)cv[i].p, cv[i].n, med, minp, flag);
    }

    printf("\n  IQR fences: [%.1f%%, %.1f%%]  (outlier < %.1f%% or > %.1f%%)\n\n",
           q1, q3, lo_fence, hi_fence);
    printf("=== INTERPRETATION ===\n");
    printf("  Signal present if: TOP 10%% rate >> 10%%, and median per-prime > 50%%\n");
    printf("  LOW OUTLIERs: these primes are consistently mis-ranked by the PLN formula\n");
    printf("    → Consider excluding from anchor selection, or they reveal formula gaps\n");
    printf("  HIGH OUTLIERs: suspiciously high rank — verify data integrity\n");
    printf("  Random baseline: every tier should be ~uniformly filled at 10/25/50%%\n");
}

/* ================================================================
 * GOLDEN ORBIT VALIDATION
 * Prints orbit_r for every known Mersenne prime exponent to test
 * whether they cluster near |R| = 1/phi ≈ 0.61803.
 * ================================================================ */
static void validate_pln(void) {
    static const uint64_t known[] = {
        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 = (int)(sizeof(known)/sizeof(known[0]));

    printf("=== PLN GOLDEN ORBIT VALIDATION ===\n");
    printf("Testing |R| = orbit_r(omega) for all %d known Mersenne prime exponents\n", n);
    printf("Golden orbit hypothesis: Mersenne primes cluster near |R| = 1/phi = 0.61803\n\n");
    printf("  %-12s  %-5s  %-16s  %-8s  %-10s  %s\n",
           "p", "D", "omega", "|R|", "|R|-1/phi", "near?");
    printf("  %-12s  %-5s  %-16s  %-8s  %-10s  %s\n",
           "------------", "-----", "----------------", "--------", "----------", "-----");

    double sum_r = 0.0, sum_sq = 0.0;
    int close = 0;

    for (int i = 0; i < n; i++) {
        PLNInfo m = compute_pln(known[i]);
        double dr_gold = (double)m.r - 0.6180339887;
        sum_r  += (double)m.r;
        sum_sq += (double)(m.r * m.r);
        if (fabsl(m.r - 0.6180339887L) < 0.05L) close++;
        printf("  %-12llu  D%-4d %-16.10Lf  %-8.6Lf  %+.6f  %s\n",
               (unsigned long long)known[i], m.D, m.omega, m.r,
               dr_gold, fabsl(dr_gold) < 0.05 ? "YES" : ".");
    }

    double mean_r  = sum_r / n;
    double mean_sq = sum_sq / n;
    double std_r   = sqrt(mean_sq - mean_r*mean_r);

    printf("\n  Count : %d known Mersenne prime exponents\n", n);
    printf("  Mean |R|            : %.6f  (golden = 0.618034)\n", mean_r);
    printf("  Std dev of |R|      : %.6f  (random would be ~0.3)\n", std_r);
    printf("  Within 0.05 of 1/phi: %d/%d = %.1f%%\n\n", close, n, 100.0*close/n);
    printf("Interpretation:\n");
    printf("  mean |R| ≈ 0.618 AND std << 0.3  →  golden orbit signal is REAL\n");
    printf("  mean |R| far from 0.618 or std ≈ 0.3  →  no clustering, no signal\n");
    printf("  The 0.25-weight golden term in pln_resonance() is justified iff signal is real.\n\n");
}

/* ================================================================
 * USAGE
 * ================================================================ */
static void print_usage(void) {
    printf("Usage:\n");
    printf("  mersenne_sieve --sieve <anchor_p> <radius> <steps> <factor_bits>\n");
    printf("  mersenne_sieve --sieve-file <file> [factor_bits]\n");
    printf("  mersenne_sieve --test-p <p1> [p2 p3 ...]\n");
    printf("  mersenne_sieve --factor-limit <p> <factor_bits>\n");
    printf("  mersenne_sieve --validate-pln\n");
    printf("  mersenne_sieve --cross-validate               (uses all 51 known Mersenne primes)\n");
    printf("  mersenne_sieve --cross-validate <p1> <p2> ... (all 51 + extra exponents)\n");
    printf("  mersenne_sieve --help\n\n");
    printf("Examples:\n");
    printf("  mersenne_sieve --sieve 136279841 2.0 100000 64\n");
    printf("  mersenne_sieve --sieve-file mersenne_candidates.txt 64\n");
    printf("  mersenne_sieve --test-p 272559689 272559697 272559701\n");
    printf("  mersenne_sieve --factor-limit 272559689 64\n\n");
    printf("Pipeline with pln_final:\n");
    printf("  pln --prime-file-candidates M136279841.txt 20 0.005 12\n");
    printf("  mersenne_sieve --sieve-file mersenne_candidates.txt 64\n\n");
    printf("Certainty levels:\n");
    printf("  COMPOSITE(p composite):  2^p-1 is definitely composite (p must be prime)\n");
    printf("  COMPOSITE(factor found): 2^p-1 is definitely composite\n");
    printf("  SURVIVED sieve to 2^B:   no small factor — not a primality proof\n");
    printf("  Lucas-Lehmer (GPU):      100%% deterministic — the only real test\n\n");
    printf("Notes on the PLN framework:\n");
    printf("  PLN scores measure geometric proximity in phi-lattice space.\n");
    printf("  They are heuristic rankings for prioritizing GPU test submissions.\n");
    printf("  They have no known mathematical connection to Mersenne primality.\n");
    printf("  The only valid primality tests here are:\n");
    printf("    - is_prime_u64 (for the exponent p itself)\n");
    printf("    - trial_factor (composite detection only)\n");
    printf("    - Lucas-Lehmer on GPU (definitive, not in this tool)\n");
}

int main(int argc, char *argv[]) {
    printf("MERSENNE CANDIDATE SIEVE — honest certainty edition\n\n");

    if (argc < 2 || strcmp(argv[1],"--help")==0 || strcmp(argv[1],"-h")==0) {
        print_usage(); return 0;
    }

    if (strcmp(argv[1],"--sieve")==0 && argc >= 6) {
        uint64_t anchor_p = strtoull(argv[2], NULL, 10);
        long double radius = strtold(argv[3], NULL);
        int steps = atoi(argv[4]);
        double fbits = atof(argv[5]);
        if (fbits <= 0 || fbits > 63) fbits = 60.0;
        run_sieve(anchor_p, radius, steps, fbits);
        return 0;
    }

    if (strcmp(argv[1],"--sieve-file")==0 && argc >= 3) {
        double fbits = argc > 3 ? atof(argv[3]) : 60.0;
        if (fbits <= 0 || fbits > 63) fbits = 60.0;
        sieve_from_file(argv[2], fbits);
        return 0;
    }

    if (strcmp(argv[1],"--test-p")==0 && argc >= 3) {
        uint64_t ps[64];
        int count = 0;
        for (int i = 2; i < argc && count < 64; i++)
            ps[count++] = strtoull(argv[i], NULL, 10);
        test_specific(ps, count, 54.0);
        return 0;
    }

    if (strcmp(argv[1],"--factor-limit")==0 && argc >= 4) {
        uint64_t p = strtoull(argv[2], NULL, 10);
        double fbits = atof(argv[3]);
        uint64_t ps[1] = {p};
        test_specific(ps, 1, fbits);
        return 0;
    }

    if (strcmp(argv[1],"--validate-pln")==0) {
        validate_pln();
        return 0;
    }

    if (strcmp(argv[1],"--cross-validate")==0) {
        /* Optional extra exponents after the flag; uses all 51 + any extras */
        uint64_t extra[64]; int nextra = 0;
        for (int i = 2; i < argc && nextra < 64; i++) {
            uint64_t p = strtoull(argv[i], NULL, 10);
            if (p > 1) extra[nextra++] = p;
        }
        cross_validate(extra, nextra);
        return 0;
    }

    print_usage(); return 1;
}