/*
 * bigg_overfit.c — PATH 2: BigG/Fudge10 Overfit Mersenne Scorer
 * ==============================================================
 *
 * DESIGN PATH SEPARATION:
 *   Path 1 (resonance)  = mersenne_sieve.c + pln_final.c  — generalize
 *   Path 2 (this file)  = overfit the BigG/Fudge10 physics — memorize
 *   Path 3 (base4096)   = precision scaling — base4096_hunter.py
 *
 *   DO NOT MIX: this file deliberately over-fits to known Mersenne primes
 *   using the empirically-validated D_n(r) power-law formula from
 *   hdgl_analog_v30.c (the BigG unification project, <1% CODATA error
 *   on h, G, k_B). The resonance path in mersenne_sieve.c must remain
 *   independent and ignorant of these constants.
 *
 * Core idea:
 *   For each candidate exponent p, compute the full HDGL lattice fingerprint
 *   using the EXACT constants from BigG/Fudge10:
 *     - NumericLattice upper/lower fields, analog dims, sibling harmonics
 *     - D_n(r) = sqrt(phi * F_n * 2^n * P_n * omega) * r^k
 *     - Base(inf) seeds
 *   Then compare that fingerprint against a table fitted to all 51 known
 *   Mersenne exponents. The fit is DELIBERATELY overfit — the point is
 *   maximum sensitivity to the known patterns, not generalization.
 *
 * Usage:
 *   bigg_overfit --score <p1> [p2 p3 ...]     Score candidates
 *   bigg_overfit --fit                         Show fit quality on known primes
 *   bigg_overfit --rank <center> <delta_pct> <N>  Rank top N candidates
 *   bigg_overfit --help
 *
 * Compile:
 *   clang -O3 -std=c99 -D_CRT_SECURE_NO_WARNINGS -o bigg_overfit.exe bigg_overfit.c
 */

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

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

/* ================================================================
 * CONSTANTS — V3.0b extended (16 dims), from hdgl_analog_v30b.c
 * ================================================================ */
#define PHI      1.6180339887498948482
#define LN_PHI   0.4812118250596034748
#define LN2      0.6931471805599453094
#define LOG10_2  0.30102999566398119521
#define SQRT5    2.2360679774997896964
#define GOLDEN_ANG 3.8832220774530017295  /* 2*pi/phi — golden angle */

/* V3.0b: expanded to 16 dimensions */
#define NUM_DN 16
static const double FIB_TABLE[NUM_DN] = {
    1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987
};
static const double PRIME_TABLE[NUM_DN] = {
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53
};

/* NumericLattice constants — verbatim from hdgl_analog_v30.c */
static const double UPPER_FIELD[7] = {
    170.6180339887, 150.9442719100, 12.6180339887,
    8.8541019662, 4.2360679775, 3.6180339887, 1.6180339887
};
static const double ANALOG_DIMS[13] = {
    8.3141592654, 7.8541019662, 6.4721359549, 5.6180339887,
    4.8541019662, 3.6180339887, 2.6180339887, 1.6180339887,
    1.0000000000, 7.8541019662, 11.0901699437, 17.9442719100,
    29.0344465435
};
static const double LOWER_FIELD[8] = {
    0.0000000001, 0.0344465435, 0.0557280900, 0.0901699437,
    0.1458980338, 0.2360679775, 0.3819660113, 0.6180339887
};
static const double SIBLING_HARMONICS[8] = {
    0.0901699437, 0.1458980338, 0.2360679775, 0.3090169944,
    0.3819660113, 0.4721359549, 0.6545084972, 0.8729833462
};
static const double BASE_INF_SEEDS[64] = {
    0.6180339887, 1.6180339887, 2.6180339887, 3.6180339887, 4.8541019662,
    5.6180339887, 6.4721359549, 7.8541019662, 8.3141592654, 0.0901699437,
    0.1458980338, 0.2360679775, 0.3090169944, 0.3819660113, 0.4721359549,
    0.6545084972, 0.8729833462, 1.0000000000, 1.2360679775, 1.6180339887,
    2.2360679775, 2.6180339887, 3.1415926535, 3.6180339887, 4.2360679775,
    4.8541019662, 5.6180339887, 6.4721359549, 7.2360679775, 7.8541019662,
    8.6180339887, 9.2360679775, 9.8541019662, 10.6180339887, 11.0901699437,
    11.9442719100, 12.6180339887, 13.6180339887, 14.2360679775, 14.8541019662,
    15.6180339887, 16.4721359549, 17.2360679775, 17.9442719100, 18.6180339887,
    19.2360679775, 19.8541019662, 20.6180339887, 21.0901699437, 21.9442719100,
    22.6180339887, 23.6180339887, 24.2360679775, 24.8541019662, 25.6180339887,
    26.4721359549, 27.2360679775, 27.9442719100, 28.6180339887, 29.0344465435,
    29.6180339887, 30.2360679775, 30.8541019662, 31.6180339887
};

/* ================================================================
 * ALL 52 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_KNOWN (int)(sizeof(ALL_MERSENNE_P)/sizeof(ALL_MERSENNE_P[0]))

/* ================================================================
 * D_n(r) — V3.0b: Binet gen_fib + prismatic per-dim Omega
 *
 * gen_fib(x): generalised Fibonacci via Binet for real x.
 *   F_x = (phi^x - cos(pi*x)*phi^{-x}) / sqrt(5)
 *
 * prismatic_omega(omega_p, n): per-dimension field tension.
 *   Omega_n = 0.5 + 0.5*sin(omega_p + n*GOLDEN_ANG)
 *   From V3.0b prismatic_recursion: omega_val = 0.5+0.5*sin(t+idx*0.01)
 *   We substitute lat->time -> omega_p (PLN coord) and 0.01 -> GOLDEN_ANG
 *   so each dimension samples a unique phase of omega_p.
 *
 * compute_Dn_r(n, r, omega_p):
 *   x = n + frac(omega_p)   — non-integer Fibonacci index
 *   D_n = sqrt(phi * |F_x| * 2^n * P_n * Omega_n) * |r|^k
 *   k = (n+1)/8 (progressive dimensionality, unchanged from V3.0)
 * ================================================================ */
static double gen_fib(double x) {
    return (pow(PHI, x) - cos(M_PI * x) * pow(PHI, -x)) / SQRT5;
}

static double prismatic_omega_dim(double omega_p, int n) {
    return 0.5 + 0.5 * sin(omega_p + (double)n * GOLDEN_ANG);
}

/* Unified lattice operator |L_n| — emergent-1 framework.
 * |L_n| = sqrt((coeff*|r|^k + 1_eff*cos(πβ))² + (1_eff*sin(πβ))²)
 * 1_eff = 1 + |cos(πβφ)| · ln(P_n) / φ^x  — emergent "1" at prism scale.
 * r^k preserved for Mahalanobis geometry; ±1,0 emerge at macro (β→0). */
static double compute_Dn_r(int n, double r, double omega_p) {
    if (n < 1 || n > NUM_DN) return 0.0;
    int    idx    = n - 1;
    double beta   = omega_p - floor(omega_p);      /* frac(ω_p) */
    double x      = (double)n + beta;
    double k      = (double)(n + 1) / 8.0;
    double F_x    = fabs(gen_fib(x));
    double Omega_n = prismatic_omega_dim(omega_p, n);
    double P_n     = PRIME_TABLE[idx];
    double two_n   = pow(2.0, (double)n);
    double phi_kx  = pow(PHI, k * x);
    double coeff   = sqrt(PHI * F_x * two_n * P_n * Omega_n * phi_kx);

    /* Emergent 1: 1_eff = 1 + δ,  δ = |cos(πβφ)| · ln(P_n) / φ^x */
    double delta   = fabs(cos(M_PI * beta * PHI)) * log(P_n) / pow(PHI, x);
    double one_eff = 1.0 + delta;

    /* Magnitude scaling: coeff * |r|^k */
    double scale   = coeff * pow(fabs(r), k);

    /* Phase: 1_eff · e^{iπβ} */
    double phase_re = one_eff * cos(M_PI * beta);
    double phase_im = one_eff * sin(M_PI * beta);

    double re = scale + phase_re;
    return sqrt(re * re + phase_im * phase_im);
}

/* ================================================================
 * PLN COORDINATES (same as mersenne_sieve.c — for position only)
 * ================================================================ */
typedef struct {
    uint64_t p;
    double   omega;
    int      D;          /* 1–8 */
    double   theta;      /* frac(omega) */
    double   r;          /* orbit radius */
} PLNCoord;

static double orbit_r(double omega) {
    double frac = omega - floor(omega);
    double theta = frac * 6.2831853071795864769;
    double D_cont = fmod(omega, 8.0) + 1.0;
    double lat = PHI*PHI*(1.0 - pow(PHI, -D_cont)) + 1.0/D_cont;
    if (lat < 1.0) lat = 1.0;
    if (lat > PHI*PHI) lat = PHI*PHI;
    double Rre = cos(theta) - 1.0 + lat;
    double Rim = sin(theta);
    return sqrt(Rre*Rre + Rim*Rim);
}

static PLNCoord compute_pln(uint64_t p) {
    PLNCoord m;
    m.p = p;
    double log2x = (double)p;
    double lnx = log2x * LN2;
    double arg = lnx / LN_PHI;
    m.omega = (arg > 1.0) ? log(arg) / LN_PHI : 0.0;
    m.D     = (int)fmod(floor(m.omega), 8.0) + 1;
    m.theta = fmod(m.omega, 1.0);
    m.r     = orbit_r(m.omega);
    return m;
}

/* ================================================================
 * FAST 64-BIT PRIMALITY (same proven Miller-Rabin)
 * ================================================================ */
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;
}
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;
    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;
}

/* ================================================================
 * BIGG FINGERPRINT
 *
 * The overfit fingerprint for a candidate exponent p.
 * Uses all 8 D_n values, lattice projections, sibling harmonics,
 * and base-infinity seed correlations.
 *
 * Returns a feature vector of fixed length.
 * ================================================================ */
#define FP_LEN 48   /* 16 D_n dims + 32 number-theoretic/lattice features */

typedef struct {
    uint64_t p;
    double   features[FP_LEN];
    double   score;      /* overfit score (higher = more Mersenne-like) */
} BigGFingerprint;

static BigGFingerprint compute_fingerprint(uint64_t p) {
    BigGFingerprint fp;
    memset(&fp, 0, sizeof(fp));
    fp.p = p;

    PLNCoord pln = compute_pln(p);
    int fi = 0;

    /* Feature 0-15: D_n(r) for all 16 dimensions (V3.0b extended)
     * Each uses prismatic_omega_dim(omega, n) and gen_fib(n+frac(omega))
     * to produce genuinely independent per-dimension values. */
    for (int n = 1; n <= NUM_DN && fi < FP_LEN; n++) {
        fp.features[fi++] = compute_Dn_r(n, pln.r, pln.omega);
    }

    /* Feature 8: omega */
    fp.features[fi++] = pln.omega;

    /* Feature 9: theta (fractional omega) */
    fp.features[fi++] = pln.theta;

    /* Feature 10: orbit_r */
    fp.features[fi++] = pln.r;

    /* Feature 11: D (dimension 1-8) */
    fp.features[fi++] = (double)pln.D;

    /* Feature 12–14: lattice projections onto upper field */
    fp.features[fi++] = fmod(pln.omega * UPPER_FIELD[0], 1.0);
    fp.features[fi++] = fmod(pln.omega * UPPER_FIELD[6], 1.0);
    fp.features[fi++] = fmod(pln.r * UPPER_FIELD[3], 1.0);

    /* Feature 15–17: lower field correlations */
    fp.features[fi++] = fmod(pln.theta + LOWER_FIELD[7], 1.0);  /* +0.618 wrap */
    fp.features[fi++] = fmod(pln.omega * LOWER_FIELD[5], 1.0);  /* *0.236 wrap */
    fp.features[fi++] = fabs(sin(pln.omega * LOWER_FIELD[6] * M_PI));

    /* Feature 18–21: sibling harmonic projections */
    for (int i = 0; i < 4 && fi < FP_LEN; i++) {
        fp.features[fi++] = fabs(cos(pln.omega * SIBLING_HARMONICS[i*2] * M_PI));
    }

    /* Feature 22–25: analog dimension phase locks */
    for (int i = 0; i < 4 && fi < FP_LEN; i++) {
        double phase = fmod(pln.omega, ANALOG_DIMS[i]);
        fp.features[fi++] = phase / ANALOG_DIMS[i]; /* normalize 0..1 */
    }

    /* Feature 26–29: base-inf seed distances (nearest 4) */
    double min4[4] = {1e30, 1e30, 1e30, 1e30};
    for (int s = 0; s < 64; s++) {
        double dist = fabs(fmod(pln.omega, 32.0) - BASE_INF_SEEDS[s]);
        /* insert into sorted min4 */
        for (int j = 0; j < 4; j++) {
            if (dist < min4[j]) {
                for (int k = 3; k > j; k--) min4[k] = min4[k-1];
                min4[j] = dist;
                break;
            }
        }
    }
    for (int i = 0; i < 4 && fi < FP_LEN; i++) {
        fp.features[fi++] = min4[i];
    }

    /* Feature 30: D_n cross-product ratio (D_1*D_8 / D_4*D_5) */
    double d1 = compute_Dn_r(1, pln.r, pln.omega);
    double d8 = compute_Dn_r(8, pln.r, pln.omega);
    double d4 = compute_Dn_r(4, pln.r, pln.omega);
    double d5 = compute_Dn_r(5, pln.r, pln.omega);
    fp.features[fi++] = (d4*d5 > 1e-15) ? (d1*d8)/(d4*d5) : 0.0;

    /* Feature 31: golden ratio resonance depth */
    double golden_dist = fabs(pln.r - 1.0/PHI);
    fp.features[fi++] = exp(-golden_dist * golden_dist * 10.0);

    /* Features 32–39: rapidly-varying number-theoretic fingerprint
     * These change at per-candidate scale even for large p, providing
     * fine-grained discrimination the PLN features cannot. */
    fp.features[fi++] = fmod((double)p * PHI, 1.0);           /* 32: golden residue */
    fp.features[fi++] = fmod((double)p * PHI * PHI, 1.0);     /* 33: phi^2 residue */
    fp.features[fi++] = fmod((double)p * 1.41421356237, 1.0); /* 34: sqrt(2) residue */
    fp.features[fi++] = fmod((double)p * M_PI, 1.0);          /* 35: pi residue */
    fp.features[fi++] = (double)(p % 6) / 6.0;                /* 36: sextic residue */
    fp.features[fi++] = (double)(p % 8) / 8.0;                /* 37: octal residue */
    fp.features[fi++] = sqrt((double)p) - floor(sqrt((double)p)); /* 38: frac(sqrt(p)) */
    fp.features[fi++] = log2((double)p) - floor(log2((double)p)); /* 39: frac(log2(p)) */

    return fp;
}

/* ================================================================
 * OVERFIT SCORING — Nearest-Neighbor Kernel Density Estimator
 *
 * Store fingerprints of all 52 known Mersenne primes.
 * Score(candidate) = max similarity to ANY known prime's fingerprint.
 * This is true memorization — a candidate can only score high if it
 * closely matches a specific known prime, not the global average.
 *
 * Per-feature normalization uses global std across known primes.
 * Features with near-zero variance are damped (they don't discriminate
 * between known primes and can't discriminate candidates either).
 *
 * This is DELIBERATELY overfit. No regularization. No held-out set.
 * ================================================================ */

typedef struct {
    BigGFingerprint fps[N_KNOWN];   /* all 52 known fingerprints */
    double feat_std[FP_LEN];       /* per-feature std for normalization */
    double feat_weight[FP_LEN];    /* per-feature weight (inverse std, damped) */
} OverfitModel;

static OverfitModel fit_model(void) {
    OverfitModel model;

    /* Compute fingerprints for all known primes */
    for (int i = 0; i < N_KNOWN; i++) {
        model.fps[i] = compute_fingerprint(ALL_MERSENNE_P[i]);
    }

    /* Per-feature std for normalization */
    for (int f = 0; f < FP_LEN; f++) {
        double sum = 0, sum2 = 0;
        for (int i = 0; i < N_KNOWN; i++) {
            sum += model.fps[i].features[f];
            sum2 += model.fps[i].features[f] * model.fps[i].features[f];
        }
        double mean = sum / N_KNOWN;
        double var = sum2 / N_KNOWN - mean * mean;
        model.feat_std[f] = (var > 1e-20) ? sqrt(var) : 1.0;

        /* Weight: features with very low variance among known primes
         * get damped — they're constant for all primes and can't
         * discriminate. Features with moderate variance are useful. */
        if (var < 1e-10) {
            model.feat_weight[f] = 0.0;  /* dead feature */
        } else {
            model.feat_weight[f] = 1.0 / model.feat_std[f];
        }
    }

    /* Normalize weights to sum to 1 */
    double wsum = 0;
    for (int f = 0; f < FP_LEN; f++) wsum += model.feat_weight[f];
    if (wsum > 0) {
        for (int f = 0; f < FP_LEN; f++) model.feat_weight[f] /= wsum;
    }

    return model;
}

static double score_candidate(const OverfitModel *model, const BigGFingerprint *fp) {
    /* Nearest-neighbor: find max similarity to any known Mersenne prime.
     * Similarity = exp(-weighted_distance^2 / bandwidth)
     * bandwidth = 2.0 makes scores drop off meaningfully */
    double max_sim = 0.0;
    for (int i = 0; i < N_KNOWN; i++) {
        double wdist2 = 0.0;
        for (int f = 0; f < FP_LEN; f++) {
            if (model->feat_weight[f] < 1e-15) continue;
            double diff = fp->features[f] - model->fps[i].features[f];
            double z = diff / model->feat_std[f];
            wdist2 += model->feat_weight[f] * z * z;
        }
        double sim = exp(-wdist2 / 2.0);
        if (sim > max_sim) max_sim = sim;
    }
    return max_sim;
}

/* ================================================================
 * COMMANDS
 * ================================================================ */

static void cmd_fit(void) {
    printf("=== BigG OVERFIT FIT QUALITY (all %d known Mersenne primes) ===\n\n", N_KNOWN);

    OverfitModel model = fit_model();

    printf("  %-12s  D   omega          |R|       score\n", "p");
    printf("  %-12s  --  ---------------  --------  ---------\n", "------------");

    double scores[N_KNOWN];
    for (int i = 0; i < N_KNOWN; i++) {
        BigGFingerprint fp = compute_fingerprint(ALL_MERSENNE_P[i]);
        fp.score = score_candidate(&model, &fp);
        scores[i] = fp.score;
        PLNCoord pln = compute_pln(ALL_MERSENNE_P[i]);
        printf("  %-12llu  D%d  %15.10f  %8.6f  %9.6f\n",
               (unsigned long long)ALL_MERSENNE_P[i],
               pln.D, pln.omega, pln.r, fp.score);
    }

    /* Stats */
    double sum = 0, min_s = 1e30, max_s = -1e30;
    for (int i = 0; i < N_KNOWN; i++) {
        sum += scores[i];
        if (scores[i] < min_s) min_s = scores[i];
        if (scores[i] > max_s) max_s = scores[i];
    }
    printf("\n  Mean score: %.6f\n", sum / N_KNOWN);
    printf("  Min  score: %.6f\n", min_s);
    printf("  Max  score: %.6f\n", max_s);
    printf("  Range     : %.6f\n", max_s - min_s);

    /* Feature importance */
    printf("\n  Feature weights (top 10 most discriminating):\n");
    int idx[FP_LEN]; for (int i = 0; i < FP_LEN; i++) idx[i] = i;
    for (int i = 0; i < FP_LEN - 1; i++)
        for (int j = i + 1; j < FP_LEN; j++)
            if (model.feat_weight[idx[j]] > model.feat_weight[idx[i]])
                { int t = idx[i]; idx[i] = idx[j]; idx[j] = t; }
    static const char *feat_names[FP_LEN] = {
        "D1(r)","D2(r)","D3(r)","D4(r)","D5(r)","D6(r)","D7(r)","D8(r)",
        "omega","theta","orbit_r","D_dim",
        "upper_0_proj","upper_6_proj","upper_3_proj",
        "lower_7_phase","lower_5_phase","lower_6_sin",
        "sib_harm_0","sib_harm_2","sib_harm_4","sib_harm_6",
        "analog_0_lock","analog_1_lock","analog_2_lock","analog_3_lock",
        "inf_near_0","inf_near_1","inf_near_2","inf_near_3",
        "Dn_cross_ratio","golden_depth",
        "golden_resid","phi2_resid","sqrt2_resid","pi_resid",
        "mod6","mod8","frac_sqrt","frac_log2"
    };
    for (int i = 0; i < 10 && i < FP_LEN; i++) {
        int fi = idx[i];
        printf("    %2d. %-18s  weight=%.4f  std=%.6f\n",
               i+1, feat_names[fi], model.feat_weight[fi],
               model.feat_std[fi]);
    }
}

static void cmd_score(uint64_t *exponents, int n) {
    OverfitModel model = fit_model();

    printf("=== BigG OVERFIT SCORES ===\n\n");
    printf("  %-12s  prime?  D   omega          |R|       score    known?\n", "p");
    printf("  %-12s  ------  --  ---------------  --------  -------  ------\n", "------------");

    for (int i = 0; i < n; i++) {
        uint64_t p = exponents[i];
        bool prime = is_prime_u64(p);
        BigGFingerprint fp = compute_fingerprint(p);
        fp.score = score_candidate(&model, &fp);
        PLNCoord pln = compute_pln(p);

        bool known = false;
        for (int k = 0; k < N_KNOWN; k++)
            if (ALL_MERSENNE_P[k] == p) { known = true; break; }

        printf("  %-12llu  %-6s  D%d  %15.10f  %8.6f  %7.4f  %s\n",
               (unsigned long long)p,
               prime ? "yes" : "NO",
               pln.D, pln.omega, pln.r, fp.score,
               known ? "KNOWN" : "");
    }
}

static void cmd_rank(uint64_t center, double delta_pct, int top_n) {
    OverfitModel model = fit_model();

    /* Compute score for center as reference */
    BigGFingerprint center_fp = compute_fingerprint(center);
    double center_score = score_candidate(&model, &center_fp);

    uint64_t p_lo = (uint64_t)(center * (1.0 - delta_pct / 100.0));
    uint64_t p_hi = (uint64_t)(center * (1.0 + delta_pct / 100.0));
    if (p_lo < 3) p_lo = 3;
    if (p_lo % 2 == 0) p_lo++;

    printf("=== BigG OVERFIT RANKING ===\n");
    printf("  Center: %llu   Window: ±%.1f%%  [%llu, %llu]\n",
           (unsigned long long)center, delta_pct,
           (unsigned long long)p_lo, (unsigned long long)p_hi);
    printf("  Center score: %.6f\n", center_score);
    printf("  Scanning for top %d candidates...\n\n", top_n);

    /* Heap: keep top N */
    typedef struct { uint64_t p; double score; } Ranked;
    Ranked *heap = calloc(top_n + 1, sizeof(Ranked));
    int heap_n = 0;

    /* Stride — sample ~200K primes */
    uint64_t range = p_hi - p_lo;
    uint64_t stride = range / 200000;
    if (stride < 2) stride = 2;
    if (stride % 2 != 0) stride++;
    /* coprime to small primes */
    static const int SP[] = {3,5,7,11,13,17,19,23,29,31,37,41,43,47};
    bool adj; do { adj = false;
        for (int si = 0; si < 14; si++)
            if (stride % SP[si] == 0) { stride += 2; adj = true; break; }
    } while (adj);

    int tested = 0;
    for (uint64_t p = p_lo; p <= p_hi; p += stride) {
        if (!is_prime_u64(p)) continue;

        BigGFingerprint fp = compute_fingerprint(p);
        fp.score = score_candidate(&model, &fp);
        tested++;

        /* Min-heap insert */
        if (heap_n < top_n) {
            heap[heap_n] = (Ranked){p, fp.score};
            heap_n++;
            /* bubble up is unnecessary for small N, insertion sort at end */
        } else if (fp.score > heap[0].score) {
            /* Find minimum in heap (simple linear for small top_n) */
            int min_i = 0;
            for (int j = 1; j < heap_n; j++)
                if (heap[j].score < heap[min_i].score) min_i = j;
            heap[min_i] = (Ranked){p, fp.score};
        }
    }

    /* Sort heap descending */
    for (int i = 0; i < heap_n - 1; i++)
        for (int j = i + 1; j < heap_n; j++)
            if (heap[j].score > heap[i].score)
                { Ranked t = heap[i]; heap[i] = heap[j]; heap[j] = t; }

    printf("  Tested %d prime exponents (stride=%llu)\n\n",
           tested, (unsigned long long)stride);
    printf("  %-4s  %-12s  D   omega          |R|       score\n", "rank", "p");
    printf("  %-4s  %-12s  --  ---------------  --------  -------\n", "----", "------------");

    for (int i = 0; i < heap_n; i++) {
        PLNCoord pln = compute_pln(heap[i].p);
        bool known = false;
        for (int k = 0; k < N_KNOWN; k++)
            if (ALL_MERSENNE_P[k] == heap[i].p) { known = true; break; }
        printf("  %-4d  %-12llu  D%d  %15.10f  %8.6f  %7.4f%s\n",
               i + 1, (unsigned long long)heap[i].p,
               pln.D, pln.omega, pln.r, heap[i].score,
               known ? "  ** KNOWN MERSENNE **" : "");
    }

    free(heap);
}

/* ================================================================
 * MAIN
 * ================================================================ */
static void print_usage(void) {
    printf("BigG Overfit Mersenne Scorer — PATH 2 (deliberate overfit)\n\n");
    printf("Usage:\n");
    printf("  bigg_overfit --fit                              Show fit on all 52 known primes\n");
    printf("  bigg_overfit --score <p1> [p2 p3 ...]           Score specific exponents\n");
    printf("  bigg_overfit --rank <center> <delta%%> <top_N>   Rank top N in ±delta%% window\n");
    printf("  bigg_overfit --help\n\n");
    printf("This tool uses the BigG/Fudge10 D_n(r) power-law formula to deliberately\n");
    printf("overfit to known Mersenne prime structure. NOT a primality test.\n");
    printf("Keep separate from Path 1 (resonance) — see design-paths.md.\n");
}

int main(int argc, char *argv[]) {
    if (argc < 2) { print_usage(); return 1; }

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

    if (strcmp(argv[1], "--fit") == 0) {
        cmd_fit();
        return 0;
    }

    if (strcmp(argv[1], "--score") == 0 && argc >= 3) {
        int cap = argc - 2 + 1;
        uint64_t *exponents = malloc(cap * sizeof(uint64_t));
        int n = 0;
        for (int i = 2; i < argc; i++) {
            uint64_t p = strtoull(argv[i], NULL, 10);
            if (p > 1) exponents[n++] = p;
        }
        cmd_score(exponents, n);
        free(exponents);
        return 0;
    }

    if (strcmp(argv[1], "--score-file") == 0 && argc >= 3) {
        FILE *f = fopen(argv[2], "r");
        if (!f) { fprintf(stderr, "Cannot open file: %s\n", argv[2]); return 1; }
        int cap = 16384, n = 0;
        uint64_t *exponents = malloc(cap * sizeof(uint64_t));
        char line[64];
        while (fgets(line, sizeof(line), f)) {
            char *s = line;
            while (*s == ' ' || *s == '\t') s++;
            if (*s == '#' || *s == '\n' || *s == '\r' || *s == '\0') continue;
            uint64_t p = strtoull(s, NULL, 10);
            if (p > 1) {
                if (n >= cap) { cap *= 2; exponents = realloc(exponents, cap * sizeof(uint64_t)); }
                exponents[n++] = p;
            }
        }
        fclose(f);
        cmd_score(exponents, n);
        free(exponents);
        return 0;
    }

    if (strcmp(argv[1], "--rank") == 0 && argc >= 5) {
        uint64_t center = strtoull(argv[2], NULL, 10);
        double delta = atof(argv[3]);
        int top_n = atoi(argv[4]);
        if (center < 3 || delta <= 0 || top_n < 1) {
            printf("Invalid arguments. center > 2, delta > 0, top_N >= 1\n");
            return 1;
        }
        if (top_n > 1000) top_n = 1000;
        cmd_rank(center, delta, top_n);
        return 0;
    }

    print_usage();
    return 1;
}
