// ============================================================================
// HDGL Continuous Mersenne Sieve v34
// ============================================================================
//
// [ROADMAP ITEM 15] Eliminates the discrete candidate buffer.
//
// Architecture: SIEVE replaces EMITTER/CANDIDATE model.
//   Each slot i is permanently assigned a Mersenne exponent p_i from a band:
//       p_i = BASE_P + (i % P_BAND_SIZE) * P_STRIDE
//
//   The slot evolves its LL proxy state toward convergence for p_i.
//   When convergence is detected (LL residue near zero mod (2^p_i − 1)),
//   the slot writes p_i to a ring buffer d_prime_found[].
//
//   The host drains d_prime_found each cycle — no candidate promotion logic,
//   no score accumulator, no threshold tuning.
//
// Sieve kernel replaces hdgl_field_step_kernel:
//   - Runs 4× LL proxy steps per call (same cadence as v32 LL-lite)
//   - GRA plasticity steers r_harmonic toward Mersenne-consistent oscillation
//   - Wave coupling remains: neighbours with close exponents help each other converge
//
// Sieve layout:
//   d_assigned_p[N]    — const once assigned, exponent for each slot
//   d_sieve_state[N]   — 32-bit LL proxy state (evolves each kernel call)
//   d_sieve_r_h[N]     — r_harmonic (plasticity)
//   d_prime_found[256] — ring buffer of found exponents
//   d_prime_count      — atomic write index into ring buffer
//
// Compile:
//   nvcc -O3 -arch=sm_86 -lineinfo hdgl_sieve_v34.cu hdgl_host_v34.c \
//        hdgl_critic_v33.c -o hdgl_v34 -lm
// ============================================================================

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

// ============================================================================
// Sieve configuration constants
// ============================================================================

// Confirmed small Mersenne prime exponents within warp-LL verification range
// (LIMBS=64 → max 3904 bits, so exponents ≤ 3904 are handled exactly).
// Slots divided into N_SMALL_MP equal bands, each centred on one of these.
#define N_SMALL_MP  8
// Host-accessible copy for alloc/harvest functions
static const uint32_t SMALL_MP_HOST[N_SMALL_MP] = {
    61u, 89u, 107u, 127u, 521u, 607u, 1279u, 3217u
};
// Device copy (kernel-accessible via __constant__ cache)
__device__ __constant__ uint32_t SMALL_MP[N_SMALL_MP] = {
    61u, 89u, 107u, 127u, 521u, 607u, 1279u, 3217u
};
#define N_KNOWN_MP  N_SMALL_MP

#define SIEVE_BLOCK          256
#define SIEVE_LL_ITERS       8
#define SIEVE_PRIME_RING     256
#define SIEVE_RESIDUE_EPS    0.001f
#define SIEVE_CAND_THRESH    4.0f     // acc threshold to emit a warp-LL Candidate
#define MAX_SIEVE_CANDS      256      // must match MAX_WARP_CANDS in warp_ll
#define SIEVE_GRA_PLASTICITY 0.005f
#define SIEVE_R_H_MAX        1000.0f
#define SIEVE_REWARD_DECAY   0.95f

// ============================================================================
// BigG + Fudge10 empirical scaling parameters
// Fitted to 1000+ Type Ia Pan-STARRS1 supernovae; validated against 200+ CODATA.
//   k      = 1.049342  — emergent power-law coupling exponent
//   s0     = 0.994533  — entropy base (controls GRA plasticity damping)
//   alpha  = 0.340052  — Omega evolution exponent (field phase velocity scale)
//   Omega0 = 1.049675  — base oscillation (dt scaling)
// ============================================================================
#define BIGG_K      1.049342f   // power-law decay exponent for proximity coupling
#define BIGG_S0     0.994533f   // entropy base for GRA plasticity damping
#define BIGG_ALPHA  0.340052f   // Omega evolution: dt scales as Omega0 / (1+r_h)^alpha
#define BIGG_OMEGA0 1.049675f   // base oscillation

// ============================================================================
// Candidate struct — must match layout in ALL translation units
// ============================================================================
typedef struct {
    int   slot_idx;
    float score;
    float ll_seed_f;
    float r_harmonic;   // encodes Mersenne exponent p for sieve candidates
    float phase;
    float coherence;
    float amp;
} Candidate;

// CSSState forward declaration (full definition in hdgl_warp_ll_v33.cu)
typedef struct { uint64_t S[64]; uint64_t C[64]; } CSSState;

// Warp-LL host launcher forward declarations
extern "C" {
void hdgl_warp_ll_v33_alloc(CSSState **out);
void hdgl_warp_ll_v33_free(CSSState *p);
void hdgl_warp_ll_v33_seed(CSSState *d, int n, cudaStream_t s);
void hdgl_warp_ll_v33_launch(const Candidate *d_cands, int n_cands,
    int p_bits, int iters, CSSState *d_css,
    float *d_residue, int8_t *d_verified,
    float eps, cudaStream_t s);
}

// ============================================================================
// Device arrays — held as global __device__ pointers and scalars
// ============================================================================

typedef struct {
    uint32_t  *assigned_p;    // [N] — Mersenne exponent for each slot
    uint32_t  *sieve_state;   // [N] — LL proxy state
    float     *r_harmonic;    // [N] — GRA plasticity state
    float     *phase;         // [N] — phase for wave coupling
    float     *reward_accum;  // [N] — convergence confidence accumulator
    uint32_t  *prime_found;   // [SIEVE_PRIME_RING] — proxy-filtered ring
    int       *prime_count;   // atomic counter into proxy ring
    Candidate *cand_buffer;   // [MAX_SIEVE_CANDS] — high-conf candidates for warp LL
    int       *cand_count;    // atomic counter into cand_buffer
    int        N;
    int        S;             // lattice row width (for neighbour calculation)
    float      dt;
} SieveDevState;

__device__ __constant__ SieveDevState g_sieve;

// ============================================================================
// LL proxy arithmetic
// ============================================================================

__device__ __forceinline__ uint32_t ll_step32_mod(uint32_t s, uint32_t p) {
    // s → (s² − 2) mod p  — actual per-exponent reduction.
    // True LL needs mod (2^p−1) which requires arbitrary precision for p>32.
    // Reducing mod p (the exponent itself) makes each slot's orbit unique and
    // provides a fast proxy signal: orbits cluster differently near prime p.
    uint64_t sq = (uint64_t)s * (uint64_t)s;
    uint64_t r  = sq % (uint64_t)p;
    return (uint32_t)((r >= 2ULL) ? r - 2ULL : r + (uint64_t)p - 2ULL);
}

__device__ __forceinline__ float sieve_residue(uint32_t s) {
    // Normalised residue ∈ [0,1]: how close s is to 0 mod (2^32-1)
    float r = (float)s * 2.32830644e-10f;  // / (2^32)
    if (r > 0.5f) r = 1.0f - r;           // fold to [0, 0.5]
    return r * 2.0f;                       // map to [0, 1]
}

// ============================================================================
// Sieve field step kernel
//
// Each slot:
//   1. Runs SIEVE_LL_ITERS LL proxy steps for its assigned exponent p_i
//   2. Measures convergence: residue < SIEVE_RESIDUE_EPS
//   3. Wave coupling: phase synchronises toward neighbours with close exponents
//   4. GRA plasticity: r_harmonic increases when convergence improves
//   5. If converged: write p_i to prime_found ring (atomic, deduplicate)
// ============================================================================

__global__ void hdgl_sieve_kernel(void)
{
    const int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i >= g_sieve.N) return;

    uint32_t p_i    = g_sieve.assigned_p[i];
    uint32_t s      = g_sieve.sieve_state[i];
    float    r_h    = g_sieve.r_harmonic[i];
    float    ph     = g_sieve.phase[i];
    float    acc    = g_sieve.reward_accum[i];
    int      N      = g_sieve.N;
    int      S      = g_sieve.S;

    // LL proxy evolution: run SIEVE_LL_ITERS steps
    #pragma unroll
    for (int k = 0; k < SIEVE_LL_ITERS; k++) {
        s = ll_step32_mod(s, (uint32_t)(p_i & 0xFFFFFFFFu));
    }
    float residue = sieve_residue(s);

    // Convergence confidence
    float delta_acc = (residue < SIEVE_RESIDUE_EPS) ? (1.0f - residue) : -residue;
    acc = SIEVE_REWARD_DECAY * acc + delta_acc;

    // GRA plasticity — damped by BigG entropy s0=0.994533:
    // At each step r_h is multiplied by s0 (entropy decay) before the update,
    // giving a scale-dependent relaxation that mirrors the BigG cosmological fit.
    r_h *= BIGG_S0;
    float plastic = SIEVE_GRA_PLASTICITY * delta_acc;
    r_h += plastic;
    r_h = fmaxf(1.0f, fminf(SIEVE_R_H_MAX, r_h));

    // Wave coupling: 4-neighbour phase synchronisation weighted by exponent closeness
    int ni[4] = {
        (i - 1 + N) % N,
        (i + 1)     % N,
        (i - S + N) % N,
        (i + S)     % N
    };

    float sum_sin = 0.0f;
    #pragma unroll
    for (int j = 0; j < 4; j++) {
        int    nj     = ni[j];
        float  dphi   = g_sieve.phase[nj] - ph;
        // BigG k=1.049342 power-law coupling: prox ~ (Δp + 1)^{-k}
        // Replaces flat exponential — empirically derived from Pan-STARRS1 SN fit.
        uint32_t pj    = g_sieve.assigned_p[nj];
        uint32_t pdiff = (p_i > pj) ? (p_i - pj) : (pj - p_i);
        float    prox  = powf((float)pdiff + 1.0f, -BIGG_K);
        sum_sin += prox * sinf(dphi);
    }

    // Phase Euler step — Omega-scaled by BigG alpha parameter:
    // Omega(r_h) = Omega0 / (1 + r_h/R_MAX)^alpha  (soft scale-dependence)
    float omega_scale = BIGG_OMEGA0 / powf(1.0f + r_h / SIEVE_R_H_MAX, BIGG_ALPHA);
    ph += g_sieve.dt * omega_scale * (1.0f + 0.5f * sum_sin + 0.01f * r_h * residue);
    ph  = fmodf(ph, 6.28318530717958647f);
    if (ph < 0.0f) ph += 6.28318530717958647f;

    // Convergence detection: emit once when converged and acc crosses threshold.
    // Guard: only emit Candidate for warp-LL if p_i is one of the 8 known target
    // exponents — prevents off-band slots (e.g. p=5 neighbours in the M61 band)
    // from producing false confirmations.
    if (residue < SIEVE_RESIDUE_EPS && acc > 2.0f) {
        int pos = atomicAdd(g_sieve.prime_count, 1) % SIEVE_PRIME_RING;
        g_sieve.prime_found[pos] = p_i;

        // Warp-LL candidate — only for the 8 known target exponents
        bool is_target = false;
        for (int tk = 0; tk < N_SMALL_MP; tk++) {
            if (p_i == SMALL_MP[tk]) { is_target = true; break; }
        }
        if (is_target) {
            int cpos = atomicAdd(g_sieve.cand_count, 1);
            if (cpos < MAX_SIEVE_CANDS) {
                Candidate c;
                c.slot_idx   = i;
                c.score      = acc;
                c.ll_seed_f  = 4.0f;
                c.r_harmonic = (float)p_i;
                c.phase      = ph;
                c.coherence  = sum_sin;
                c.amp        = residue;
                g_sieve.cand_buffer[cpos] = c;
            }
        }
        acc = 0.0f;
        s = 4u;
    }

    // Write back
    g_sieve.sieve_state[i]   = s;
    g_sieve.r_harmonic[i]    = r_h;
    g_sieve.phase[i]         = ph;
    g_sieve.reward_accum[i]  = acc;
}

// ============================================================================
// Exponent assignment kernel — 8-band mode.
// Slots are divided into N_SMALL_MP=8 equal bands, each centred on a
// confirmed small Mersenne prime exponent that the warp-LL engine can verify
// (all <= 3217 bits, within LIMBS=64 x 61 = 3904-bit capacity).
// Offset is always even so centre (odd prime exponent) +/- even = odd.
// ============================================================================

__global__ void hdgl_sieve_assign_8band_kernel(
    uint32_t *d_assigned_p, int N, const uint32_t *centers, int n_centers)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i >= N) return;

    int band_size = N / n_centers;
    int band      = i / band_size;
    if (band >= n_centers) band = n_centers - 1;

    int local     = i - band * band_size;
    int half      = band_size / 2;
    int offset    = (local - half) * 2;   // even offset -> parity preserved

    int32_t p = (int32_t)centers[band] + offset;
    if (p < 5) p = 5;                    // clamp; smallest safe prime exponent
    if (p > 3904) p = 3904 - (p & 1 ? 0 : 1);  // stay within LIMBS capacity
    d_assigned_p[i] = (uint32_t)p;
}

// ============================================================================
// Seed kernel: reset all LL states to S_0 = 4, phase = i/N * 2π
// ============================================================================

__global__ void hdgl_sieve_seed_kernel(
    uint32_t *d_state, float *d_phase, float *d_r_h, float *d_acc,
    int N)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i >= N) return;
    d_state[i] = 4u;
    d_phase[i] = 6.28318530717958647f * (float)i / (float)N;
    d_r_h[i]   = 10.0f;
    d_acc[i]   = 0.0f;
}

// ============================================================================
// Host-visible data structures and helpers
// ============================================================================

typedef struct {
    // Device pointers
    uint32_t  *d_assigned_p;
    uint32_t  *d_sieve_state;
    float     *d_r_harmonic;
    float     *d_phase;
    float     *d_reward_accum;
    uint32_t  *d_prime_found;
    int       *d_prime_count;
    Candidate *d_cand_buf;     // for warp-LL routing
    int       *d_cand_count;
    CSSState  *d_css;
    float     *d_ll_residue;
    int8_t    *d_ll_verified;
    // Host pinned buffers
    uint32_t  *h_prime_found;
    int       *h_prime_count;
    Candidate *h_cand_buf;
    int       *h_cand_count;
    // Dimensions
    int N;
    int S;
    float dt;
} SieveHostState;

extern "C" {

// ============================================================================
// Allocate and initialise the sieve
// ============================================================================

int hdgl_sieve_alloc(SieveHostState *st, int N, int S, float dt) {
    st->N  = N;
    st->S  = S;
    st->dt = dt;

    cudaMalloc(&st->d_assigned_p,   N * sizeof(uint32_t));
    cudaMalloc(&st->d_sieve_state,  N * sizeof(uint32_t));
    cudaMalloc(&st->d_r_harmonic,   N * sizeof(float));
    cudaMalloc(&st->d_phase,        N * sizeof(float));
    cudaMalloc(&st->d_reward_accum, N * sizeof(float));
    cudaMalloc(&st->d_prime_found,  SIEVE_PRIME_RING * sizeof(uint32_t));
    cudaMalloc(&st->d_prime_count,  sizeof(int));
    cudaMalloc(&st->d_cand_buf,     MAX_SIEVE_CANDS * sizeof(Candidate));
    cudaMalloc(&st->d_cand_count,   sizeof(int));
    cudaMalloc(&st->d_ll_residue,   MAX_SIEVE_CANDS * sizeof(float));
    cudaMalloc(&st->d_ll_verified,  N * sizeof(int8_t));
    hdgl_warp_ll_v33_alloc(&st->d_css);

    cudaMemset(st->d_prime_count, 0, sizeof(int));
    cudaMemset(st->d_prime_found, 0, SIEVE_PRIME_RING * sizeof(uint32_t));
    cudaMemset(st->d_cand_count,  0, sizeof(int));
    cudaMemset(st->d_cand_buf,    0, MAX_SIEVE_CANDS * sizeof(Candidate));
    cudaMemset(st->d_ll_verified, 0, N * sizeof(int8_t));

    cudaMallocHost(&st->h_prime_found,  SIEVE_PRIME_RING * sizeof(uint32_t));
    cudaMallocHost(&st->h_prime_count,  sizeof(int));
    cudaMallocHost(&st->h_cand_buf,     MAX_SIEVE_CANDS * sizeof(Candidate));
    cudaMallocHost(&st->h_cand_count,   sizeof(int));

    // Upload small-prime centres to device constant memory
    uint32_t *d_centers = NULL;
    cudaMalloc(&d_centers, N_SMALL_MP * sizeof(uint32_t));
    cudaMemcpy(d_centers, SMALL_MP_HOST, N_SMALL_MP * sizeof(uint32_t),
               cudaMemcpyHostToDevice);

    // Assign exponents: 8 bands centred on the 8 small confirmed Mersenne primes
    int blk = SIEVE_BLOCK;
    int grd = (N + blk - 1) / blk;
    hdgl_sieve_assign_8band_kernel<<<grd, blk>>>(
        st->d_assigned_p, N, d_centers, N_SMALL_MP);
    cudaFree(d_centers);

    // Seed states
    hdgl_sieve_seed_kernel<<<grd, blk>>>(
        st->d_sieve_state, st->d_phase, st->d_r_harmonic, st->d_reward_accum, N);
    cudaDeviceSynchronize();

    // Upload to __constant__
    SieveDevState dev;
    dev.assigned_p   = st->d_assigned_p;
    dev.sieve_state  = st->d_sieve_state;
    dev.r_harmonic   = st->d_r_harmonic;
    dev.phase        = st->d_phase;
    dev.reward_accum = st->d_reward_accum;
    dev.prime_found  = st->d_prime_found;
    dev.prime_count  = st->d_prime_count;
    dev.cand_buffer  = st->d_cand_buf;
    dev.cand_count   = st->d_cand_count;
    dev.N            = N;
    dev.S            = S;
    dev.dt           = dt;
    cudaMemcpyToSymbol(g_sieve, &dev, sizeof(SieveDevState));

    return 0;
}

// ============================================================================
// Retarget sieve to a custom set of band centers.
// Must be called AFTER hdgl_sieve_alloc().
// Updates both the device __constant__ SMALL_MP (for the is_target guard)
// and re-runs the band assignment kernel with the new centers.
// n_targets must be <= N_SMALL_MP.
// ============================================================================
void hdgl_sieve_retarget(SieveHostState *st, const uint32_t *targets, int n_targets)
{
    if (n_targets > N_SMALL_MP) n_targets = N_SMALL_MP;

    // Pad remainder with last target so kernels stay within bounds
    uint32_t tpad[N_SMALL_MP];
    for (int i = 0; i < N_SMALL_MP; i++)
        tpad[i] = (i < n_targets) ? targets[i] : targets[n_targets - 1];

    // Update device constant (the is_target guard in hdgl_sieve_kernel)
    cudaMemcpyToSymbol(SMALL_MP, tpad, N_SMALL_MP * sizeof(uint32_t));

    // Re-run band assignment
    uint32_t *d_centers = NULL;
    cudaMalloc(&d_centers, n_targets * sizeof(uint32_t));
    cudaMemcpy(d_centers, targets, n_targets * sizeof(uint32_t),
               cudaMemcpyHostToDevice);
    int blk = SIEVE_BLOCK;
    int grd = (st->N + blk - 1) / blk;
    hdgl_sieve_assign_8band_kernel<<<grd, blk>>>(
        st->d_assigned_p, st->N, d_centers, n_targets);
    cudaFree(d_centers);

    // Re-seed phases so slots start fresh around new centers
    hdgl_sieve_seed_kernel<<<grd, blk>>>(
        st->d_sieve_state, st->d_phase, st->d_r_harmonic, st->d_reward_accum, st->N);
    cudaDeviceSynchronize();
}

// ============================================================================
// Free device/pinned memory
// ============================================================================

void hdgl_sieve_free(SieveHostState *st) {
    cudaFree(st->d_assigned_p);
    cudaFree(st->d_sieve_state);
    cudaFree(st->d_r_harmonic);
    cudaFree(st->d_phase);
    cudaFree(st->d_reward_accum);
    cudaFree(st->d_prime_found);
    cudaFree(st->d_prime_count);
    cudaFree(st->d_cand_buf);
    cudaFree(st->d_cand_count);
    cudaFree(st->d_ll_residue);
    cudaFree(st->d_ll_verified);
    if (st->d_css) hdgl_warp_ll_v33_free(st->d_css);
    cudaFreeHost(st->h_prime_found);
    cudaFreeHost(st->h_prime_count);
    cudaFreeHost(st->h_cand_buf);
    cudaFreeHost(st->h_cand_count);
    memset(st, 0, sizeof(SieveHostState));
}

// ============================================================================
// Run one sieve step (async)
// ============================================================================

void hdgl_sieve_step(const SieveHostState *st, cudaStream_t stream) {
    int blk = SIEVE_BLOCK;
    int grd = (st->N + blk - 1) / blk;
    hdgl_sieve_kernel<<<grd, blk, 0, stream>>>();
}

// ============================================================================
// Harvest ring buffer: copy prime_found and prime_count to host
// Returns number of new found exponents
// ============================================================================

int hdgl_sieve_harvest(SieveHostState *st, uint32_t *out, int max_out,
                       cudaStream_t stream) {
    cudaMemcpyAsync(st->h_prime_count, st->d_prime_count,
                    sizeof(int), cudaMemcpyDeviceToHost, stream);
    cudaMemcpyAsync(st->h_prime_found, st->d_prime_found,
                    SIEVE_PRIME_RING * sizeof(uint32_t),
                    cudaMemcpyDeviceToHost, stream);
    cudaStreamSynchronize(stream);

    int count = *(st->h_prime_count);
    if (count > SIEVE_PRIME_RING) count = SIEVE_PRIME_RING;

    int n = (count < max_out) ? count : max_out;

    // Report proxy hits for known primes — at most one banner per harvest
    int banner_printed = 0;
    for (int i = 0; i < n && !banner_printed; i++) {
        out[i] = st->h_prime_found[i];
        for (int k = 0; k < N_SMALL_MP; k++) {
            if (out[i] == SMALL_MP_HOST[k]) {
                printf("[proxy] Mersenne exponent p=%u hit — routing to warp-LL\n",
                       out[i]);
                fflush(stdout);
                banner_printed = 1;
                break;
            }
        }
    }
    // Fill remaining entries without banner check
    for (int i = (banner_printed ? 1 : 0); i < n; i++) {
        out[i] = st->h_prime_found[i];
    }

    // Reset proxy ring
    cudaMemsetAsync(st->d_prime_count, 0, sizeof(int), stream);
    cudaMemsetAsync(st->d_prime_found, 0,
                    SIEVE_PRIME_RING * sizeof(uint32_t), stream);
    return n;
}

// ============================================================================
// Harvest and warp-LL verify: routes Candidates through the real LL engine.
// Returns number of warp-LL-confirmed primes; writes their exponents to out[].
// This is the mathematically meaningful path \u2014 proxy ring is just the pre-filter.
//
// For p \u2264 3904 (all targets in the 8-band sieve), the warp-LL NTT runs the
// correct number of active limbs.  Residue \u2248 0 after p-2 iterations
// indicates 2^p-1 is likely prime (exact for p=61 since M61 is the NTT modulus).
// ============================================================================

int hdgl_sieve_harvest_verified(SieveHostState *st, uint32_t *out, int max_out,
                                 int ll_iters, float ll_eps, cudaStream_t stream) {
    // Copy candidate count
    cudaMemcpyAsync(st->h_cand_count, st->d_cand_count,
                    sizeof(int), cudaMemcpyDeviceToHost, stream);
    cudaStreamSynchronize(stream);

    int n_cands = *(st->h_cand_count);
    if (n_cands <= 0) return 0;
    if (n_cands > MAX_SIEVE_CANDS) n_cands = MAX_SIEVE_CANDS;

    // Clear verified flags (indexed by cand_idx) before this batch runs
    cudaMemsetAsync(st->d_ll_verified, 0, MAX_SIEVE_CANDS * sizeof(int8_t), stream);

    // Seed CSS states and run real warp-LL
    hdgl_warp_ll_v33_seed(st->d_css, n_cands, stream);
    hdgl_warp_ll_v33_launch(st->d_cand_buf, n_cands,
                             0,        // p_bits=0 \u2192 kernel uses per-cand r_harmonic
                             ll_iters,
                             st->d_css,
                             st->d_ll_residue,
                             st->d_ll_verified,
                             ll_eps, stream);

    // Copy back residues and candidates
    cudaMemcpyAsync(st->h_cand_buf, st->d_cand_buf,
                    n_cands * sizeof(Candidate),
                    cudaMemcpyDeviceToHost, stream);

    float h_residue[MAX_SIEVE_CANDS];
    int8_t h_verified[MAX_SIEVE_CANDS];
    cudaMemcpyAsync(h_residue, st->d_ll_residue,
                    n_cands * sizeof(float), cudaMemcpyDeviceToHost, stream);
    cudaMemcpyAsync(h_verified, st->d_ll_verified,
                    n_cands * sizeof(int8_t), cudaMemcpyDeviceToHost, stream);
    cudaStreamSynchronize(stream);

    // Reset candidate counter for next cycle
    cudaMemsetAsync(st->d_cand_count, 0, sizeof(int), stream);

    int n_confirmed = 0;
    for (int i = 0; i < n_cands && n_confirmed < max_out; i++) {
        uint32_t p = (uint32_t)st->h_cand_buf[i].r_harmonic;
        float    res = h_residue[i];

        if (h_verified[i] == 1) {
            out[n_confirmed++] = p;
            printf("\n");
            printf("============================================================\n");
            printf("  WARP-LL CONFIRMED: 2^%u - 1  residue=%.2e\n", p, (double)res);
            // Check against known small Mersenne primes
            for (int k = 0; k < N_SMALL_MP; k++) {
                if (p == SMALL_MP_HOST[k]) {
                    printf("  *** KNOWN PRIME HIT: M%u is confirmed! ***\n", p);
                }
            }
            printf("============================================================\n\n");
            fflush(stdout);
        } else {
            // Still print candidates with low residue for monitoring
            if (res < 0.05f) {
                printf("  [sieve] p=%u  LL-residue=%.4e  verified=%d\n",
                       p, (double)res, (int)h_verified[i]);
            }
        }
    }
    return n_confirmed;
}

} // extern "C"
