/*
 * dn_resonance_checker.c
 * PLN-omega-scaled Dn(r) Mahalanobis scorer for Mersenne prime candidates.
 *
 * Corrects two V3.0 scale failures:
 *   1. r = omega / OMEGA_MAX  (never saturates for known Mersenne range)
 *   2. No exponential energy decay  (no energy death at large p)
 *
 * Uses the 16-dimension V3.0b Fibonacci/prime table.
 * Calibrated against all 52 known Mersenne prime exponents.
 *
 * The 16 D_n(r) components with fixed Omega are a 1-D manifold:
 *   D_n(r) = A_n * |r|^k_n   where r = omega/8 and A_n is constant per dimension.
 * The LOO test measures honest discriminative power of this omega-based signature.
 * Next extension: per-dimension Omega (solved from known primes) or analog coupling.
 *
 * Build:
 *   clang -O3 -std=c99 -D_CRT_SECURE_NO_WARNINGS -o dn_resonance_checker.exe dn_resonance_checker.c
 *
 * Usage:
 *   dn_resonance_checker.exe --calibrate           # D_n stats for 52 known primes
 *   dn_resonance_checker.exe --loo                 # LOO cross-validation (52 trials)
 *   dn_resonance_checker.exe --score 137843969     # full z-score breakdown
 *   echo 137843969 | dn_resonance_checker.exe      # score from stdin
 *   dn_resonance_checker.exe 1 2 3 ...             # score from args
 *   type candidates.txt | dn_resonance_checker.exe # batch score
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdint.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

/* ──────────────────────────────────────────────────────────────────────────── */
/* Constants                                                                    */
/* ──────────────────────────────────────────────────────────────────────────── */
#define PHI         1.61803398874989484820
#define LN_PHI      0.48121182505960344750
#define LN2         0.69314718055994530942
#define SQRT5       2.23606797749978969641
#define GOLDEN_ANG  3.88322207745300172950  /* 2*pi/phi — golden angle, V3.0b irrational dim spacing */
#define OMEGA_MAX   8.0     /* prism: all 52 known Mersenne in omega in [−0.81, 7.57] */
#define NUM_DIMS    16      /* V3.0b extended table */
#define NUM_KNOWN   52
#define LOO_SAMPLES 500
static double LOCAL_BW      = 0.5;  /* omega-space Gaussian BW: Re(s)=1/2 critical strip.
                                      * Override with --bw <value> on the command line. */
static double LOCAL_BW_BETA = 0.10; /* beta-space Gaussian BW (beta = frac(omega)).
                                      * 2D kernel w = exp(-0.5*(dOmega/BW)^2 - 0.5*(dBeta/BW_BETA)^2)
                                      * focuses calibration on training primes with matching phase
                                      * (fixes edge-beta failures: primes with beta near 0 or 1
                                      *  appeared as outliers when 1D calibration averaged over all betas).
                                      * LOO sweep: 1D-only = 33.4%, 2D BW_BETA=0.10 = 18.6% mean rank.
                                      * Override with --bw_beta <value> on the command line. */

/* 16 known Riemann zeta zero imaginary parts (non-trivial zeros of ζ(1/2+iγ))
 * These are the frequencies at which the prime-counting function π(x) oscillates.
 * RH claims all non-trivial zeros lie on Re(s)=1/2 — beta_crit ≈ 1/2 encodes this.
 * BW=0.5 found empirically is the same 1/2: the critical strip half-width.
 * Each dimension n gets one frequency γ_n (a frequency comb spanning 14–67 rad). */
static const double GAMMA_RZ[NUM_DIMS] = {
    14.134725, 21.022039, 25.010857, 30.424876, 32.935061, 37.586176,
    40.918719, 43.327073, 48.005150, 49.773832, 52.970321, 56.446247,
    59.347044, 60.831779, 65.112544, 67.079810
};

/* 16-entry Fibonacci F_1..F_16 */
static const double FIB16[NUM_DIMS] = {
    1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987
};

/* 16-entry primes P_1..P_16 */
static const double PRM16[NUM_DIMS] = {
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53
};

/* 52 known Mersenne prime exponents (M_p = 2^p − 1 is prime) */
static const long long KP[NUM_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
};

/* ──────────────────────────────────────────────────────────────────────────── */
/* XorShift-64 RNG (seeded per LOO trial for reproducibility)                  */
/* ──────────────────────────────────────────────────────────────────────────── */
static uint64_t xr_state = 123456789012345ULL;
static void     xr_seed(uint64_t s) { xr_state = s ? s : 1ULL; }
static uint64_t xr_next(void) {
    xr_state ^= xr_state << 13;
    xr_state ^= xr_state >> 7;
    xr_state ^= xr_state << 17;
    return xr_state;
}

/* ──────────────────────────────────────────────────────────────────────────── */
/* PLN omega: log_phi(log_phi(p * ln2))                                        */
/* Maps all 52 known Mersenne exponents to omega in [−0.81, 7.57].             */
/* ──────────────────────────────────────────────────────────────────────────── */
static double pln_omega(long long p) {
    double inner = log((double)p * LN2) / LN_PHI; /* log_phi(p*ln2) */
    if (inner > 0.0)
        return log(inner) / LN_PHI;
    /* inner <= 0 only for p < 2; treat as degenerate */
    return log(fabs(inner) + 1e-300) / LN_PHI * (inner < 0.0 ? -1.0 : 1.0);
}

/* ──────────────────────────────────────────────────────────────────────────── */
/* gen_fib(x): generalised Fibonacci via Binet formula for real x              */
/*   F_x = (phi^x - cos(pi*x)*phi^{-x}) / sqrt(5)                             */
/* For integer n recovers F_n exactly. The cos(pi*x) term oscillates           */
/* differently for each dimension when x = n + frac(omega), breaking the       */
/* 1-D manifold degeneracy present in the integer-only version.                */
/* ──────────────────────────────────────────────────────────────────────────── */
static double gen_fib(double x) {
    double phi_x    = pow(PHI,  x);
    double phi_neg_x = pow(PHI, -x);
    return (phi_x - cos(M_PI * x) * phi_neg_x) / SQRT5;
}

/* ──────────────────────────────────────────────────────────────────────────── */
/* prismatic_omega(omega_p, n): per-dimension field tension Omega_n             */
/*   Omega_n = 0.5 + 0.5 * sin(omega_p + n * GOLDEN_ANG)                      */
/*                                                                             */
/* Derived from V3.0b prismatic_recursion:                                     */
/*   omega_val = 0.5 + 0.5*sin(lat->time + idx*0.01)                          */
/* We substitute lat->time -> omega_p (PLN coordinate of the candidate) and    */
/* idx*0.01 -> n*GOLDEN_ANG (golden-angle spacing for irrational separation).  */
/*                                                                             */
/* Effect: each dimension n samples a different phase of the sin curve.        */
/* Two candidates differing by Δω at the prism scale get a different           */
/* per-dimension pattern of Omega_n values — genuine multi-dimensionality.     */
/* ──────────────────────────────────────────────────────────────────────────── */
static double prismatic_omega(double omega_p, int n) {
    return 0.5 + 0.5 * sin(omega_p + (double)n * GOLDEN_ANG);
}

/* ──────────────────────────────────────────────────────────────────────────── */
/* compute_Dn: Unified lattice operator |L_n| — emergent-1 / D_v6 framework   */
/*                                                                             */
/*   |L_n| = sqrt((coeff·|r|^k + 1_eff·cos(πβ))² + (1_eff·sin(πβ))²)        */
/*                                                                             */
/*   x      = n + β,  β = frac(ω_p)    — non-integer Fibonacci index         */
/*   F_x    = gen_fib(x)                — Binet generalisation               */
/*   Ω_n    = prismatic_omega(ω_p, n)   — per-dim field tension (V3.0b)      */
/*   φ^{kx} = pow(PHI, k*x)             — φ-recursive scaling               */
/*   1_eff  = 1 + |cos(πβφ)| · ln(P_n)/φ^x — emergent lattice "1"          */
/*                                                                             */
/* Connection to Riemann Hypothesis (see GAMMA_RZ table above):               */
/*   β = frac(ω_p) ≈ 0.5 encodes the critical line Re(s) = 1/2.             */
/*   The optimal calibration bandwidth BW = 0.5 (empirically found via LOO    */
/*   sweep) IS the Riemann half-strip width: the D_v6 framework encodes this  */
/*   with n_crit≈-1/2, beta_crit≈1/2 (from hdgl_analog_v30b.c / HDGL).      */
/*   γ_n as per-dim phase frequencies (tested, see git history) inflated       */
/*   local variance at BW=0.5 — they are preserved in GAMMA_RZ for           */
/*   Riemann resonance visualization (D_v6 dust oscilloscope, horizon spiral). */
/* ──────────────────────────────────────────────────────────────────────────── */
static double compute_Dn(int n, double r, double omega_p) {
    int    i      = n - 1;
    double beta   = omega_p - floor(omega_p);      /* frac(ω_p) — phase, ≈1/2 on critical line */
    double x      = (double)n + beta;               /* non-integer Fibonacci index */
    double k      = (double)(n + 1) / 8.0;
    double F_x    = fabs(gen_fib(x));
    double Omega_n = prismatic_omega(omega_p, n);
    double P_n     = PRM16[i];
    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 embedding: 1_eff · e^{iπβ}  — macro ±1 when β ∈ {0,1}          */
    /* The Riemann zero frequencies γ_n (14–67 rad) live at the CALIBRATION  */
    /* level, not the feature level: BW=0.5 = Re(s)=1/2 (critical strip).    */
    /* γ_n * β oscillates 2–11× per unit β — adding them as per-dim phases   */
    /* inflates local variance at BW=0.5 and degrades LOO (tested: 33.4%→35.6%).  */
    /* GAMMA_RZ is preserved for Riemann-resonance visualization (see D_v6).  */
    double phase_re = one_eff * cos(M_PI * beta);
    double phase_im = one_eff * sin(M_PI * beta);

    /* |L_n| = |scale + 1_eff·e^{iπβ}| */
    double re = scale + phase_re;
    return sqrt(re * re + phase_im * phase_im);
    (void)FIB16[i];
}

/* ──────────────────────────────────────────────────────────────────────────── */
/* 16-dim fingerprint for exponent p                                           */
/* ──────────────────────────────────────────────────────────────────────────── */
static void fingerprint(long long p, double vec[NUM_DIMS]) {
    double om = pln_omega(p);
    double r  = om / OMEGA_MAX;
    for (int n = 1; n <= NUM_DIMS; n++)
        vec[n - 1] = compute_Dn(n, r, om);
}

/* ──────────────────────────────────────────────────────────────────────────── */
/* local_calibrate: weighted mean & std for dimensions, centred on query_omega */
/*                                                                             */
/* w_i = exp(-0.5 * ((query_omega - omega_i) / LOCAL_BW)^2)                   */
/*                                                                             */
/* As the prism travels up the spiral, primes closest in omega-space dominate  */
/* the fit.  Small primes (omega ~ -0.8 .. 2) lose weight entirely when the    */
/* query is near omega=7.5, and vice-versa.  excl_idx = -1 to include all.    */
/* ──────────────────────────────────────────────────────────────────────────── */
static void local_calibrate(double query_omega, int excl_idx,
                            double mean[NUM_DIMS], double stddev[NUM_DIMS]) {
    double wsum[NUM_DIMS], wmean[NUM_DIMS], wvar[NUM_DIMS];
    double wtotal = 0.0;
    double weights[NUM_KNOWN];

    /* 2D Gaussian kernel: weight by proximity in both omega AND beta=frac(omega).
     * This anchors calibration to training primes sharing the same PLN phase,
     * preventing the 1D-average-beta centroid from misrepresenting edge-beta primes. */
    double query_beta = query_omega - floor(query_omega);

    /* Pass 1: compute 2D weights and weighted mean */
    for (int d = 0; d < NUM_DIMS; d++) { wsum[d] = 0.0; }
    for (int i = 0; i < NUM_KNOWN; i++) {
        if (i == excl_idx) { weights[i] = 0.0; continue; }
        double om_i   = pln_omega(KP[i]);
        double beta_i = om_i - floor(om_i);
        double dw = (query_omega - om_i) / LOCAL_BW;
        double db = (query_beta  - beta_i) / LOCAL_BW_BETA;
        weights[i] = exp(-0.5 * (dw * dw + db * db));
        if (weights[i] < 1e-12) weights[i] = 0.0;
        wtotal += weights[i];
        double v[NUM_DIMS];
        fingerprint(KP[i], v);
        for (int d = 0; d < NUM_DIMS; d++) wsum[d] += weights[i] * v[d];
    }
    /* Fallback 1: 2D kernel collapsed to zero — try 1D omega-only kernel */
    if (wtotal < 1e-12) {
        wtotal = 0.0;
        for (int d = 0; d < NUM_DIMS; d++) wsum[d] = 0.0;
        for (int i = 0; i < NUM_KNOWN; i++) {
            if (i == excl_idx) { weights[i] = 0.0; continue; }
            double om_i = pln_omega(KP[i]);
            double dw   = (query_omega - om_i) / LOCAL_BW;
            weights[i]  = exp(-0.5 * dw * dw);
            if (weights[i] < 1e-12) weights[i] = 0.0;
            wtotal += weights[i];
            double v[NUM_DIMS]; fingerprint(KP[i], v);
            for (int d = 0; d < NUM_DIMS; d++) wsum[d] += weights[i] * v[d];
        }
    }
    /* Fallback 2: all weights zero (pathological) — uniform */
    if (wtotal < 1e-12) {
        wtotal = 0.0;
        for (int i = 0; i < NUM_KNOWN; i++) {
            if (i == excl_idx) continue;
            weights[i] = 1.0;
            wtotal += 1.0;
        }
        for (int d = 0; d < NUM_DIMS; d++) {
            wsum[d] = 0.0;
            double v[NUM_DIMS];
            fingerprint(KP[d < NUM_KNOWN ? d : 0], v); /* recompute cleanly below */
        }
        /* Redo sum with uniform weights */
        for (int d = 0; d < NUM_DIMS; d++) wsum[d] = 0.0;
        for (int i = 0; i < NUM_KNOWN; i++) {
            if (i == excl_idx) continue;
            double v[NUM_DIMS]; fingerprint(KP[i], v);
            for (int d = 0; d < NUM_DIMS; d++) wsum[d] += v[d];
        }
    }
    for (int d = 0; d < NUM_DIMS; d++) {
        wmean[d] = wsum[d] / wtotal;
        wvar[d]  = 0.0;
    }

    /* Pass 2: weighted variance */
    for (int i = 0; i < NUM_KNOWN; i++) {
        if (weights[i] < 1e-12) continue;
        double v[NUM_DIMS];
        fingerprint(KP[i], v);
        for (int d = 0; d < NUM_DIMS; d++) {
            double diff = v[d] - wmean[d];
            wvar[d] += weights[i] * diff * diff;
        }
    }
    for (int d = 0; d < NUM_DIMS; d++) {
        mean[d]   = wmean[d];
        /* Match Python: use std=1.0 when variance is near zero (degenerate
         * single-point cluster after 1D/2D fallback), preventing z-score explosion. */
        double var_d = wvar[d] / wtotal;
        stddev[d] = (var_d > 1e-20) ? sqrt(var_d) : 1.0;
    }
}

/* ──────────────────────────────────────────────────────────────────────────── */
/* RMS Mahalanobis distance (diagonal covariance)                              */
/* z[d] = (D_d(p) - mean[d]) / stddev[d];  dist = sqrt(mean(z^2))             */
/* ──────────────────────────────────────────────────────────────────────────── */
static double maha_dist(long long p, const double mean[NUM_DIMS],
                        const double stddev[NUM_DIMS], double zs[NUM_DIMS]) {
    double vec[NUM_DIMS], d2 = 0.0;
    fingerprint(p, vec);
    for (int i = 0; i < NUM_DIMS; i++) {
        double z = (vec[i] - mean[i]) / stddev[i];
        if (zs) zs[i] = z;
        d2 += z * z;
    }
    return sqrt(d2 / NUM_DIMS);
}

/* ──────────────────────────────────────────────────────────────────────────── */
/* --calibrate mode                                                             */
/* ──────────────────────────────────────────────────────────────────────────── */
static void run_calibrate(void) {
    /* For the display, calibrate centered at p52 (our prism center) */
    double om52 = pln_omega(136279841LL);
    double mean[NUM_DIMS], stddev[NUM_DIMS];
    local_calibrate(om52, -1, mean, stddev);

    printf("=== Dn(r) CALIBRATION: locally-weighted, centred at p=136279841 (omega=%.4f) ===\n", om52);
    printf("  Prism:  r = omega / %.1f\n", OMEGA_MAX);
    printf("  omega:  log_phi(log_phi(p * ln2))   range across 52 primes: [-0.81, 7.57]\n");
    printf("  BW:     %.2f omega units (Gaussian proximity weighting)\n", (double)LOCAL_BW);
    printf("  Dims:   %d   (Fibonacci F_1..F_16, Primes P_1..P_16)\n\n", NUM_DIMS);
    printf("  dim   k=(n+1)/8   F_x(p52)     Omega_n(p52)  1_eff(p52)   mean(local)  std(local)   cv\n");
    printf("  ----  ----------  -----------  ------------  ----------  -----------  -----------  --------\n");
    for (int d = 0; d < NUM_DIMS; d++) {
        int    n     = d + 1;
        double k     = (double)(n + 1) / 8.0;
        double beta52 = om52 - floor(om52);
        double x52    = (double)n + beta52;
        double F_x52  = gen_fib(x52);
        if (F_x52 < 0.0) F_x52 = -F_x52;
        double Omn52  = prismatic_omega(om52, n);
        double delta52 = fabs(cos(M_PI * beta52 * PHI)) * log(PRM16[d]) / pow(PHI, x52);
        double one_eff52 = 1.0 + delta52;
        double cv    = mean[d] > 1e-12 ? stddev[d] / mean[d] : 0.0;
        printf("  D%-4d  %10.4f  %11.5f  %12.6f  %10.6f  %11.5f  %11.5f  %8.4f\n",
               n, k, F_x52, Omn52, one_eff52, mean[d], stddev[d], cv);
    }

    printf("\n  %-12s  %8s  %8s  %8s  %8s  %s\n",
           "p", "omega", "r=o/8", "maha_d", "w(p52)", "");
    printf("  %-12s  %8s  %8s  %8s  %8s\n",
           "------------", "--------", "--------", "--------", "--------");
    for (int i = 0; i < NUM_KNOWN; i++) {
        double om  = pln_omega(KP[i]);
        double dw  = (om52 - om) / LOCAL_BW;
        double wi  = exp(-0.5 * dw * dw);
        /* Local stats centred at this prime's own omega for self-distance */
        double lm[NUM_DIMS], ls[NUM_DIMS];
        local_calibrate(om, i, lm, ls);
        double d   = maha_dist(KP[i], lm, ls, NULL);
        printf("  %-12lld  %8.4f  %8.4f  %8.4f  %8.4f%s\n",
               KP[i], om, om / OMEGA_MAX, d, wi,
               d > 2.0 ? "  <-- outlier" : d > 1.5 ? "  marginal" : "");
    }
}

/* ──────────────────────────────────────────────────────────────────────────── */
/* is_known_mersenne: 1 if p is one of the 52 known Mersenne prime exponents   */
/* ──────────────────────────────────────────────────────────────────────────── */
static int is_known_mersenne(long long p) {
    for (int i = 0; i < NUM_KNOWN; i++)
        if (KP[i] == p) return 1;
    return 0;
}

/* ──────────────────────────────────────────────────────────────────────────── */
/* --loo mode                                                                   */
/* ──────────────────────────────────────────────────────────────────────────── */
static void run_loo(void) {
    printf("=== LOO CROSS-VALIDATION: 16-dim PLN Dn(r) fingerprint (local calibration) ===\n");
    printf("  %d random odd integers per trial (seeded from left-out prime)\n", LOO_SAMPLES);
    printf("  Local BW_omega = %.2f, BW_beta = %.2f (2D omega x beta kernel)\n\n",
           (double)LOCAL_BW, (double)LOCAL_BW_BETA);
    printf("  rank%% = fraction of randoms scoring LOWER (closer to local center) than known prime\n");
    printf("  < 50%% means the known prime beats random baseline\n\n");
    printf("  %-12s  %8s  %8s  %6s  %8s  %4s\n",
           "p", "LOO_dist", "w_self", "n_eff", "rank_pct", "pass");
    printf("  %-12s  %8s  %8s  %6s  %8s  %4s\n",
           "------------", "--------", "--------", "------", "--------", "----");

    double sum_rank = 0.0;
    int    n_pass   = 0;
    double band_sum[3]  = {0.0, 0.0, 0.0};  /* primitive / middle / astronomical */
    int    band_cnt[3]  = {0, 0, 0};
    int    band_pass[3] = {0, 0, 0};

    for (int leave = 0; leave < NUM_KNOWN; leave++) {
        long long p   = KP[leave];
        double    om  = pln_omega(p);
        xr_seed((uint64_t)p ^ 0xDEADBEEFCAFEBABEULL);

        /* Local calibration centred at this prime's omega, excluding itself */
        double lm[NUM_DIMS], ls[NUM_DIMS];
        local_calibrate(om, leave, lm, ls);

        /* Effective number of training primes (sum of weights normalised) */
        double wsum = 0.0;
        double w_self = 0.0; /* weight this prime WOULD have had (diagnostic) */
        for (int i = 0; i < NUM_KNOWN; i++) {
            double dw = (om - pln_omega(KP[i])) / LOCAL_BW;
            double wi = exp(-0.5 * dw * dw);
            wsum += (i != leave) ? wi : 0.0;
            if (i == leave) w_self = wi;
        }
        double n_eff = wsum;  /* effective N */

        /* Distance of left-out prime from its local neighborhood */
        double d_left = maha_dist(p, lm, ls, NULL);

        /* Comparison range: symmetric in omega ± LOCAL_BW, mapped back to p */
        long long lo = (p / 2 < 3) ? 3LL : p / 2;
        long long hi = (p * 2 > 500000000LL) ? 500000000LL : p * 2;
        if (hi - lo < 1000LL) { lo = (lo > 500) ? lo - 500 : 3LL; hi = lo + 1000LL; }
        long long range = hi - lo;

        int n_better = 0, n_valid = 0;
        for (int s = 0; s < LOO_SAMPLES; s++) {
            long long rp = lo + (long long)(xr_next() % (uint64_t)range);
            if (rp < 3) rp = 3;
            if (rp % 2 == 0) rp++;
            if (is_known_mersenne(rp)) continue;  /* don't count known Mersenne as random */
            n_valid++;
            if (maha_dist(rp, lm, ls, NULL) < d_left) n_better++;
        }

        double rank_pct = n_valid > 0 ? 100.0 * n_better / n_valid : 50.0;
        sum_rank += rank_pct;
        if (rank_pct < 50.0) n_pass++;
        int band = (om < 4.65) ? 0 : (om < 6.0) ? 1 : 2;
        band_sum[band]  += rank_pct;
        band_cnt[band]++;
        if (rank_pct < 50.0) band_pass[band]++;

        printf("  %-12lld  %8.4f  %8.4f  %6.1f  %7.1f%%  %s\n",
               p, d_left, w_self, n_eff, rank_pct,
               rank_pct < 50.0 ? "YES" : "no");
    }

    double mean_rank = sum_rank / NUM_KNOWN;
    printf("\n  --- LOO SUMMARY ---\n");
    printf("  Mean rank%%:         %.1f%%  (<50%% = better than random)\n", mean_rank);
    printf("  Beat random (YES):  %d / %d  (%.0f%%)\n",
           n_pass, NUM_KNOWN, 100.0 * n_pass / NUM_KNOWN);
    printf("  Verdict: %s\n",
           mean_rank < 35.0 ? "STRONG discriminative power" :
           mean_rank < 45.0 ? "Moderate signal above random" :
           mean_rank < 50.0 ? "Marginal signal above random" :
                              "No discriminative power (at/below random baseline)");
    printf("  Local BW_omega=%.2f, BW_beta=%.2f: 2D Gaussian kernel in omega x beta space.\n",
           (double)LOCAL_BW, (double)LOCAL_BW_BETA);
    static const char *bname[3] = {
        "primitive  p<=127      (omega<4.65)",
        "middle     521-132049  (4.65<=omega<6.0)",
        "astron.    p>=216091   (omega>=6.0)"
    };
    printf("\n  --- STRATIFIED BY OMEGA BAND ---\n");
    for (int b = 0; b < 3; b++) {
        if (band_cnt[b] > 0)
            printf("  %-44s  n=%2d  mean=%5.1f%%  YES=%d/%d\n",
                   bname[b], band_cnt[b],
                   band_sum[b] / band_cnt[b], band_pass[b], band_cnt[b]);
    }
}

/* ──────────────────────────────────────────────────────────────────────────── */
/* --score mode                                                                 */
/* ──────────────────────────────────────────────────────────────────────────── */
static void score_one(long long p, const double mean[NUM_DIMS],
                      const double stddev[NUM_DIMS]) {
    double zs[NUM_DIMS];
    double dist = maha_dist(p, mean, stddev, zs);
    double om   = pln_omega(p);
    printf("Candidate p = %lld\n", p);
    printf("  PLN omega = %.6f   r = %.6f\n", om, om / OMEGA_MAX);
    printf("  Mahalanobis dist (RMS z-score) = %.4f\n", dist);
    printf("  Verdict: %s\n",
           dist < 1.0 ? "RESONANT  (within 1-sigma of known-prime cloud)" :
           dist < 2.0 ? "marginal  (1–2 sigma)" :
                        "distant   (>2 sigma from cloud)");
    printf("  Per-dimension z-scores:\n");
    for (int d = 0; d < NUM_DIMS; d++) {
        printf("    D%-2d  k=%.4f  z=%+7.4f%s\n",
               d + 1, (double)(d + 2) / 8.0, zs[d],
               fabs(zs[d]) > 2.0 ? "  <-- outlier" : "");
    }
}

/* ──────────────────────────────────────────────────────────────────────────── */
/* main                                                                         */
/* ──────────────────────────────────────────────────────────────────────────── */
int main(int argc, char *argv[]) {
    /* Pre-scan for --bw and --bw_beta bandwidth overrides */
    for (int i = 1; i < argc; i++) {
        if (strcmp(argv[i], "--bw") == 0 && i + 1 < argc) {
            LOCAL_BW = atof(argv[i + 1]);
            i++;
        }
        if (strcmp(argv[i], "--bw_beta") == 0 && i + 1 < argc) {
            LOCAL_BW_BETA = atof(argv[i + 1]);
            i++;
        }
    }
    if (argc >= 2 && strcmp(argv[1], "--calibrate") == 0) {
        run_calibrate();
        return 0;
    }
    if (argc >= 2 && strcmp(argv[1], "--loo") == 0) {
        run_loo();
        return 0;
    }
    if (argc >= 3 && strcmp(argv[1], "--score") == 0) {
        long long p  = atoll(argv[2]);
        if (p >= 2) {
            double om = pln_omega(p);
            double mean[NUM_DIMS], stddev[NUM_DIMS];
            local_calibrate(om, -1, mean, stddev);
            printf("  (local calibration centred at omega=%.4f, BW=%.2f)\n",
                   om, (double)LOCAL_BW);
            score_one(p, mean, stddev);
        }
        return 0;
    }

    /* Batch scoring: each prime gets its own local calibration */
    int has_header = 0;
#define PRINT_HDR() if (!has_header) { \
    printf("%-14s  %8s  %8s  %s\n", "p", "dist", "omega", "verdict"); \
    printf("%-14s  %8s  %8s  %s\n", \
           "--------------", "--------", "--------", "-------"); \
    has_header = 1; }

    if (argc > 1) {
        PRINT_HDR();
        for (int i = 1; i < argc; i++) {
            long long p = atoll(argv[i]);
            if (p < 2) continue;
            double om = pln_omega(p);
            double mean[NUM_DIMS], stddev[NUM_DIMS];
            local_calibrate(om, -1, mean, stddev);
            double d  = maha_dist(p, mean, stddev, NULL);
            printf("%-14lld  %8.4f  %8.4f  %s\n", p, d, om,
                   d < 1.0 ? "RESONANT" : d < 2.0 ? "marginal" : "distant");
        }
    } else {
        PRINT_HDR();
        char line[256];
        while (fgets(line, sizeof(line), stdin)) {
            long long p = atoll(line);
            if (p < 2) continue;
            double om = pln_omega(p);
            double mean[NUM_DIMS], stddev[NUM_DIMS];
            local_calibrate(om, -1, mean, stddev);
            double d  = maha_dist(p, mean, stddev, NULL);
            printf("%-14lld  %8.4f  %8.4f  %s\n", p, d, om,
                   d < 1.0 ? "RESONANT" : d < 2.0 ? "marginal" : "distant");
        }
    }
    return 0;
}
