// ============================================================================
// HDGL Large Prime Search Host
// ============================================================================
//
// Orchestrates a continuous resonance-guided search for Mersenne prime
// candidates ABOVE the current largest known prime (M136279841).
//
// Architecture:
//   - Retargets the 8 sieve bands to regions above M136279841
//   - Runs the full GRA field with BigG power-law coupling
//   - Harvests resonance-passing candidates from the proxy ring
//   - Reports candidates sorted by resonance score
//   - NO warp-LL verification (LL for p~136M needs GIMPS-scale compute)
//   - Critic trains on resonance signal alone
//
// Band centers (8 × 2048 slots, each band ±1024 around its center):
//   {136280003, 136310003, 136340003, 136370003,
//    136400003, 136430003, 136460003, 136490003}
//   Coverage: ~136,279,000 – ~136,491,000 (with ±1024 margins)
//
// Compile:
//   nvcc -O3 -arch=sm_75 hdgl_analog_v33.cu hdgl_warp_ll_v33.cu ^
//        hdgl_sieve_v34.cu hdgl_critic_v33.c hdgl_search_host.cu ^
//        -o hdgl_search
//
// Run:
//   .\hdgl_search.exe [cycles]    (default: 10000)
// ============================================================================

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

#ifdef _WIN32
#define WIN32_LEAN_AND_MEAN
#include <windows.h>
typedef struct { LARGE_INTEGER v; } sh_timespec_t;
static LARGE_INTEGER g_qpf;
static void sh_clock_init(void) { QueryPerformanceFrequency(&g_qpf); }
static void sh_gettime(sh_timespec_t *t) { QueryPerformanceCounter(&t->v); }
static double sh_elapsed(const sh_timespec_t *a, const sh_timespec_t *b) {
    return (double)(b->v.QuadPart - a->v.QuadPart) / (double)g_qpf.QuadPart;
}
#else
#include <time.h>
typedef struct timespec sh_timespec_t;
static void sh_clock_init(void) {}
static void sh_gettime(sh_timespec_t *t) { clock_gettime(CLOCK_MONOTONIC, t); }
static double sh_elapsed(const sh_timespec_t *a, const sh_timespec_t *b) {
    return (b->tv_sec - a->tv_sec) + (b->tv_nsec - a->tv_nsec) * 1e-9;
}
#endif

#include "hdgl_critic_v33.h"

// ============================================================================
// Shared type definitions (must match all TUs exactly)
// ============================================================================

typedef struct {
    int   slot_idx;
    float score;
    float ll_seed_f;
    float r_harmonic;
    float phase;
    float coherence;
    float amp;
} Candidate;

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 CRITIC_W_TOTAL  57
#define MAX_SIEVE_CANDS 256

extern "C" {
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);
void hdgl_sieve_retarget(SieveHostState *st, const uint32_t *targets, int n_targets);
}

// ============================================================================
// Configuration
// ============================================================================

#define SIEVE_N          16384
#define SIEVE_S          128
#define SIEVE_DT         0.001f
#define STEPS_PER_CYCLE  32
#define DEFAULT_CYCLES   10000
#define REPORT_EVERY     200

// Large-search band centers — 8 bands × 2048 slots, ~30000 apart for
// ~240k exponent coverage starting just above M136279841.
#define N_LARGE_TARGETS  8
static const uint32_t LARGE_TARGETS[N_LARGE_TARGETS] = {
    136280003u, 136310003u, 136340003u, 136370003u,
    136400003u, 136430003u, 136460003u, 136490003u
};

// Track top-scoring candidates across the whole run
#define TOP_CAND_MAX 1024
typedef struct { uint32_t p; float score; int hits; } TopCand;

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

static volatile int g_running = 1;
static void sig_handler(int s) { (void)s; g_running = 0; }

static TopCand g_top[TOP_CAND_MAX];
static int     g_n_top = 0;

// Record or update a candidate hit
static void record_candidate(uint32_t p, float score)
{
    for (int i = 0; i < g_n_top; i++) {
        if (g_top[i].p == p) {
            g_top[i].hits++;
            if (score < g_top[i].score) g_top[i].score = score;
            return;
        }
    }
    if (g_n_top < TOP_CAND_MAX) {
        g_top[g_n_top].p     = p;
        g_top[g_n_top].score = score;
        g_top[g_n_top].hits  = 1;
        g_n_top++;
    }
}

static int cmp_top(const void *a, const void *b)
{
    // Sort by hits desc, then score asc
    const TopCand *ta = (const TopCand*)a;
    const TopCand *tb = (const TopCand*)b;
    if (tb->hits != ta->hits) return tb->hits - ta->hits;
    return (ta->score > tb->score) - (ta->score < tb->score);
}

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

int main(int argc, char **argv)
{
    int total_cycles = (argc > 1) ? atoi(argv[1]) : DEFAULT_CYCLES;
    if (total_cycles <= 0) total_cycles = DEFAULT_CYCLES;

    signal(SIGINT,  sig_handler);
    signal(SIGTERM, sig_handler);
    sh_clock_init();

    printf("====================================================================\n");
    printf("  HDGL Large Prime Search — resonance-guided above M136279841\n");
    printf("  Targets: 8 bands × 2048 slots above exponent 136,279,841\n");
    printf("  Coverage: ~136,279,000 – ~136,491,000\n");
    printf("  Mode: φ-resonance scoring only (LL verification requires GIMPS)\n");
    printf("  Cycles: %d   Press Ctrl-C to stop early.\n", total_cycles);
    printf("====================================================================\n\n");

    // Init critic
    critic_init();
    critic_load_weights("hdgl_critic_search.bin");

    // Allocate sieve with default targets, then retarget
    SieveHostState sieve;
    memset(&sieve, 0, sizeof(sieve));
    if (hdgl_sieve_alloc(&sieve, SIEVE_N, SIEVE_S, SIEVE_DT) != 0) {
        fprintf(stderr, "Failed to allocate sieve\n");
        return 1;
    }

    // Point all 8 bands at the search range above M136279841
    hdgl_sieve_retarget(&sieve, LARGE_TARGETS, N_LARGE_TARGETS);

    printf("Sieve allocated and retargeted.\n");
    printf("Bands:\n");
    for (int i = 0; i < N_LARGE_TARGETS; i++) {
        printf("  Band %d: center p=%u  (slots [%u, %u])\n",
               i, LARGE_TARGETS[i],
               LARGE_TARGETS[i] - 1024, LARGE_TARGETS[i] + 1024);
    }
    printf("\n");

    cudaStream_t stream_sieve;
    CUDA_CHECK(cudaStreamCreate(&stream_sieve));

    uint32_t found_proxy[256];
    sh_timespec_t t0, t1;
    sh_gettime(&t0);
    long long total_steps = 0;
    long long total_cands = 0;

    for (int cycle = 0; cycle < total_cycles && g_running; cycle++) {

        // Step sieve
        for (int s = 0; s < STEPS_PER_CYCLE; s++)
            hdgl_sieve_step(&sieve, stream_sieve);
        CUDA_CHECK(cudaStreamSynchronize(stream_sieve));
        total_steps += STEPS_PER_CYCLE;

        // Harvest resonance-passing candidates from proxy ring
        int n = hdgl_sieve_harvest(&sieve, found_proxy, 256, stream_sieve);
        total_cands += n;

        // Record all resonance hits
        for (int i = 0; i < n; i++) {
            uint32_t p = found_proxy[i];
            if (p > 136279841u)
                record_candidate(p, 0.0f);  // score updated below
        }

        // Critic training: use resonance-only signal
        // Feature: [0=residue_unknown, p_norm, 0, 0, 0]
        if (n > 0) {
            float feat[5];
            for (int i = 0; i < n && i < 8; i++) {
                feat[0] = 0.5f;
                feat[1] = (float)found_proxy[i] / 140000000.0f;
                feat[2] = 0.5f;
                feat[3] = 0.5f;
                feat[4] = 0.5f;
                critic_observe(feat, 0.5f);   // neutral signal (no LL to anchor to)
            }
        }

        // Progress report
        if ((cycle + 1) % REPORT_EVERY == 0) {
            sh_gettime(&t1);
            double elapsed = sh_elapsed(&t0, &t1);
            double rate = (double)(total_steps) / elapsed;
            sh_gettime(&t0);
            total_steps = 0;

            printf("[cycle %5d/%d]  %.1fk steps/s  "
                   "proxy_cands=%lld  unique_hits=%d\n",
                   cycle + 1, total_cycles,
                   rate / 1000.0, total_cands, g_n_top);
        }
    }

    // Save critic state for future runs
    critic_save_weights("hdgl_critic_search.bin");

    // ---- Final report -------------------------------------------------------

    qsort(g_top, g_n_top, sizeof(TopCand), cmp_top);

    printf("\n====================================================================\n");
    printf("  SEARCH COMPLETE\n");
    printf("  Unique resonance candidates above M136279841: %d\n", g_n_top);
    printf("====================================================================\n");

    if (g_n_top > 0) {
        printf("\nTop candidates by resonance hit frequency:\n");
        printf("  %-14s  %-8s  %s\n", "p", "hits", "note");
        printf("  %-14s  %-8s  %s\n", "-----------", "----", "----");
        int show = (g_n_top < 30) ? g_n_top : 30;
        for (int i = 0; i < show; i++) {
            printf("  %-14u  %-8d  %s\n",
                   g_top[i].p, g_top[i].hits,
                   (g_top[i].hits >= 5) ? "<-- strong candidate" : "");
        }
        printf("\nTo verify: submit top candidates to GIMPS (mersenne.org)\n");
        printf("Full Lucas-Lehmer for p~136M requires ~months on a GPU cluster.\n");
    } else {
        printf("No candidates found in this run. Try more cycles.\n");
    }

    hdgl_sieve_free(&sieve);
    cudaStreamDestroy(stream_sieve);
    return 0;
}
