// ============================================================================
// HDGL Phase-LL Sieve Kernel v35
// Collapses LL into phase dynamics + Kuramoto coherence
// ============================================================================

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

#define TWO_PI 6.28318530717958647f

// Tunables
#define PHASE_ITERS        8
#define K_COUPLING         0.6f
#define GAMMA_DAMP         0.12f
#define ORACLE_THRESH      0.92f
#define PERSIST_THRESH     3.0f
#define LL_VERIFY_STEPS    64
#define RING_SIZE          256

// ============================================================================
// Device State
// ============================================================================

typedef struct {
    uint32_t *assigned_p;
    float    *phase;
    float    *phase_vel;
    float    *persist;
    uint32_t *prime_found;
    int      *prime_count;
    int       N;
    int       S;
    float     dt;
} PhaseDev;

__device__ __constant__ PhaseDev g_dev;

// ============================================================================
// Math helpers
// ============================================================================

__device__ __forceinline__ float fracf(float x) {
    return x - floorf(x);
}

// Λₙ = log_φ(n log n)
__device__ __forceinline__ float lambda_n(uint32_t n)
{
    const float phi = 1.61803398875f;
    float ln   = logf((float)n);
    float lnln = logf(ln);
    return (ln + lnln) / logf(phi);
}

// intrinsic phase from Λₙ
__device__ __forceinline__ float theta_from_lambda(float Lambda)
{
    return TWO_PI * fracf(Lambda);
}

// ============================================================================
// Phase-LL surrogate (quadratic map in phase space)
// ============================================================================

__device__ __forceinline__ float phase_ll_step(float theta, float theta_p)
{
    float t = 2.0f * theta;

    // wrap
    t = t - TWO_PI * floorf(t / TWO_PI);

    // nonlinear correction (encodes -2 and modular folding)
    float corr =
        -0.08f * sinf(theta - theta_p)
        + 0.02f * sinf(2.0f * theta);

    return t + corr;
}

// proximity to LL fixed point (θ ≈ 0 mod 2π)
__device__ __forceinline__ float phase_error(float theta)
{
    return fabsf(sinf(0.5f * theta));
}

// ============================================================================
// Minimal LL verifier (32-bit folded)
// ============================================================================

__device__ __forceinline__ uint32_t ll_step32(uint32_t s)
{
    uint64_t sq = (uint64_t)s * s;
    uint64_t r  = (sq - 2ULL);

    uint32_t lo = (uint32_t)r;
    uint32_t hi = (uint32_t)(r >> 32);

    uint32_t x = lo + hi;
    if (x < lo) x++; // carry fold

    return x;
}

// ============================================================================
// Main Kernel
// ============================================================================

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

    // --- Load state ---
    uint32_t p = g_dev.assigned_p[i];
    float theta = g_dev.phase[i];
    float vel   = g_dev.phase_vel[i];
    float acc   = g_dev.persist[i];

    int N = g_dev.N;
    int S = g_dev.S;
    float dt = g_dev.dt;

    // --- Spectral phase anchor ---
    float Lambda = lambda_n(p);
    float theta_p = theta_from_lambda(Lambda);

    // --- Phase-LL evolution ---
    #pragma unroll
    for (int k = 0; k < PHASE_ITERS; k++) {
        theta = phase_ll_step(theta, theta_p);
    }

    // --- Neighbour indices ---
    int ni[4] = {
        (i - 1 + N) % N,
        (i + 1) % N,
        (i - S + N) % N,
        (i + S) % N
    };

    // --- Kuramoto coupling ---
    float coupling = 0.0f;
    float sum_cos = 0.0f;
    float sum_sin = 0.0f;

    #pragma unroll
    for (int j = 0; j < 4; j++) {
        int nj = ni[j];
        float thj = g_dev.phase[nj];

        coupling += sinf(thj - theta);
        sum_cos  += cosf(thj);
        sum_sin  += sinf(thj);
    }

    // --- Second-order phase dynamics ---
    float omega_n = theta_p;  // intrinsic frequency

    float acc_theta = omega_n + K_COUPLING * coupling - GAMMA_DAMP * vel;

    vel   += dt * acc_theta;
    theta += dt * vel;

    // wrap
    theta = theta - TWO_PI * floorf(theta / TWO_PI);

    // --- Local order parameter ---
    float R = sqrtf(sum_cos * sum_cos + sum_sin * sum_sin) * 0.25f;

    // --- Phase LL error ---
    float err = phase_error(theta);

    // --- Oracle ---
    float oracle = R * __expf(-5.0f * err);

    // --- Persistence ---
    acc = 0.95f * acc + oracle;

    // ========================================================================
    // LL VERIFY (rare path)
    // ========================================================================

    if (oracle > ORACLE_THRESH && acc > PERSIST_THRESH)
    {
        uint32_t s = 4u;

        #pragma unroll
        for (int k = 0; k < LL_VERIFY_STEPS; k++) {
            s = ll_step32(s);
            if (s == 0u) break;
        }

        if (s == 0u) {
            int pos = atomicAdd(g_dev.prime_count, 1) % RING_SIZE;
            g_dev.prime_found[pos] = p;
        }

        // reset persistence to avoid spam
        acc = 0.0f;
    }

    // --- Write back ---
    g_dev.phase[i]     = theta;
    g_dev.phase_vel[i] = vel;
    g_dev.persist[i]   = acc;
}