// ============================================================================
// HDGL Benchmark & Test Harness v33
// ============================================================================
//
// Tests four subsystems independently, then runs an integrated benchmark:
//
//   TEST 1 — Critic (CPU):
//     • Init → forward pass on known input → verify output range
//     • 100 observe+update cycles → verify weights change
//     • TD target computation
//
//   TEST 2 — NTT arithmetic (GPU):
//     • Load a known 64-limb number, NTT-square it, verify result vs direct
//
//   TEST 3 — Sieve kernel (GPU):
//     • Allocate 4096 slots, run 100 sieve steps, verify ring buffer is
//       populated (all candidates come back in [0,1] residue range)
//
//   TEST 4 — Field kernel (GPU):
//     • Allocate N=65536 slots, run 10 field steps, verify no NaN/Inf
//
//   BENCHMARK — Field kernel throughput:
//     • N = 256K, 512K, 1M slots
//     • 500 field steps each
//     • Report ns/step and GSlots/s
//
// Compile:
//   nvcc -O3 -arch=sm_75 hdgl_analog_v33.cu hdgl_warp_ll_v33.cu \
//        hdgl_sieve_v34.cu hdgl_critic_v33.c hdgl_bench_v33.cu -o hdgl_bench -lm
// ============================================================================

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

// Windows-compatible high-resolution timer
#ifdef _WIN32
#define WIN32_LEAN_AND_MEAN
#include <windows.h>
static LARGE_INTEGER g_qpf_bench;
static void bench_clock_init(void) { QueryPerformanceFrequency(&g_qpf_bench); }
static double now_sec(void) {
    LARGE_INTEGER t; QueryPerformanceCounter(&t);
    return (double)t.QuadPart / (double)g_qpf_bench.QuadPart;
}
#else
#include <time.h>
static void bench_clock_init(void) {}
static double now_sec(void) {
    struct timespec t; clock_gettime(CLOCK_MONOTONIC, &t);
    return t.tv_sec + t.tv_nsec * 1e-9;
}
#endif

// ============================================================================
// Forward declarations
// ============================================================================

typedef struct {
    int   slot_idx;
    float score;
    float ll_seed_f;
    float r_harmonic;
    float phase;
    float coherence;   // |sum_sin| at emission — feature for critic
    float amp;         // local amplitude at emission — feature for critic
} Candidate;

typedef struct {
    float    *A_re, *A_im, *phase, *phase_vel, *r_harmonic;
    uint32_t *ll_state;
    float    *reward_accum;
    float    *w_cos, *w_sin, *w_sigma;
    int8_t   *ll_verified;
    Candidate *candidates;
    int      *cand_count;
    int       N, S;
    float     omega, dt, w_amp_self, w_amp_neigh;
} DevSoA;

typedef struct { uint64_t S[64]; uint64_t C[64]; } CSSState;

typedef struct {
    uint32_t  *d_assigned_p, *d_sieve_state;
    float     *d_r_harmonic, *d_phase, *d_reward_accum;
    uint32_t  *d_prime_found;
    int       *d_prime_count;
    Candidate *d_cand_buf;
    int       *d_cand_count;
    CSSState  *d_css;
    float     *d_ll_residue;
    int8_t    *d_ll_verified;
    uint32_t  *h_prime_found;
    int       *h_prime_count;
    Candidate *h_cand_buf;
    int       *h_cand_count;
    int N, S; float dt;
} SieveHostState;

#define SPECTRAL_N     4
#define CRITIC_W_TOTAL 57
#define MAX_CAND_BUF   256

// All cross-TU functions are compiled as C or exported with extern "C",
// so declare them with extern "C" here to match linkage.
// (critic_* functions come from hdgl_critic_v33.h which already has the guard)
extern "C" {
    // hdgl_analog_v33.cu
    void hdgl_v33_upload_soa(const DevSoA *h);
    void hdgl_v33_upload_critic(const float *w, int n);
    void hdgl_v33_field_step(int N, int block, cudaStream_t s);
    void hdgl_v33_weight_sync(cudaStream_t s);
    void hdgl_v33_read_global_weights(float oc[], float os[], float osg[]);

    // hdgl_warp_ll_v33.cu
    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 *dc, int nc, int pb, int it,
                                  CSSState *dc2, float *dr, int8_t *dv,
                                  float eps, cudaStream_t s);

    // hdgl_sieve_v34.cu
    int  hdgl_sieve_alloc(SieveHostState *st, int N, int S, float dt);
    void hdgl_sieve_free(SieveHostState *st);
    void hdgl_sieve_step(const SieveHostState *st, cudaStream_t s);
    int  hdgl_sieve_harvest(SieveHostState *st, uint32_t *out, int max, cudaStream_t s);
}  // extern "C"

// ============================================================================
// Helpers
// ============================================================================

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

// now_sec() defined conditionally near the top of this file

static int g_pass = 0, g_fail = 0;

static void check(const char *label, int cond) {
    if (cond) { printf("  [PASS] %s\n", label); g_pass++; }
    else       { printf("  [FAIL] %s\n", label); g_fail++; }
}

// ============================================================================
// TEST 1 — Critic (pure CPU)
// ============================================================================

static void test_critic(void) {
    printf("\n=== TEST 1: Critic (CPU) ===\n");
    critic_init();

    // Known input: all zeros → output should be close to b2 (small Xavier bias)
    float s_zero[CRITIC_IN] = {0};
    float out0 = critic_forward(s_zero);
    check("forward(zeros) is finite", isfinite(out0));

    // Known input: all ones
    float s_ones[CRITIC_IN] = {1.0f, 1.0f, 1.0f, 1.0f, 1.0f};
    float out1 = critic_forward(s_ones);
    check("forward(ones) is finite", isfinite(out1));

    // 200 observe + update cycles: weights should change
    float w_before[CRITIC_W_TOTAL], w_after[CRITIC_W_TOTAL];
    critic_pack_weights(w_before);

    for (int i = 0; i < 200; i++) {
        float feat[CRITIC_IN] = {
            (float)i / 200.0f, 0.5f, 0.8f, 0.1f, 0.05f
        };
        float target = 1.0f / (1.0f + (float)i * 0.01f);
        critic_observe(feat, target);
    }
    critic_pack_weights(w_after);

    float max_diff = 0.0f;
    for (int i = 0; i < CRITIC_W_TOTAL; i++) {
        float d = fabsf(w_after[i] - w_before[i]);
        if (d > max_diff) max_diff = d;
    }
    check("weights changed after 200 updates", max_diff > 1e-6f);
    printf("  max_weight_delta = %.4e\n", (double)max_diff);

    // TD target: r + gamma * V(s')
    float s_next[CRITIC_IN] = {0.1f, 0.4f, 0.9f, 0.2f, 0.1f};
    float td = critic_td_target(0.5f, s_next);
    check("td_target is finite", isfinite(td));
    printf("  td_target(r=0.5, s') = %.4f\n", (double)td);

    // Verify pack/unpack round-trip
    float packed[CRITIC_W_TOTAL];
    critic_pack_weights(packed);
    check("critic_weight_count() == CRITIC_W_TOTAL",
          critic_weight_count() == CRITIC_W_TOTAL);

    critic_print_stats();
}

// ============================================================================
// TEST 2 — NTT + Warp LL (GPU)
// ============================================================================

static void test_warp_ll(void) {
    printf("\n=== TEST 2: Warp LL v33 (GPU) ===\n");

    CSSState *d_css = NULL;
    float    *d_res = NULL;
    int8_t   *d_ver = NULL;
    Candidate *d_cands = NULL;

    hdgl_warp_ll_v33_alloc(&d_css);
    CUDA_CHECK(cudaMalloc(&d_res,   8 * sizeof(float)));
    CUDA_CHECK(cudaMalloc(&d_ver,   8 * sizeof(int8_t)));
    CUDA_CHECK(cudaMalloc(&d_cands, 8 * sizeof(Candidate)));

    // Set up 4 candidates: first 2 target M61 (exact direct path),
    // last 2 target M127 (multi-limb NTT path — approximate but should be finite)
    Candidate h_cands[4];
    for (int i = 0; i < 4; i++) {
        h_cands[i].slot_idx   = i;
        h_cands[i].score      = 1.0f;
        h_cands[i].ll_seed_f  = 4.0f;
        h_cands[i].r_harmonic = (i < 2) ? 61.0f : 127.0f;
        h_cands[i].phase      = 0.0f;
    }
    CUDA_CHECK(cudaMemcpy(d_cands, h_cands, 4 * sizeof(Candidate), cudaMemcpyHostToDevice));

    cudaStream_t s;
    CUDA_CHECK(cudaStreamCreate(&s));

    // p_bits=0 → use per-candidate r_harmonic.
    // M61 needs exactly 59 LL iterations (p-2) to reach residue=0.
    hdgl_warp_ll_v33_seed(d_css, 4, s);
    hdgl_warp_ll_v33_launch(d_cands, 4, /*p_bits*/0, /*iters*/59, d_css, d_res, d_ver, 1e-6f, s);
    CUDA_CHECK(cudaStreamSynchronize(s));

    float h_res[4];
    int8_t h_ver[4];
    CUDA_CHECK(cudaMemcpy(h_res, d_res, 4 * sizeof(float), cudaMemcpyDeviceToHost));
    CUDA_CHECK(cudaMemcpy(h_ver, d_ver, 4 * sizeof(int8_t), cudaMemcpyDeviceToHost));

    int all_finite = 1, all_normalised = 1, m61_confirmed = 1;
    for (int i = 0; i < 4; i++) {
        printf("  cand[%d] (p=%d): residue=%.4e  verified=%d\n",
               i, (i < 2) ? 61 : 127, (double)h_res[i], (int)h_ver[i]);
        if (!isfinite(h_res[i])) all_finite = 0;
        if (h_res[i] < 0.0f || h_res[i] > 1.0f) all_normalised = 0;
        if (i < 2 && h_ver[i] != 1) m61_confirmed = 0;
    }
    check("all residues are finite",    all_finite);
    check("all residues in [0,1]",      all_normalised);
    check("M61 verified by direct LL",  m61_confirmed);

    cudaStreamDestroy(s);
    hdgl_warp_ll_v33_free(d_css);
    cudaFree(d_res);
    cudaFree(d_ver);
    cudaFree(d_cands);
}

// ============================================================================
// TEST 3 — Sieve kernel (GPU)
// ============================================================================

static void test_sieve(void) {
    printf("\n=== TEST 3: Sieve v34 (GPU) ===\n");

    SieveHostState sieve;
    int N = 4096, S = 64;
    hdgl_sieve_alloc(&sieve, N, S, 0.001f);

    cudaStream_t s;
    CUDA_CHECK(cudaStreamCreate(&s));

    // Run 200 sieve steps
    for (int i = 0; i < 200; i++) hdgl_sieve_step(&sieve, s);
    CUDA_CHECK(cudaStreamSynchronize(s));

    // Harvest — just verify no crash and count is reasonable
    uint32_t found[256];
    int n = hdgl_sieve_harvest(&sieve, found, 256, s);
    printf("  harvested %d candidate exponents after 200 steps\n", n);
    check("harvest returns non-negative count", n >= 0);
    check("harvest count within ring capacity",  n <= 256);

    if (n > 0) {
        // Exponents now span 8 small-prime bands: [61, 3217] ± band_offset
        check("first exponent is in verifiable range",
              found[0] >= 5u && found[0] <= 4096u);
        printf("  first exponent: p = %u\n", found[0]);
    }

    // Verify phase array has no NaN after 200 steps
    float *h_phase = (float *)malloc(N * sizeof(float));
    CUDA_CHECK(cudaMemcpy(h_phase, sieve.d_phase, N * sizeof(float), cudaMemcpyDeviceToHost));
    int nan_count = 0;
    for (int i = 0; i < N; i++) if (!isfinite(h_phase[i])) nan_count++;
    check("no NaN/Inf in phase after 200 steps", nan_count == 0);
    printf("  nan_count_in_phase = %d / %d\n", nan_count, N);
    free(h_phase);

    cudaStreamDestroy(s);
    hdgl_sieve_free(&sieve);
}

// ============================================================================
// TEST 4 — Field kernel (GPU)
// ============================================================================

static void test_field(void) {
    printf("\n=== TEST 4: Field kernel v33 (GPU) ===\n");

    const int N = 65536, S = 256;
    float *d_Ar, *d_Ai, *d_ph, *d_phv, *d_rh, *d_acc;
    float *d_wc, *d_ws, *d_wsg;
    uint32_t *d_ll;
    int8_t *d_ver;
    Candidate *d_cands;
    int *d_cnt;

    CUDA_CHECK(cudaMalloc(&d_Ar,    N * sizeof(float)));
    CUDA_CHECK(cudaMalloc(&d_Ai,    N * sizeof(float)));
    CUDA_CHECK(cudaMalloc(&d_ph,    N * sizeof(float)));
    CUDA_CHECK(cudaMalloc(&d_phv,   N * sizeof(float)));
    CUDA_CHECK(cudaMalloc(&d_rh,    N * sizeof(float)));
    CUDA_CHECK(cudaMalloc(&d_acc,   N * sizeof(float)));
    CUDA_CHECK(cudaMalloc(&d_wc,    N * SPECTRAL_N * sizeof(float)));
    CUDA_CHECK(cudaMalloc(&d_ws,    N * SPECTRAL_N * sizeof(float)));
    CUDA_CHECK(cudaMalloc(&d_wsg,   N * SPECTRAL_N * sizeof(float)));
    CUDA_CHECK(cudaMalloc(&d_ll,    N * sizeof(uint32_t)));
    CUDA_CHECK(cudaMalloc(&d_ver,   N * sizeof(int8_t)));
    CUDA_CHECK(cudaMalloc(&d_cands, MAX_CAND_BUF * sizeof(Candidate)));
    CUDA_CHECK(cudaMalloc(&d_cnt,   sizeof(int)));

    // Initialise with simple values
    float *tmp = (float *)malloc(N * SPECTRAL_N * sizeof(float));
    float sigma_init[4] = {3.14159f, 1.5708f, 0.7854f, 0.3927f};

    for (int i = 0; i < N; i++) ((float *)tmp)[i] = 1.0f;
    CUDA_CHECK(cudaMemcpy(d_Ar, tmp, N * sizeof(float), cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemset(d_Ai, 0, N * sizeof(float)));

    for (int i = 0; i < N; i++)
        tmp[i] = 6.28318530f * (float)i / (float)N;
    CUDA_CHECK(cudaMemcpy(d_ph, tmp, N * sizeof(float), cudaMemcpyHostToDevice));

    for (int i = 0; i < N; i++) tmp[i] = 0.314159f;
    CUDA_CHECK(cudaMemcpy(d_phv, tmp, N * sizeof(float), cudaMemcpyHostToDevice));

    // r_harmonic spread over [2, 100]: covers φ-resonance nodes that give S<0.25
    for (int i = 0; i < N; i++) tmp[i] = 2.0f + 98.0f * ((float)rand() / (float)RAND_MAX);
    CUDA_CHECK(cudaMemcpy(d_rh, tmp, N * sizeof(float), cudaMemcpyHostToDevice));

    CUDA_CHECK(cudaMemset(d_acc, 0, N * sizeof(float)));
    CUDA_CHECK(cudaMemset(d_ver, 0, N * sizeof(int8_t)));
    CUDA_CHECK(cudaMemset(d_cnt, 0, sizeof(int)));

    // LL state: non-zero seeds
    uint32_t *ll_tmp = (uint32_t *)malloc(N * sizeof(uint32_t));
    for (int i = 0; i < N; i++) ll_tmp[i] = (uint32_t)(i * 2654435761u + 7u) | 1u;
    CUDA_CHECK(cudaMemcpy(d_ll, ll_tmp, N * sizeof(uint32_t), cudaMemcpyHostToDevice));
    free(ll_tmp);

    for (int i = 0; i < N; i++)
        for (int k = 0; k < SPECTRAL_N; k++) {
            tmp[i * SPECTRAL_N + k] = 0.05f * ((float)rand() / RAND_MAX - 0.5f);
        }
    CUDA_CHECK(cudaMemcpy(d_wc, tmp, N * SPECTRAL_N * sizeof(float), cudaMemcpyHostToDevice));

    for (int i = 0; i < N; i++)
        for (int k = 0; k < SPECTRAL_N; k++)
            tmp[i * SPECTRAL_N + k] = 0.05f * ((float)rand() / RAND_MAX - 0.5f);
    CUDA_CHECK(cudaMemcpy(d_ws, tmp, N * SPECTRAL_N * sizeof(float), cudaMemcpyHostToDevice));

    for (int i = 0; i < N; i++)
        for (int k = 0; k < SPECTRAL_N; k++)
            tmp[i * SPECTRAL_N + k] = sigma_init[k];
    CUDA_CHECK(cudaMemcpy(d_wsg, tmp, N * SPECTRAL_N * sizeof(float), cudaMemcpyHostToDevice));

    free(tmp);

    // Upload SoA
    float critic_w[CRITIC_W_TOTAL] = {0};
    critic_init();
    critic_pack_weights(critic_w);

    DevSoA soa = {
        d_Ar, d_Ai, d_ph, d_phv, d_rh, d_ll, d_acc,
        d_wc, d_ws, d_wsg, d_ver, d_cands, d_cnt,
        N, S, 0.314159f, 0.001f, 1.0f, 0.5f
    };
    hdgl_v33_upload_soa(&soa);
    hdgl_v33_upload_critic(critic_w, CRITIC_W_TOTAL);

    cudaStream_t stream;
    CUDA_CHECK(cudaStreamCreate(&stream));

    // Run 20 field steps — resonant slots need ~10 steps to accumulate acc > CAND_ACCUM_THRESH
    for (int i = 0; i < 20; i++) hdgl_v33_field_step(N, 256, stream);
    CUDA_CHECK(cudaStreamSynchronize(stream));

    // Check for NaN/Inf in A_re
    float *h_Ar = (float *)malloc(N * sizeof(float));
    CUDA_CHECK(cudaMemcpy(h_Ar, d_Ar, N * sizeof(float), cudaMemcpyDeviceToHost));
    int nan_cnt = 0, inf_cnt = 0;
    float sum_abs = 0.0f;
    for (int i = 0; i < N; i++) {
        if (isnan(h_Ar[i]))  nan_cnt++;
        else if (isinf(h_Ar[i])) inf_cnt++;
        else sum_abs += fabsf(h_Ar[i]);
    }
    check("no NaN in A_re after 20 field steps",  nan_cnt == 0);
    check("no Inf in A_re after 20 field steps",  inf_cnt == 0);
    printf("  nan=%d  inf=%d  mean|A_re|=%.4f\n",
           nan_cnt, inf_cnt, (double)(sum_abs / N));
    free(h_Ar);

    // Check candidate count is reasonable
    int h_cnt = -1;
    CUDA_CHECK(cudaMemcpy(&h_cnt, d_cnt, sizeof(int), cudaMemcpyDeviceToHost));
    check("candidate count >= 0",   h_cnt >= 0);
    check("candidate count > 0",    h_cnt > 0);   // gate must fire with spread r_harmonic
    check("candidate count <= N",   h_cnt <= N);   // sanity: can't exceed total slots
    printf("  candidates promoted = %d  (buffer capacity = %d)\n", h_cnt, MAX_CAND_BUF);

    cudaStreamDestroy(stream);
    cudaFree(d_Ar);  cudaFree(d_Ai); cudaFree(d_ph);  cudaFree(d_phv);
    cudaFree(d_rh);  cudaFree(d_acc);cudaFree(d_wc);  cudaFree(d_ws);
    cudaFree(d_wsg); cudaFree(d_ll); cudaFree(d_ver); cudaFree(d_cands);
    cudaFree(d_cnt);
}

// ============================================================================
// BENCHMARK — Field kernel throughput
// ============================================================================

static void benchmark_field(void) {
    printf("\n=== BENCHMARK: Field kernel throughput ===\n");
    printf("  %-12s  %-10s  %-14s  %-12s\n",
           "N", "steps", "wall_ms", "GSlots/s");
    printf("  ------------------------------------------------------------\n");

    static const int N_vals[] = {262144, 524288, 1048576};
    const int STEPS = 500;
    const int S = 1024;
    const int BLK = 256;

    for (int ni = 0; ni < 3; ni++) {
        int N = N_vals[ni];

        // Alloc device arrays
        float *d_Ar, *d_Ai, *d_ph, *d_phv, *d_rh, *d_acc;
        float *d_wc, *d_ws, *d_wsg;
        uint32_t *d_ll;
        int8_t *d_ver;
        Candidate *d_cands;
        int *d_cnt;

        CUDA_CHECK(cudaMalloc(&d_Ar,  N * sizeof(float)));
        CUDA_CHECK(cudaMalloc(&d_Ai,  N * sizeof(float)));
        CUDA_CHECK(cudaMalloc(&d_ph,  N * sizeof(float)));
        CUDA_CHECK(cudaMalloc(&d_phv, N * sizeof(float)));
        CUDA_CHECK(cudaMalloc(&d_rh,  N * sizeof(float)));
        CUDA_CHECK(cudaMalloc(&d_acc, N * sizeof(float)));
        CUDA_CHECK(cudaMalloc(&d_wc,  N * SPECTRAL_N * sizeof(float)));
        CUDA_CHECK(cudaMalloc(&d_ws,  N * SPECTRAL_N * sizeof(float)));
        CUDA_CHECK(cudaMalloc(&d_wsg, N * SPECTRAL_N * sizeof(float)));
        CUDA_CHECK(cudaMalloc(&d_ll,  N * sizeof(uint32_t)));
        CUDA_CHECK(cudaMalloc(&d_ver, N * sizeof(int8_t)));
        CUDA_CHECK(cudaMalloc(&d_cands, MAX_CAND_BUF * sizeof(Candidate)));
        CUDA_CHECK(cudaMalloc(&d_cnt,   sizeof(int)));

        // Simple init via memset
        CUDA_CHECK(cudaMemset(d_Ar,  0x3f, N * sizeof(float)));   // ≈ 1.0 in fp32
        CUDA_CHECK(cudaMemset(d_Ai,  0, N * sizeof(float)));
        CUDA_CHECK(cudaMemset(d_ph,  0, N * sizeof(float)));
        CUDA_CHECK(cudaMemset(d_phv, 0x3e, N * sizeof(float)));   // ≈ 0.5
        CUDA_CHECK(cudaMemset(d_rh,  0x41, N * sizeof(float)));   // ≈ some fp val
        CUDA_CHECK(cudaMemset(d_acc, 0, N * sizeof(float)));
        CUDA_CHECK(cudaMemset(d_wc,  0, N * SPECTRAL_N * sizeof(float)));
        CUDA_CHECK(cudaMemset(d_ws,  0, N * SPECTRAL_N * sizeof(float)));
        CUDA_CHECK(cudaMemset(d_ll,  0x01, N * sizeof(uint32_t))); // non-zero
        CUDA_CHECK(cudaMemset(d_ver, 0, N * sizeof(int8_t)));
        CUDA_CHECK(cudaMemset(d_cnt, 0, sizeof(int)));

        // Init w_sigma to float constants via pattern
        // 0x3f 0xc9 0x0f 0xda = pi in IEEE 754 little-endian (x86: 0xDA0F49C0)
        // Use a simple loop via float array
        {
            float *tmp = (float *)malloc(N * SPECTRAL_N * sizeof(float));
            float sigma[4] = {3.14159f, 1.5708f, 0.7854f, 0.3927f};
            for (int i = 0; i < N; i++)
                for (int k = 0; k < 4; k++) tmp[i*4+k] = sigma[k];
            CUDA_CHECK(cudaMemcpy(d_wsg, tmp, N * SPECTRAL_N * sizeof(float), cudaMemcpyHostToDevice));
            free(tmp);
        }

        DevSoA soa = {
            d_Ar, d_Ai, d_ph, d_phv, d_rh, d_ll, d_acc,
            d_wc, d_ws, d_wsg, d_ver, d_cands, d_cnt,
            N, S, 0.314159f, 0.001f, 1.0f, 0.5f
        };
        hdgl_v33_upload_soa(&soa);

        float cw[CRITIC_W_TOTAL] = {0};
        hdgl_v33_upload_critic(cw, CRITIC_W_TOTAL);

        cudaStream_t stream;
        CUDA_CHECK(cudaStreamCreate(&stream));

        // Warmup
        for (int i = 0; i < 5; i++) hdgl_v33_field_step(N, BLK, stream);
        CUDA_CHECK(cudaStreamSynchronize(stream));

        // Timed run
        double t0 = now_sec();
        for (int i = 0; i < STEPS; i++) hdgl_v33_field_step(N, BLK, stream);
        CUDA_CHECK(cudaStreamSynchronize(stream));
        double t1 = now_sec();

        double elapsed_ms = (t1 - t0) * 1000.0;
        double gslots = (double)STEPS * (double)N / (t1 - t0) / 1e9;

        printf("  %-12d  %-10d  %-14.2f  %-12.3f\n",
               N, STEPS, elapsed_ms, gslots);

        cudaStreamDestroy(stream);
        cudaFree(d_Ar);  cudaFree(d_Ai); cudaFree(d_ph);  cudaFree(d_phv);
        cudaFree(d_rh);  cudaFree(d_acc);cudaFree(d_wc);  cudaFree(d_ws);
        cudaFree(d_wsg); cudaFree(d_ll); cudaFree(d_ver); cudaFree(d_cands);
        cudaFree(d_cnt);
    }
}

// ============================================================================
// main
// ============================================================================

int main(void) {
    printf("HDGL v33 Test & Benchmark Harness\n");
    printf("CUDA device: ");
    cudaDeviceProp prop;
    CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
    printf("%s  cc=%d.%d  smem=%zu KB  sharedPerBlock=%zu KB\n",
           prop.name, prop.major, prop.minor,
           prop.totalGlobalMem / (1024*1024),
           prop.sharedMemPerBlock / 1024);

    bench_clock_init();
    srand(42);

    test_critic();
    test_warp_ll();
    test_sieve();
    test_field();
    benchmark_field();

    printf("\n========================================\n");
    printf("RESULTS: %d passed, %d failed\n", g_pass, g_fail);
    printf("========================================\n");
    return (g_fail > 0) ? 1 : 0;
}
