/* secp256k1-batch-mul.cu -- bend form A: batched scalar*base-point on
 *                            secp256k1 via CUDA.
 *
 * Plan: examples/cuda-fanout/plans/form-A-secp256k1-batched-mul.md
 * Catalog: CATALOG.md form A (cuda-secp256k1-batched-mul).
 *
 * Wire (BSCP input, little-endian u32 headers, big-endian field bytes):
 *   "BSCP"           4 B magic
 *   u32 op_id        0x01 = scalar*base
 *   u32 n            number of scalars
 *   u8[32]  base_x   base point X, big-endian (secp256k1 native)
 *   u8[32]  base_y   base point Y, big-endian
 *   u8[n*32] scalars big-endian 256-bit ints
 *
 * Wire (BSCR output):
 *   "BSCR"           4 B magic on success, "BERR" on failure
 *   u32 status       0 ok / 1 bad-magic / 2 bad-n / 3 cuda-fail /
 *                    4 truncated / 5 bad-op / 6 base-not-on-curve
 *   u32 n            echoed scalar count
 *   u8[n*64] points  N affine points, each (X32 || Y32) big-endian.
 *                    Infinity encoded as all-zero (64 NUL bytes).
 *
 * Math: per thread compute k*P via Jacobian double-and-add
 *   (Q = O; for bit msb..lsb: Q = 2*Q; if bit set: Q = Q + P)
 * Field arithmetic uses GPUMath.h verbatim from FixedPaul/VanitySearch-Bitcrack
 * (vendored under vendor/vanity-search-bitcrack/, AGPL-3.0, original
 * Copyright (c) 2019 Jean Luc PONS preserved).
 *
 * After the per-thread Jacobian result lands in (X:Y:Z), we convert to
 * affine via Z^-1, Z^-2, Z^-3. Day-1 ships per-thread _ModInv; Day-2
 * may swap in _ModInvGrouped across N threads for amortisation.
 *
 * Daemon protocol mirrors shake256-fanout / cgbn-batch-worker:
 *   ready                      ; once at startup
 *   process-bin <in> <out>     ; binary request
 *   done <out>                 ; success
 *   error <reason>             ; failure
 *   quit / bye                 ; shutdown
 *
 * CLI: --binary <in> <out>     one-shot mode for tests
 *      --daemon                stdin loop, used by gpu-worker.lsp
 *
 * License: AGPLv3.
 */

#include <cuda_runtime.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <time.h>

#include "vendor/vanity-search-bitcrack/GPUMath.h"

#define CUDA_CHECK(stmt) do { \
    cudaError_t err = (stmt); \
    if (err != cudaSuccess) { \
        fprintf(stderr, "CUDA error %s at %s:%d: %s\n", \
                #stmt, __FILE__, __LINE__, cudaGetErrorString(err)); \
        return -1; \
    } \
} while (0)

enum {
    OP_SCALAR_BASE_MUL = 0x01,
};

enum {
    STATUS_OK            = 0,
    STATUS_BAD_MAGIC     = 1,
    STATUS_BAD_N         = 2,
    STATUS_CUDA_FAIL     = 3,
    STATUS_TRUNCATED     = 4,
    STATUS_BAD_OP        = 5,
    STATUS_NOT_ON_CURVE  = 6,
};

/* ------------------------------------------------------------------ */
/* Modular helpers we need but GPUMath.h does not provide.            */

/* Add two field elements mod p. r = a + b mod p. a,b assumed < p.
 * GPUMath.h ships ModSub but no ModAdd, so we synthesise one via a
 * conditional subtract after a 256-bit addition.
 */
__device__ __forceinline__ void ModAdd256(uint64_t* r, uint64_t* a, uint64_t* b) {
    /* a, b assumed < p. Sum fits in 257 bits. Subtract p once if the
     * 257-bit sum >= p; the conditional is (carry == 1) OR
     * (carry == 0 AND sum >= p), which equivalently is "subtract-p
     * produces no borrow once we extend the carry into bit 256".
     */
    uint64_t s[4];
    uint64_t carry;
    UADDO(s[0], a[0], b[0]);
    UADDC(s[1], a[1], b[1]);
    UADDC(s[2], a[2], b[2]);
    UADDC(s[3], a[3], b[3]);
    UADD(carry, 0ULL, 0ULL);
    /* Trial: s - p */
    uint64_t t[4]; uint64_t borrow;
    USUBO(t[0], s[0], 0xFFFFFFFEFFFFFC2FULL);
    USUBC(t[1], s[1], 0xFFFFFFFFFFFFFFFFULL);
    USUBC(t[2], s[2], 0xFFFFFFFFFFFFFFFFULL);
    USUBC(t[3], s[3], 0xFFFFFFFFFFFFFFFFULL);
    USUB(borrow, 0ULL, 0ULL);
    /* Take subtracted form if (carry produced) OR (no borrow). */
    if (carry || !borrow) {
        r[0] = t[0]; r[1] = t[1]; r[2] = t[2]; r[3] = t[3];
    } else {
        r[0] = s[0]; r[1] = s[1]; r[2] = s[2]; r[3] = s[3];
    }
}

/* in-place: r = 2*r mod p */
__device__ __forceinline__ void ModDouble256(uint64_t* r) {
    ModAdd256(r, r, r);
}

/* in-place: r = 3*r mod p */
__device__ __forceinline__ void ModTriple256(uint64_t* r, uint64_t* a) {
    uint64_t t[4];
    ModAdd256(t, a, a);
    ModAdd256(r, t, a);
}

/* Convenience wrapper around GPUMath.h's 5-limb _ModInv. The on-curve
 * inputs are 4 limbs (< p), so we pad to 5 limbs (limb[4]=0) and
 * accept the in-place result.
 */
__device__ __forceinline__ void ModInv256(uint64_t* r) {
    uint64_t v[5];
    v[0] = r[0]; v[1] = r[1]; v[2] = r[2]; v[3] = r[3]; v[4] = 0;
    _ModInv(v);
    r[0] = v[0]; r[1] = v[1]; r[2] = v[2]; r[3] = v[3];
}

/* Convenience wrapper around GPUMath.h's 5-limb _ModMult / _ModSqr.
 * Limb 4 is set to 0 for the inputs since they are < p; the result
 * limb 4 is also < 1 because _ModMult does the final fold.
 */
/* Final-reduce a 4-limb value possibly in [0, 2p) down to [0, p). */
__device__ __forceinline__ void cond_sub_p(uint64_t* r) {
    uint64_t t[4]; uint64_t borrow;
    USUBO(t[0], r[0], 0xFFFFFFFEFFFFFC2FULL);
    USUBC(t[1], r[1], 0xFFFFFFFFFFFFFFFFULL);
    USUBC(t[2], r[2], 0xFFFFFFFFFFFFFFFFULL);
    USUBC(t[3], r[3], 0xFFFFFFFFFFFFFFFFULL);
    USUB(borrow, 0ULL, 0ULL);
    if (!borrow) { r[0]=t[0]; r[1]=t[1]; r[2]=t[2]; r[3]=t[3]; }
}

__device__ __forceinline__ void ModMult4(uint64_t* r, uint64_t* a, uint64_t* b) {
    /* _ModMult ignores aa[4]/bb[4] but writes only rr[0..3]. */
    uint64_t aa[5], bb[5], rr[5] = {0,0,0,0,0};
    aa[0] = a[0]; aa[1] = a[1]; aa[2] = a[2]; aa[3] = a[3]; aa[4] = 0;
    bb[0] = b[0]; bb[1] = b[1]; bb[2] = b[2]; bb[3] = b[3]; bb[4] = 0;
    _ModMult(rr, aa, bb);
    r[0] = rr[0]; r[1] = rr[1]; r[2] = rr[2]; r[3] = rr[3];
    cond_sub_p(r);
}

__device__ __forceinline__ void ModSqr4(uint64_t* r, uint64_t* a) {
    uint64_t aa[5], rr[5] = {0,0,0,0,0};
    aa[0] = a[0]; aa[1] = a[1]; aa[2] = a[2]; aa[3] = a[3]; aa[4] = 0;
    _ModSqr(rr, aa);
    r[0] = rr[0]; r[1] = rr[1]; r[2] = rr[2]; r[3] = rr[3];
    cond_sub_p(r);
}

/* ------------------------------------------------------------------ */
/* Jacobian point ops on secp256k1 (a = 0). Coordinates X,Y,Z each
 * 4-limb little-endian field elements. Identity = (any, any, 0).
 */

__device__ __forceinline__ int is_zero4(uint64_t* a) {
    return (a[0] | a[1] | a[2] | a[3]) == 0;
}

/* Jacobian double (a=0):
 *   A = X^2; B = Y^2; C = B^2
 *   D = 2*((X+B)^2 - A - C)
 *   E = 3*A; F = E^2
 *   X3 = F - 2*D
 *   Y3 = E*(D - X3) - 8*C
 *   Z3 = 2*Y*Z
 */
__device__ void jac_double(uint64_t* X3, uint64_t* Y3, uint64_t* Z3,
                            uint64_t* X1, uint64_t* Y1, uint64_t* Z1) {
    if (is_zero4(Z1)) {
        X3[0]=X3[1]=X3[2]=X3[3]=0;
        Y3[0]=Y3[1]=Y3[2]=Y3[3]=0;
        Z3[0]=Z3[1]=Z3[2]=Z3[3]=0;
        return;
    }
    uint64_t A[4], B[4], C[4], D[4], E[4], F[4], XB[4], t[4];
    ModSqr4(A, X1);             /* A = X^2 */
    ModSqr4(B, Y1);             /* B = Y^2 */
    ModSqr4(C, B);              /* C = B^2 */
    ModAdd256(XB, X1, B);       /* XB = X + B */
    ModSqr4(t, XB);             /* t = (X+B)^2 */
    ModSub256(t, A);            /* t -= A */
    ModSub256(t, C);            /* t -= C */
    ModDouble256(t);            /* D = 2*t */
    D[0]=t[0]; D[1]=t[1]; D[2]=t[2]; D[3]=t[3];
    ModTriple256(E, A);         /* E = 3*A */
    ModSqr4(F, E);              /* F = E^2 */
    /* X3 = F - 2*D */
    uint64_t twoD[4];
    twoD[0]=D[0]; twoD[1]=D[1]; twoD[2]=D[2]; twoD[3]=D[3];
    ModDouble256(twoD);
    ModSub256(X3, F, twoD);
    /* Y3 = E*(D - X3) - 8*C */
    uint64_t DmX[4], eightC[4], EDX[4];
    ModSub256(DmX, D, X3);
    ModMult4(EDX, E, DmX);
    eightC[0]=C[0]; eightC[1]=C[1]; eightC[2]=C[2]; eightC[3]=C[3];
    ModDouble256(eightC);
    ModDouble256(eightC);
    ModDouble256(eightC);
    ModSub256(Y3, EDX, eightC);
    /* Z3 = 2*Y*Z */
    uint64_t YZ[4];
    ModMult4(YZ, Y1, Z1);
    ModAdd256(Z3, YZ, YZ);
}

/* Jacobian + affine addition (a=0). P1 in Jacobian, P2 in affine.
 * (X2, Y2) affine. Handles P1=O and P1+/-P2 special cases.
 */
__device__ void jac_add_affine(uint64_t* X3, uint64_t* Y3, uint64_t* Z3,
                                uint64_t* X1, uint64_t* Y1, uint64_t* Z1,
                                uint64_t* X2, uint64_t* Y2) {
    if (is_zero4(Z1)) {
        X3[0]=X2[0]; X3[1]=X2[1]; X3[2]=X2[2]; X3[3]=X2[3];
        Y3[0]=Y2[0]; Y3[1]=Y2[1]; Y3[2]=Y2[2]; Y3[3]=Y2[3];
        Z3[0]=1; Z3[1]=0; Z3[2]=0; Z3[3]=0;
        return;
    }
    uint64_t Z1Z1[4], U2[4], S2[4], H[4], r[4], HH[4], HHH[4], V[4];
    uint64_t Z1Z1Z1[4], t1[4], t2[4];
    ModSqr4(Z1Z1, Z1);              /* Z1Z1 = Z1^2 */
    ModMult4(U2, X2, Z1Z1);         /* U2 = X2*Z1^2 */
    ModMult4(Z1Z1Z1, Z1, Z1Z1);     /* Z1^3 */
    ModMult4(S2, Y2, Z1Z1Z1);       /* S2 = Y2*Z1^3 */
    ModSub256(H, U2, X1);           /* H = U2 - X1 */
    ModSub256(r, S2, Y1);           /* r = S2 - Y1 */
    if (is_zero4(H)) {
        /* Same X. If r == 0 -> point double; else -> infinity */
        if (is_zero4(r)) {
            jac_double(X3, Y3, Z3, X1, Y1, Z1);
            return;
        }
        X3[0]=X3[1]=X3[2]=X3[3]=0;
        Y3[0]=Y3[1]=Y3[2]=Y3[3]=0;
        Z3[0]=Z3[1]=Z3[2]=Z3[3]=0;
        return;
    }
    ModSqr4(HH, H);                 /* HH = H^2 */
    ModMult4(HHH, H, HH);           /* HHH = H^3 */
    ModMult4(V, X1, HH);            /* V = X1*HH */
    /* X3 = r^2 - HHH - 2*V */
    uint64_t r2[4], twoV[4];
    ModSqr4(r2, r);
    twoV[0]=V[0]; twoV[1]=V[1]; twoV[2]=V[2]; twoV[3]=V[3];
    ModDouble256(twoV);
    ModSub256(t1, r2, HHH);
    ModSub256(X3, t1, twoV);
    /* Y3 = r*(V - X3) - Y1*HHH */
    ModSub256(t1, V, X3);
    ModMult4(t2, r, t1);
    ModMult4(t1, Y1, HHH);
    ModSub256(Y3, t2, t1);
    /* Z3 = Z1 * H */
    ModMult4(Z3, Z1, H);
}

/* ------------------------------------------------------------------ */
/* Scan a big-endian 32-byte scalar bit-by-bit, top to bottom.        */

__device__ __forceinline__ int scalar_bit_be(const uint8_t* k, int i) {
    /* bit 0 = LSB. i in [0, 255]. k is BE: k[0] = top byte. */
    int byte_idx = 31 - (i >> 3);
    int bit_idx  = i & 7;
    return (k[byte_idx] >> bit_idx) & 1;
}

/* Per-thread: scalar k (big-endian) * base point (Px, Py, affine) ->
 * Jacobian (X, Y, Z). Output in 4-limb LE.
 */
__device__ void scalar_mul(uint64_t* outX, uint64_t* outY, uint64_t* outZ,
                           const uint8_t* k_be,
                           uint64_t* Px, uint64_t* Py) {
    /* Q = O (Z=0) */
    uint64_t Qx[4] = {0,0,0,0};
    uint64_t Qy[4] = {0,0,0,0};
    uint64_t Qz[4] = {0,0,0,0};
    uint64_t Tx[4], Ty[4], Tz[4];
    for (int i = 255; i >= 0; i--) {
        jac_double(Tx, Ty, Tz, Qx, Qy, Qz);
        Qx[0]=Tx[0]; Qx[1]=Tx[1]; Qx[2]=Tx[2]; Qx[3]=Tx[3];
        Qy[0]=Ty[0]; Qy[1]=Ty[1]; Qy[2]=Ty[2]; Qy[3]=Ty[3];
        Qz[0]=Tz[0]; Qz[1]=Tz[1]; Qz[2]=Tz[2]; Qz[3]=Tz[3];
        if (scalar_bit_be(k_be, i)) {
            jac_add_affine(Tx, Ty, Tz, Qx, Qy, Qz, Px, Py);
            Qx[0]=Tx[0]; Qx[1]=Tx[1]; Qx[2]=Tx[2]; Qx[3]=Tx[3];
            Qy[0]=Ty[0]; Qy[1]=Ty[1]; Qy[2]=Ty[2]; Qy[3]=Ty[3];
            Qz[0]=Tz[0]; Qz[1]=Tz[1]; Qz[2]=Tz[2]; Qz[3]=Tz[3];
        }
    }
    outX[0]=Qx[0]; outX[1]=Qx[1]; outX[2]=Qx[2]; outX[3]=Qx[3];
    outY[0]=Qy[0]; outY[1]=Qy[1]; outY[2]=Qy[2]; outY[3]=Qy[3];
    outZ[0]=Qz[0]; outZ[1]=Qz[1]; outZ[2]=Qz[2]; outZ[3]=Qz[3];
}

/* ------------------------------------------------------------------ */
/* Emit a 4-limb LE field element as 32 BE bytes at outp[0..32).      */
__device__ __forceinline__ void limbs_to_be32(uint8_t* outp, const uint64_t* a) {
    for (int li = 0; li < 4; li++) {
        uint64_t v = a[3 - li];
        for (int b = 0; b < 8; b++) {
            outp[li*8 + b] = (uint8_t)(v >> (56 - 8*b));
        }
    }
}

/* ------------------------------------------------------------------ */
/* Kernel v1: one thread per scalar, per-thread _ModInv (Day-1).      */

__global__ void secp_mul_batch_kernel_v1(
    const uint64_t* base_xy,   /* 8 limbs total: [Px(4), Py(4)] LE */
    const uint8_t*  scalars,   /* n * 32 bytes BE */
    uint8_t*        out_points,/* n * 64 bytes BE: x32 || y32 */
    uint32_t        n)
{
    uint32_t idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= n) return;
    const uint8_t* k_be = scalars + (size_t)idx * 32;

    uint64_t Px[4] = { base_xy[0], base_xy[1], base_xy[2], base_xy[3] };
    uint64_t Py[4] = { base_xy[4], base_xy[5], base_xy[6], base_xy[7] };

    uint64_t Jx[4], Jy[4], Jz[4];
    scalar_mul(Jx, Jy, Jz, k_be, Px, Py);

    uint8_t* outp = out_points + (size_t)idx * 64;

    if (is_zero4(Jz)) {
        /* infinity: emit 64 NUL bytes */
        for (int i = 0; i < 64; i++) outp[i] = 0;
        return;
    }

    /* Affine = (X/Z^2, Y/Z^3). Inverse Z. */
    uint64_t Zi[4]  = { Jz[0], Jz[1], Jz[2], Jz[3] };
    ModInv256(Zi);
    uint64_t Zi2[4], Zi3[4];
    ModSqr4(Zi2, Zi);
    ModMult4(Zi3, Zi2, Zi);
    uint64_t Ax[4], Ay[4];
    ModMult4(Ax, Jx, Zi2);
    ModMult4(Ay, Jy, Zi3);

    limbs_to_be32(outp, Ax);
    limbs_to_be32(outp + 32, Ay);
}

/* ------------------------------------------------------------------ */
/* Kernel v3: windowed-base ladder using precomputed table of i*P.    */
/*                                                                    */
/* Day-2 progress doc identified scalar_mul (256 doubles + ~128 adds  */
/* per scalar) as the kernel's dominant cost — not ModInv. v3 cuts    */
/* that down with a window-w ladder:                                  */
/*                                                                    */
/*   Precompute table[0..(2^w - 1)]: affine (x,y) of i*P, where P is  */
/*   the per-request base point (typically G).                        */
/*   For each input scalar k:                                          */
/*     Q = O                                                           */
/*     for window from high to low (256/w windows):                    */
/*       Q = (2^w) * Q                                                 */
/*       w_val = scalar bits [w*window .. w*window + w)                */
/*       Q = Q + table[w_val]                                          */
/*                                                                     */
/* Costs per scalar with w=4 (64 windows):                            */
/*   256 doubles (same as binary)                                     */
/*   ~63 adds (vs ~128 in binary double-and-add): saves ~65 adds.     */
/*                                                                     */
/* Table built once per request via secp_window_table_init kernel.    */
/* 16 entries × 64 bytes (Ax,Ay) = 1024 B for w=4; 16 KB for w=8.     */
/* Table loaded cooperatively into __shared__ at v3 kernel block      */
/* start so per-thread lookups hit shared, not global (avoids the     */
/* broadcast-only constant-memory penalty for divergent indices).     */
/*                                                                     */
/* Table-init algorithm (1 launch, 1 thread for w=4): walk i = 1..15  */
/* in Jacobian, converting each to affine via one ModInv. We use the  */
/* "build incrementally then batch-invert Z's" trick to cut 15 ModInv */
/* to 1, then 14 ModMult. Output is 16 affine (x,y) pairs in LE limbs */
/* (entry 0 = (0,0) sentinel for "infinity / skip add"; entries 1..15 */
/* = i*P affine).                                                     */

__device__ __forceinline__ int affine_is_zero(const uint64_t* xy) {
    return ((xy[0] | xy[1] | xy[2] | xy[3]) == 0)
        && ((xy[4] | xy[5] | xy[6] | xy[7]) == 0);
}

__global__ void v3_table_init(
    const uint64_t* base_xy,   /* 8 limbs: [Px(4), Py(4)] LE          */
    uint64_t*       table,     /* W * 8 u64: W affine (Ax,Ay) entries */
    uint32_t        W)         /* 2^window_w; we expect 16 or 256     */
{
    /* Single-thread builder — runs once per request; cost is constant  */
    /* in N. We compute (i*P) for i in 1..W-1 in Jacobian, collect all  */
    /* Z's, then do a single batch inversion across those W-1 Z's, and  */
    /* finalise each entry's affine (x,y).                              */
    if (blockIdx.x != 0 || threadIdx.x != 0) return;

    uint64_t Px[4] = { base_xy[0], base_xy[1], base_xy[2], base_xy[3] };
    uint64_t Py[4] = { base_xy[4], base_xy[5], base_xy[6], base_xy[7] };

    /* Entry 0: sentinel "infinity" — emit (0, 0). Window value 0 is    */
    /* a no-op add in the kernel.                                       */
    table[0] = 0; table[1] = 0; table[2] = 0; table[3] = 0;
    table[4] = 0; table[5] = 0; table[6] = 0; table[7] = 0;

    if (W <= 1) return;

    /* Storage caps at W <= 16 to keep stack bounded — Day-3 ships only */
    /* w=4 (W=16). Larger W must be added with a host-provided scratch. */
    const uint32_t MAX_W = 16;
    if (W > MAX_W) return;   /* unsupported in this kernel */

    uint64_t Jx[MAX_W][4], Jy[MAX_W][4], Jz[MAX_W][4];

    /* i = 1: Q = P. */
    Jx[1][0]=Px[0]; Jx[1][1]=Px[1]; Jx[1][2]=Px[2]; Jx[1][3]=Px[3];
    Jy[1][0]=Py[0]; Jy[1][1]=Py[1]; Jy[1][2]=Py[2]; Jy[1][3]=Py[3];
    Jz[1][0]=1; Jz[1][1]=0; Jz[1][2]=0; Jz[1][3]=0;

    for (uint32_t i = 2; i < W; i++) {
        jac_add_affine(Jx[i], Jy[i], Jz[i],
                       Jx[i-1], Jy[i-1], Jz[i-1],
                       Px, Py);
    }

    /* Batch-invert the W-1 Z's: classic Montgomery trick.              */
    /* prefix[i] = Z[1]*Z[2]*...*Z[i] for i in 1..W-1.                  */
    uint64_t prefix[MAX_W][4];
    prefix[1][0]=Jz[1][0]; prefix[1][1]=Jz[1][1];
    prefix[1][2]=Jz[1][2]; prefix[1][3]=Jz[1][3];
    for (uint32_t i = 2; i < W; i++) {
        ModMult4(prefix[i], prefix[i-1], Jz[i]);
    }
    /* u = (prefix[W-1])^-1 = (Π Z[i])^-1                               */
    uint64_t u[4] = { prefix[W-1][0], prefix[W-1][1],
                      prefix[W-1][2], prefix[W-1][3] };
    ModInv256(u);

    /* Walk i = W-1 .. 2: Zi[i] = u * prefix[i-1]; u <- u * Z[i]        */
    /* For i = 1: Zi[1] = u.                                             */
    uint64_t Zi[MAX_W][4];
    for (uint32_t ii = W - 1; ii >= 2; ii--) {
        ModMult4(Zi[ii], u, prefix[ii-1]);
        uint64_t newu[4];
        ModMult4(newu, u, Jz[ii]);
        u[0]=newu[0]; u[1]=newu[1]; u[2]=newu[2]; u[3]=newu[3];
    }
    Zi[1][0]=u[0]; Zi[1][1]=u[1]; Zi[1][2]=u[2]; Zi[1][3]=u[3];

    /* Finalise affine: Ax = Jx * Zi^2; Ay = Jy * Zi^3.                 */
    for (uint32_t i = 1; i < W; i++) {
        uint64_t Zi2[4], Zi3[4], Ax[4], Ay[4];
        ModSqr4(Zi2, Zi[i]);
        ModMult4(Zi3, Zi2, Zi[i]);
        ModMult4(Ax, Jx[i], Zi2);
        ModMult4(Ay, Jy[i], Zi3);
        uint64_t* slot = table + (size_t)i * 8;
        slot[0]=Ax[0]; slot[1]=Ax[1]; slot[2]=Ax[2]; slot[3]=Ax[3];
        slot[4]=Ay[0]; slot[5]=Ay[1]; slot[6]=Ay[2]; slot[7]=Ay[3];
    }
}

/* Read a w-bit window from BE scalar k_be at window index `widx`,
 * where widx 0 == top window. Returns int in [0, 2^w - 1].
 * Scalar has total_bits = 256; W = 2^w; num_windows = 256/w.
 */
__device__ __forceinline__ uint32_t scalar_window_be(
    const uint8_t* k_be, uint32_t widx, uint32_t w)
{
    /* Top window covers bits [256-w .. 256). We walk MSB->LSB.
     * bit index of top of window = 256 - widx*w - 1
     * bit index of bottom        = 256 - widx*w - w
     * For w in {1,2,4,8} and aligned to byte boundary where possible.
     */
    uint32_t top_bit = 256u - widx * w;        /* exclusive */
    uint32_t bot_bit = top_bit - w;            /* inclusive */
    /* Read bits [bot_bit, top_bit) from k_be (BE 32 bytes).            */
    uint32_t v = 0;
    for (uint32_t b = 0; b < w; b++) {
        uint32_t bit_idx = bot_bit + b;        /* LSB = 0 */
        uint32_t byte_idx = 31u - (bit_idx >> 3);
        uint32_t bit_in_byte = bit_idx & 7u;
        uint32_t bit = (k_be[byte_idx] >> bit_in_byte) & 1u;
        v |= (bit << b);
    }
    return v;
}

/* Per-thread windowed scalar mul. Reads affine table[0..W-1] from
 * shared memory (preloaded by the block). Returns Jacobian (X,Y,Z).
 */
__device__ void scalar_mul_windowed(
    uint64_t* outX, uint64_t* outY, uint64_t* outZ,
    const uint8_t* k_be,
    const uint64_t* table_shared,   /* W * 8 u64 */
    uint32_t w,
    uint32_t W)
{
    uint64_t Qx[4] = {0,0,0,0};
    uint64_t Qy[4] = {0,0,0,0};
    uint64_t Qz[4] = {0,0,0,0};
    uint64_t Tx[4], Ty[4], Tz[4];

    uint32_t num_windows = 256u / w;
    for (uint32_t widx = 0; widx < num_windows; widx++) {
        /* Q = 2^w * Q : w doubles. */
        for (uint32_t d = 0; d < w; d++) {
            jac_double(Tx, Ty, Tz, Qx, Qy, Qz);
            Qx[0]=Tx[0]; Qx[1]=Tx[1]; Qx[2]=Tx[2]; Qx[3]=Tx[3];
            Qy[0]=Ty[0]; Qy[1]=Ty[1]; Qy[2]=Ty[2]; Qy[3]=Ty[3];
            Qz[0]=Tz[0]; Qz[1]=Tz[1]; Qz[2]=Tz[2]; Qz[3]=Tz[3];
        }
        uint32_t wv = scalar_window_be(k_be, widx, w);
        if (wv == 0) continue;   /* table[0] = O; skip the add */
        const uint64_t* slot = table_shared + (size_t)wv * 8;
        uint64_t Px[4] = { slot[0], slot[1], slot[2], slot[3] };
        uint64_t Py[4] = { slot[4], slot[5], slot[6], slot[7] };
        jac_add_affine(Tx, Ty, Tz, Qx, Qy, Qz, Px, Py);
        Qx[0]=Tx[0]; Qx[1]=Tx[1]; Qx[2]=Tx[2]; Qx[3]=Tx[3];
        Qy[0]=Ty[0]; Qy[1]=Ty[1]; Qy[2]=Ty[2]; Qy[3]=Ty[3];
        Qz[0]=Tz[0]; Qz[1]=Tz[1]; Qz[2]=Tz[2]; Qz[3]=Tz[3];
    }
    outX[0]=Qx[0]; outX[1]=Qx[1]; outX[2]=Qx[2]; outX[3]=Qx[3];
    outY[0]=Qy[0]; outY[1]=Qy[1]; outY[2]=Qy[2]; outY[3]=Qy[3];
    outZ[0]=Qz[0]; outZ[1]=Qz[1]; outZ[2]=Qz[2]; outZ[3]=Qz[3];
}

/* v3 kernel: per-thread windowed ladder + per-thread ModInv to affine. */
/* Table preloaded into __shared__ once per block.                      */
/* Templated on W to let the compiler unroll the shared-memory loader.  */
template <uint32_t W>
__global__ void secp_mul_batch_kernel_v3(
    const uint64_t* base_xy,   /* unused but kept for ABI parity */
    const uint64_t* table_g,   /* W * 8 u64 in global memory */
    const uint8_t*  scalars,
    uint8_t*        out_points,
    uint32_t        n,
    uint32_t        w)
{
    __shared__ uint64_t s_table[W * 8];
    /* Cooperative load: each thread copies (W*8 / blockDim.x) entries. */
    for (uint32_t i = threadIdx.x; i < W * 8; i += blockDim.x) {
        s_table[i] = table_g[i];
    }
    __syncthreads();

    uint32_t idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= n) return;

    const uint8_t* k_be = scalars + (size_t)idx * 32;

    uint64_t Jx[4], Jy[4], Jz[4];
    scalar_mul_windowed(Jx, Jy, Jz, k_be, s_table, w, W);

    uint8_t* outp = out_points + (size_t)idx * 64;

    if (is_zero4(Jz)) {
        for (int i = 0; i < 64; i++) outp[i] = 0;
        return;
    }

    /* Affine = (X/Z^2, Y/Z^3). Per-thread inversion (Day-3 baseline). */
    uint64_t Zi[4]  = { Jz[0], Jz[1], Jz[2], Jz[3] };
    ModInv256(Zi);
    uint64_t Zi2[4], Zi3[4];
    ModSqr4(Zi2, Zi);
    ModMult4(Zi3, Zi2, Zi);
    uint64_t Ax[4], Ay[4];
    ModMult4(Ax, Jx, Zi2);
    ModMult4(Ay, Jy, Zi3);

    limbs_to_be32(outp, Ax);
    limbs_to_be32(outp + 32, Ay);
}

/* ------------------------------------------------------------------ */
/* Kernel v4: Day-4 stack — v3 windowed-G ladder feeds the v2          */
/* Montgomery batch-inversion pipeline.                                */
/*                                                                     */
/* Pipeline:                                                           */
/*   Phase A (v4_phase_a<W>): per-thread windowed scalar_mul, writes   */
/*     Jacobian (Jx, Jy, Jz) to scratch[idx*12..idx*12+12) and writes  */
/*     z[idx] = Jz (or 1 if Jz == 0) to scratch z[].                   */
/*   Phase B (v2_phase_b reused): per-block prefix product over z[],   */
/*     writes c[] and block_prod[].                                    */
/*   Phase C (v2_phase_c reused): single-thread cross-block walk,      */
/*     ONE _ModInv on the full product, derives G_left[] + u_entry[].  */
/*   Phase D (v2_phase_d reused): per-block backward walk, emits       */
/*     affine BE bytes.                                                */
/*                                                                     */
/* The v3 win (cut scalar_mul cost by ~16% of field-mults) compounds   */
/* with the v2 win (collapse N parallel ModInvs to 1 + N-1 ModMults)   */
/* because Day-3 progress doc identified: "Once v3 lands, the          */
/* bottleneck shifts BACK toward scalar_mul's residual cost +          */
/* per-thread inversion." v4 retires both bottlenecks in one pipeline. */

template <uint32_t W>
__global__ void v4_phase_a(
    const uint64_t* base_xy,    /* unused but kept for ABI parity */
    const uint64_t* table_g,    /* W * 8 u64 in global memory */
    const uint8_t*  scalars,
    uint64_t*       jac,        /* n * 12 u64: [Jx(4), Jy(4), Jz(4)]  */
    uint64_t*       z,          /* n * 4  u64: Jz, or 1 if Jz == 0    */
    uint32_t        n,
    uint32_t        w)
{
    __shared__ uint64_t s_table[W * 8];
    for (uint32_t i = threadIdx.x; i < W * 8; i += blockDim.x) {
        s_table[i] = table_g[i];
    }
    __syncthreads();

    uint32_t idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= n) return;
    const uint8_t* k_be = scalars + (size_t)idx * 32;

    uint64_t Jx[4], Jy[4], Jz[4];
    scalar_mul_windowed(Jx, Jy, Jz, k_be, s_table, w, W);

    uint64_t* slot = jac + (size_t)idx * 12;
    slot[0]=Jx[0]; slot[1]=Jx[1]; slot[2]=Jx[2]; slot[3]=Jx[3];
    slot[4]=Jy[0]; slot[5]=Jy[1]; slot[6]=Jy[2]; slot[7]=Jy[3];
    slot[8]=Jz[0]; slot[9]=Jz[1]; slot[10]=Jz[2]; slot[11]=Jz[3];

    uint64_t* zslot = z + (size_t)idx * 4;
    if (is_zero4(Jz)) {
        zslot[0] = 1; zslot[1] = 0; zslot[2] = 0; zslot[3] = 0;
    } else {
        zslot[0] = Jz[0]; zslot[1] = Jz[1]; zslot[2] = Jz[2]; zslot[3] = Jz[3];
    }
}

/* ------------------------------------------------------------------ */
/* Kernel v2: Montgomery batch inversion across N independent z's.    */
/*                                                                    */
/* Pipeline (3 launches, separated by host-side barriers):            */
/*  Phase A: per-thread scalar_mul -> writes Jacobian (Jx,Jy,Jz) to   */
/*           scratch[ idx*12 .. idx*12+12 ) ; encodes infinity by     */
/*           leaving Jz = 0. Also writes z[idx] := Jz (or 1 if 0).    */
/*  Phase B: per-block prefix product. Each thread tid in block writes*/
/*           c[ block_base + tid ] = prod( z[block_base..block_base+tid] )*/
/*           Thread 0 in block writes block_prod[blockIdx] = full prod*/
/*  Phase C: single-thread kernel computes cross-block scan +         */
/*           one _ModInv on the global product, then walks backward   */
/*           propagating per-block start inverses.                    */
/*  Phase D: per-thread: derive z_inv from c[idx] and block-start inv,*/
/*           finish affine projection, emit BE bytes.                 */
/*                                                                    */
/* The math, with z[0..N-1] independent: define                       */
/*   c[i] = z[0]*z[1]*...*z[i]                                        */
/*   u    = (c[N-1])^-1                                               */
/* Then:                                                              */
/*   z_inv[N-1] = u * c[N-2]                                          */
/*   z_inv[i]   = u * c[i-1]; u <- u*z[i]   for i=N-2..1              */
/*   z_inv[0]   = u                                                   */
/* Per-block variant: each block reduces a contiguous slice; cross-   */
/* block reduce + one global inv runs in a tiny serial kernel.        */
/*                                                                    */
/* Note on infinities: z[i] == 0 (Jz of the result) means the answer  */
/* is the identity. We substitute z[i] := 1 in the prefix product so  */
/* the inversion stays well-defined, and emit 64 NUL bytes in Phase D */
/* by looking at the original Jz value in scratch.                    */

/* Phase A: scalar-mul, writing (Jx,Jy,Jz) to scratch.                */
__global__ void v2_phase_a(
    const uint64_t* base_xy,
    const uint8_t*  scalars,
    uint64_t*       jac,        /* n * 12 u64: [Jx(4), Jy(4), Jz(4)]  */
    uint64_t*       z,          /* n * 4  u64: Jz, or 1 if Jz == 0    */
    uint32_t        n)
{
    uint32_t idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= n) return;
    const uint8_t* k_be = scalars + (size_t)idx * 32;
    uint64_t Px[4] = { base_xy[0], base_xy[1], base_xy[2], base_xy[3] };
    uint64_t Py[4] = { base_xy[4], base_xy[5], base_xy[6], base_xy[7] };

    uint64_t Jx[4], Jy[4], Jz[4];
    scalar_mul(Jx, Jy, Jz, k_be, Px, Py);

    uint64_t* slot = jac + (size_t)idx * 12;
    slot[0]=Jx[0]; slot[1]=Jx[1]; slot[2]=Jx[2]; slot[3]=Jx[3];
    slot[4]=Jy[0]; slot[5]=Jy[1]; slot[6]=Jy[2]; slot[7]=Jy[3];
    slot[8]=Jz[0]; slot[9]=Jz[1]; slot[10]=Jz[2]; slot[11]=Jz[3];

    uint64_t* zslot = z + (size_t)idx * 4;
    if (is_zero4(Jz)) {
        zslot[0] = 1; zslot[1] = 0; zslot[2] = 0; zslot[3] = 0;
    } else {
        zslot[0] = Jz[0]; zslot[1] = Jz[1]; zslot[2] = Jz[2]; zslot[3] = Jz[3];
    }
}

/* Phase B: per-block prefix product across z[block_base..block_end).
 * c[i] = z[block_base] * z[block_base+1] * ... * z[i].
 * Sequential within block (block size kept modest). Thread 0 in block
 * does the work, the other threads in the block sit idle for this
 * phase — keeps the algorithm trivially correct. (We trade some block
 * parallelism for byte-identity with the host reference; the win comes
 * from collapsing N _ModInv calls into 1, not from intra-block scan.)
 */
__global__ void v2_phase_b(
    const uint64_t* z,          /* n * 4 u64 */
    uint64_t*       c,          /* n * 4 u64: per-thread prefix product */
    uint64_t*       block_prod, /* num_blocks * 4 u64 */
    uint32_t        n,
    uint32_t        block_size)
{
    uint32_t bid = blockIdx.x;
    if (threadIdx.x != 0) return;
    uint32_t base = bid * block_size;
    if (base >= n) return;
    uint32_t end = base + block_size;
    if (end > n) end = n;
    uint64_t acc[4];
    /* acc = z[base] */
    acc[0] = z[(size_t)base*4 + 0];
    acc[1] = z[(size_t)base*4 + 1];
    acc[2] = z[(size_t)base*4 + 2];
    acc[3] = z[(size_t)base*4 + 3];
    c[(size_t)base*4 + 0] = acc[0];
    c[(size_t)base*4 + 1] = acc[1];
    c[(size_t)base*4 + 2] = acc[2];
    c[(size_t)base*4 + 3] = acc[3];
    for (uint32_t i = base + 1; i < end; i++) {
        uint64_t zi[4] = { z[(size_t)i*4+0], z[(size_t)i*4+1],
                           z[(size_t)i*4+2], z[(size_t)i*4+3] };
        uint64_t r[4];
        ModMult4(r, acc, zi);
        acc[0]=r[0]; acc[1]=r[1]; acc[2]=r[2]; acc[3]=r[3];
        c[(size_t)i*4+0] = acc[0];
        c[(size_t)i*4+1] = acc[1];
        c[(size_t)i*4+2] = acc[2];
        c[(size_t)i*4+3] = acc[3];
    }
    block_prod[(size_t)bid*4 + 0] = acc[0];
    block_prod[(size_t)bid*4 + 1] = acc[1];
    block_prod[(size_t)bid*4 + 2] = acc[2];
    block_prod[(size_t)bid*4 + 3] = acc[3];
}

/* Phase C: single-thread kernel. Build global prefix product of
 * block_prod[], invert the last, walk backward to derive each block's
 * "u_block": the inverse of the product of z[0..block_base-1] *
 * prefix(block_prod). For block b, we need:
 *   u_after_block_b = ( z[0] * z[1] * ... * z[block_end_b - 1] )^-1
 * From this, in Phase D, for thread i in block b at position p within
 * block (p = i - block_base):
 *   z_inv[i] = u_after_block_b * (block_prod_full / c[i]) ... we use
 *   a simpler form. Re-derive cleanly:
 *
 * Let B = block_size, K = num_blocks.
 * For block b, c_local[base..base+B-1] is the per-block prefix.
 * block_prod[b] = c_local[base+B-1].
 *
 * Define G[b] = product of block_prod[0..b] (global prefix).
 *   G[K-1] is the global product over all z.
 * Single inv: invG = G[K-1]^-1.
 * For block b: u_b = invG * (G[K-1] / G[b]) = invG * product_{j>b} block_prod[j].
 *   Equivalently, walk b = K-1..0:
 *     u_{K-1} = invG  (since "after all blocks" is invG times product of nothing)
 *     ... wait, we need u_b such that u_b * (product over block b's z's) = u_{b-1}.
 *
 * Cleaner: u_b is "global inverse over product of z[0..base_b+B-1]".
 *   u_b = invG * product_{j>b} block_prod[j]
 *   Walking backward:  u_{K-1} = invG
 *                      u_b     = u_{b+1} * block_prod[b+1]   for b < K-1
 * Then inside block b, for thread at position p within the block
 * (global index i = base_b + p):
 *   product of z[0..i] = G[b-1] * c_local[i]
 *                      = (product of block_prod[0..b-1]) * c_local[i]
 *   inverse of that    = invG * product_{j>=b} block_prod[j] / c_local[i]
 *                      = u_b * block_prod[b] / c_local[i]
 *                      = u_b * c_local[base_b + B - 1] / c_local[i]
 * But that's product over z[0..i], not z[i] alone. We want z_inv[i] = 1/z[i].
 *
 * Standard Montgomery batch:
 *   Phase forward (already done): c[i] = z[0]*z[1]*...*z[i].
 *   Phase backward: u = c[N-1]^-1; for i = N-1..1:
 *      z_inv[i] = u * c[i-1]
 *      u <- u * z[i]
 *    z_inv[0] = u
 *
 * Block decomposition for parallel suffix:
 *   At block boundary, the running u between blocks corresponds to
 *   the inverse of z[0..base_b-1] ... no wait, u after processing
 *   indices N-1 down to base_b is:
 *      u_in_b = (z[0]*z[1]*...*z[base_b-1])^-1   ??? let's derive.
 *
 *   Initial: u = c[N-1]^-1 = (z[0]*...*z[N-1])^-1.
 *   After step i: u_new = u_old * z[i] = (z[0]*...*z[i-1])^-1.
 *   So after processing i = N-1, N-2, ..., j (j > 0), the next u is
 *     u = (z[0]*z[1]*...*z[j-1])^-1.
 *   In particular, at the moment we "enter" block b (we are about to
 *   process indices base_b+B-1 down to base_b), the u value equals
 *     u_b_entry = (z[0]*...*z[base_b+B-1])^-1 = invG * product_{j > b} block_prod[j]
 *
 *   So: u_b_entry = invG * Π_{j=b+1..K-1} block_prod[j]
 *
 *   Walking backward from b = K-1 to 0:
 *     u_{K-1}_entry = invG     (product over empty set above K-1)
 *     u_b_entry     = u_{b+1}_entry * block_prod[b+1]
 *
 *   Then inside block b in Phase D, each thread can compute its
 *   z_inv[i] **only if** it knows the suffix product within the block
 *   from i+1 to base_b+B-1, plus u_b_entry. Specifically:
 *     z_inv[i] = u_b_entry * (c[i-1]_in_global) * Π_{j=i+1..base_b+B-1} z[j]
 *              ... no wait, simpler form:
 *
 *   Standard recurrence within block, starting with u = u_b_entry,
 *   walking i = base_b+B-1 down to base_b:
 *     z_inv[i] = u * c_global[i-1]   (where c_global[i-1] is the global
 *                                     prefix product z[0..i-1])
 *     u <- u * z[i]
 *
 *   c_global[i-1] = G[b-1] * c_local[i-1 within block b] = (product
 *   of block_prod[0..b-1]) * c_local[i-1]. Storing the per-block
 *   "block prefix to its left", `G_left[b]` = product of block_prod[0..b-1],
 *   simplifies the formula.
 *
 *   But: this requires walking *within* the block serially because z_inv[i]
 *   uses u that depends on z[i+1..base_b+B-1]. So Phase D must also be
 *   serial within a block. That's fine — for now, thread 0 in each block
 *   does the work in Phase D too, and we still get N -> K inversions.
 *
 *   Performance: per-block serial work is bounded; we recover the win
 *   because K = N / B inversions instead of N.
 *
 *   Day-2 simplification: we don't even need G_left. We can store, for
 *   each i, c_global[i] explicitly during Phase B by feeding the
 *   block's starting accumulator from G_left[b]. Phase B writes
 *   c_global[i] = G_left[b] * c_local[i]. Then Phase D's z_inv[i]
 *   uses c_global[i-1] directly. We just need to compute G_left[]
 *   ahead of Phase B.
 *
 *   Actually, splitting it that way needs Phase B to be after Phase C,
 *   which inverts the ordering. Simpler ordering:
 *     Phase A: scalar_mul.
 *     Phase B: per-block prefix of z[] -> c_local[] and block_prod[].
 *     Phase C (serial): walk blocks left-to-right computing G_left[b],
 *       then compute G = G_left[K-1] * block_prod[K-1] = full product,
 *       invert G (single _ModInv), then walk blocks right-to-left
 *       computing u_b_entry. Store G_left[] and u_b_entry[] for D.
 *     Phase D: per-block serial via thread 0 — walk i = end-1 .. base,
 *       computing c_global[i-1] = G_left[b] * c_local[i-1],
 *       z_inv[i] = u * c_global[i-1], u <- u * z[i].
 *       Special case i == base: c_global[i-1] is the product over
 *       blocks 0..b-1 minus none = G_left[b] *only if i-1 < base. We
 *       handle "i == 0" globally by setting z_inv[0] = final u.
 *
 *   This keeps all multiplications correct and only one inversion is
 *   used.
 */

__global__ void v2_phase_c(
    const uint64_t* block_prod,  /* K * 4 */
    uint64_t*       G_left,      /* K * 4 */
    uint64_t*       u_entry,     /* K * 4 */
    uint32_t        K)
{
    if (blockIdx.x != 0 || threadIdx.x != 0) return;
    /* Walk left-to-right computing G_left[b] = Π_{j<b} block_prod[j]. */
    uint64_t acc[4] = {1, 0, 0, 0};
    for (uint32_t b = 0; b < K; b++) {
        G_left[(size_t)b*4+0] = acc[0];
        G_left[(size_t)b*4+1] = acc[1];
        G_left[(size_t)b*4+2] = acc[2];
        G_left[(size_t)b*4+3] = acc[3];
        uint64_t bp[4] = { block_prod[(size_t)b*4+0],
                           block_prod[(size_t)b*4+1],
                           block_prod[(size_t)b*4+2],
                           block_prod[(size_t)b*4+3] };
        uint64_t r[4];
        ModMult4(r, acc, bp);
        acc[0]=r[0]; acc[1]=r[1]; acc[2]=r[2]; acc[3]=r[3];
    }
    /* acc = full product G = Π block_prod[b]. Invert. */
    uint64_t invG[4] = { acc[0], acc[1], acc[2], acc[3] };
    ModInv256(invG);
    /* Walk right-to-left computing u_entry[b] = invG * Π_{j>b} block_prod[j].
     * u_entry[K-1] = invG.
     * u_entry[b]   = u_entry[b+1] * block_prod[b+1].
     */
    if (K == 0) return;
    u_entry[(size_t)(K-1)*4+0] = invG[0];
    u_entry[(size_t)(K-1)*4+1] = invG[1];
    u_entry[(size_t)(K-1)*4+2] = invG[2];
    u_entry[(size_t)(K-1)*4+3] = invG[3];
    uint64_t u[4] = { invG[0], invG[1], invG[2], invG[3] };
    for (int32_t b = (int32_t)K - 2; b >= 0; b--) {
        uint64_t bp[4] = { block_prod[(size_t)(b+1)*4+0],
                           block_prod[(size_t)(b+1)*4+1],
                           block_prod[(size_t)(b+1)*4+2],
                           block_prod[(size_t)(b+1)*4+3] };
        uint64_t r[4];
        ModMult4(r, u, bp);
        u[0]=r[0]; u[1]=r[1]; u[2]=r[2]; u[3]=r[3];
        u_entry[(size_t)b*4+0] = u[0];
        u_entry[(size_t)b*4+1] = u[1];
        u_entry[(size_t)b*4+2] = u[2];
        u_entry[(size_t)b*4+3] = u[3];
    }
}

/* Phase D: per-block serial walk emits z_inv[i] and the final BE bytes.
 * One thread per block does the work.
 *
 * For thread i = base..base+B-1 (within block b), walking i = end-1 .. base:
 *   c_global_prev =
 *      if i == base and b == 0: 1
 *      if i == base and b >  0: G_left[b]  (== Π_{j<b} block_prod[j] = c_global[base-1])
 *                                   wait — c_global[base-1] = G_left[b]
 *                                   * c_local[base-1 of block b-1] ... that's
 *                                   confusing. Let's recompute.
 *
 * c_global[i] = Π z[0..i] = G_left[b] * c_local[i] where c_local[i] is the
 * per-block prefix product written by Phase B (within block b that contains i).
 *
 * For i > base of its block: c_global[i-1] = G_left[b] * c_local[i-1].
 * For i == base of block b > 0: c_global[i-1] = G_left[b] * 1 / 1 ... no:
 *   c_global[base_b - 1] is in the PREVIOUS block b-1 at position end_{b-1} - 1,
 *   which is block_prod[b-1]. And G_left[b] = G_left[b-1] * block_prod[b-1].
 *   So c_global[base_b - 1] = G_left[b-1] * block_prod[b-1] = G_left[b]. Good.
 * For i == base of block 0: c_global[-1] is the empty product = 1.
 *
 * z_inv[i] = u * c_global[i-1], u <- u * z[i].
 * Final after i == base of block 0 is z_inv[0] = u.
 *
 * Each thread (one per block) walks its block backward, leaving its
 * starting u in a per-block scratch slot for cross-block chaining? No —
 * each block uses its own u_entry[b], independent of other blocks.
 *
 * After Phase D finishes for block b, the resulting u value would be
 * u_b_entry * Π z[base_b..base_b+B-1] = u_b_entry * block_prod[b].
 * For b = 0, that equals invG * Π_{j>=1} block_prod[j] * block_prod[0] = 1.
 * So z_inv[0] = invG * Π_{j>=1} block_prod[j] = u_0_entry. Good.
 *
 * So the per-block walk is self-contained given u_entry[b] and G_left[b]
 * and c_local[].
 */

__global__ void v2_phase_d(
    const uint64_t* jac,         /* n * 12 : (Jx, Jy, Jz) */
    const uint64_t* z,           /* n * 4  : effective z = Jz or 1 */
    const uint64_t* c,           /* n * 4  : per-block prefix product */
    const uint64_t* G_left,      /* K * 4  */
    const uint64_t* u_entry,     /* K * 4  */
    uint8_t*        out_points,  /* n * 64 BE */
    uint32_t        n,
    uint32_t        block_size)
{
    uint32_t bid = blockIdx.x;
    if (threadIdx.x != 0) return;
    uint32_t base = bid * block_size;
    if (base >= n) return;
    uint32_t end = base + block_size;
    if (end > n) end = n;

    uint64_t u[4] = { u_entry[(size_t)bid*4+0],
                      u_entry[(size_t)bid*4+1],
                      u_entry[(size_t)bid*4+2],
                      u_entry[(size_t)bid*4+3] };
    uint64_t gl[4] = { G_left[(size_t)bid*4+0],
                       G_left[(size_t)bid*4+1],
                       G_left[(size_t)bid*4+2],
                       G_left[(size_t)bid*4+3] };

    /* Walk i = end-1 down to base. */
    for (int64_t ii = (int64_t)end - 1; ii >= (int64_t)base; ii--) {
        uint32_t i = (uint32_t)ii;
        /* z_inv[i] = u * c_global[i-1] where
         *   c_global[i-1] = gl                  if i == base AND bid == 0 -> 1 sentinel
         *                 = gl                  if i == base AND bid >  0
         *                 = gl * c_local[i-1]   if i >  base
         * Handle bid==0 && i==0 specially by short-circuit (we never read c[-1]).
         */
        uint64_t cgprev[4];
        if (i == base) {
            if (bid == 0) {
                /* c_global[-1] = 1 */
                cgprev[0] = 1; cgprev[1] = 0; cgprev[2] = 0; cgprev[3] = 0;
            } else {
                cgprev[0] = gl[0]; cgprev[1] = gl[1];
                cgprev[2] = gl[2]; cgprev[3] = gl[3];
            }
        } else {
            uint64_t cl[4] = { c[(size_t)(i-1)*4+0], c[(size_t)(i-1)*4+1],
                               c[(size_t)(i-1)*4+2], c[(size_t)(i-1)*4+3] };
            ModMult4(cgprev, gl, cl);
        }
        uint64_t zinv[4];
        ModMult4(zinv, u, cgprev);

        /* u <- u * z[i] */
        uint64_t zi[4] = { z[(size_t)i*4+0], z[(size_t)i*4+1],
                           z[(size_t)i*4+2], z[(size_t)i*4+3] };
        uint64_t newu[4];
        ModMult4(newu, u, zi);
        u[0]=newu[0]; u[1]=newu[1]; u[2]=newu[2]; u[3]=newu[3];

        /* Affine projection. Check infinity: original Jz of slot. */
        const uint64_t* slot = jac + (size_t)i * 12;
        uint64_t Jz[4] = { slot[8], slot[9], slot[10], slot[11] };
        uint8_t* outp = out_points + (size_t)i * 64;
        if (is_zero4(Jz)) {
            for (int j = 0; j < 64; j++) outp[j] = 0;
            continue;
        }
        uint64_t Jx[4] = { slot[0], slot[1], slot[2], slot[3] };
        uint64_t Jy[4] = { slot[4], slot[5], slot[6], slot[7] };
        uint64_t Zi2[4], Zi3[4];
        ModSqr4(Zi2, zinv);
        ModMult4(Zi3, Zi2, zinv);
        uint64_t Ax[4], Ay[4];
        ModMult4(Ax, Jx, Zi2);
        ModMult4(Ay, Jy, Zi3);
        limbs_to_be32(outp, Ax);
        limbs_to_be32(outp + 32, Ay);
    }
}

/* ------------------------------------------------------------------ */
/* Host helpers.                                                      */

static int read_file_to_buffer(const char* path, uint8_t** out_buf, size_t* out_len) {
    FILE* f = fopen(path, "rb");
    if (!f) { fprintf(stderr, "open %s failed\n", path); return -1; }
    fseek(f, 0, SEEK_END);
    long sz = ftell(f);
    fseek(f, 0, SEEK_SET);
    uint8_t* buf = (uint8_t*)malloc(sz);
    if (!buf) { fclose(f); return -1; }
    if (fread(buf, 1, sz, f) != (size_t)sz) { free(buf); fclose(f); return -1; }
    fclose(f);
    *out_buf = buf;
    *out_len = sz;
    return 0;
}

static int write_buffer_to_file(const char* path, const uint8_t* buf, size_t len) {
    FILE* f = fopen(path, "wb");
    if (!f) { fprintf(stderr, "create %s failed\n", path); return -1; }
    if (fwrite(buf, 1, len, f) != len) { fclose(f); return -1; }
    fclose(f);
    return 0;
}

static uint32_t rd32_le(const uint8_t* p) {
    return (uint32_t)p[0] | ((uint32_t)p[1] << 8)
         | ((uint32_t)p[2] << 16) | ((uint32_t)p[3] << 24);
}
static void wr32_le(uint8_t* p, uint32_t v) {
    p[0] = v & 0xff; p[1] = (v >> 8) & 0xff;
    p[2] = (v >> 16) & 0xff; p[3] = (v >> 24) & 0xff;
}

/* Big-endian 32-byte -> 4-limb little-endian u64. */
static void be32_to_limbs(uint64_t out[4], const uint8_t* in) {
    for (int li = 0; li < 4; li++) {
        uint64_t v = 0;
        for (int b = 0; b < 8; b++) {
            v = (v << 8) | in[li*8 + b];
        }
        out[3 - li] = v;
    }
}

/* Emit BSCR error payload to disk so the daemon protocol can write a
 * legal output even on input-validation failure.
 */
static int write_bscr_error(const char* out_path, uint32_t status) {
    uint8_t buf[12];
    memcpy(buf, "BSCR", 4);
    wr32_le(buf + 4, status);
    wr32_le(buf + 8, 0);
    return write_buffer_to_file(out_path, buf, sizeof(buf));
}

/* Global flags controlling kernel variant. Toggled by CLI flags. */
static int g_use_batch_inv = 0;  /* v2 default OFF — Day-2 regressed; Day-3 ladder is canonical */
static int g_window_w      = 0;  /* 0 = v1/v2 binary ladder; 4 or 8 = v3 windowed ladder */

static int process_one_bin(const char* in_path, const char* out_path) {
    uint8_t* in_buf = NULL;
    size_t in_len = 0;
    if (read_file_to_buffer(in_path, &in_buf, &in_len) != 0) return -1;

    if (in_len < 12 || memcmp(in_buf, "BSCP", 4) != 0) {
        fprintf(stderr, "bad BSCP magic / short input (%zu bytes)\n", in_len);
        free(in_buf);
        write_bscr_error(out_path, STATUS_BAD_MAGIC);
        return -1;
    }
    uint32_t op_id = rd32_le(in_buf + 4);
    uint32_t n     = rd32_le(in_buf + 8);

    if (op_id != OP_SCALAR_BASE_MUL) {
        fprintf(stderr, "unsupported op_id 0x%02x\n", op_id);
        free(in_buf);
        write_bscr_error(out_path, STATUS_BAD_OP);
        return -1;
    }

    size_t expected = 12 + 64 + (size_t)n * 32;
    if (in_len != expected) {
        fprintf(stderr, "BSCP size %zu != expected %zu (n=%u)\n", in_len, expected, n);
        free(in_buf);
        write_bscr_error(out_path, STATUS_TRUNCATED);
        return -1;
    }
    if (n == 0) {
        free(in_buf);
        uint8_t hdr[12];
        memcpy(hdr, "BSCR", 4);
        wr32_le(hdr + 4, STATUS_OK);
        wr32_le(hdr + 8, 0);
        return write_buffer_to_file(out_path, hdr, sizeof(hdr));
    }

    const uint8_t* base_x_be  = in_buf + 12;
    const uint8_t* base_y_be  = in_buf + 12 + 32;
    const uint8_t* scalars_be = in_buf + 12 + 64;

    uint64_t base_limbs[8];
    be32_to_limbs(&base_limbs[0], base_x_be);
    be32_to_limbs(&base_limbs[4], base_y_be);

    /* Device buffers */
    uint64_t* d_base = NULL;
    uint8_t*  d_scalars = NULL;
    uint8_t*  d_points  = NULL;
    CUDA_CHECK(cudaMalloc(&d_base,    8 * sizeof(uint64_t)));
    CUDA_CHECK(cudaMalloc(&d_scalars, (size_t)n * 32));
    CUDA_CHECK(cudaMalloc(&d_points,  (size_t)n * 64));

    CUDA_CHECK(cudaMemcpy(d_base, base_limbs, 8 * sizeof(uint64_t),
                          cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_scalars, scalars_be, (size_t)n * 32,
                          cudaMemcpyHostToDevice));

    uint32_t threads_per_block = 128;
    uint32_t blocks = (n + threads_per_block - 1) / threads_per_block;

    /* v2 batch-inv scratch buffers (only allocated if g_use_batch_inv). */
    uint64_t* d_jac        = NULL;
    uint64_t* d_z          = NULL;
    uint64_t* d_c          = NULL;
    uint64_t* d_block_prod = NULL;
    uint64_t* d_G_left     = NULL;
    uint64_t* d_u_entry    = NULL;

    /* Phase B/D block size: each block is handled serially by thread 0.
     * Smaller block_size = more blocks = more thread-0's working in
     * parallel = better GPU utilisation. Larger block_size = fewer
     * cross-block reductions in Phase C. Sweet spot ~256 keys/block. */
    const uint32_t v2_block_size = 256;
    uint32_t v2_K = (n + v2_block_size - 1) / v2_block_size;

    /* v3 windowed-ladder scratch (only allocated if g_window_w > 0). */
    uint64_t* d_table = NULL;
    uint32_t  v3_W    = (g_window_w == 8) ? 256u : 16u;

    struct timespec t0, t1;
    clock_gettime(CLOCK_MONOTONIC, &t0);
    if (g_window_w > 0 && g_use_batch_inv) {
        /* v4: v3 windowed scalar_mul (Phase A) + v2 batch inversion     */
        /* (Phase B/C/D). Stacked Day-4 path.                             */
        CUDA_CHECK(cudaMalloc(&d_table, (size_t)v3_W * 8 * sizeof(uint64_t)));
        CUDA_CHECK(cudaMalloc(&d_jac,        (size_t)n * 12 * sizeof(uint64_t)));
        CUDA_CHECK(cudaMalloc(&d_z,          (size_t)n * 4  * sizeof(uint64_t)));
        CUDA_CHECK(cudaMalloc(&d_c,          (size_t)n * 4  * sizeof(uint64_t)));
        CUDA_CHECK(cudaMalloc(&d_block_prod, (size_t)v2_K * 4 * sizeof(uint64_t)));
        CUDA_CHECK(cudaMalloc(&d_G_left,     (size_t)v2_K * 4 * sizeof(uint64_t)));
        CUDA_CHECK(cudaMalloc(&d_u_entry,    (size_t)v2_K * 4 * sizeof(uint64_t)));

        v3_table_init<<<1, 1>>>(d_base, d_table, v3_W);
        if (g_window_w == 4) {
            v4_phase_a<16><<<blocks, threads_per_block>>>(
                d_base, d_table, d_scalars, d_jac, d_z, n, 4u);
        } else {
            uint32_t tpb8 = 64;
            uint32_t blocks8 = (n + tpb8 - 1) / tpb8;
            v4_phase_a<256><<<blocks8, tpb8>>>(
                d_base, d_table, d_scalars, d_jac, d_z, n, 8u);
        }
        v2_phase_b<<<v2_K, 1>>>(d_z, d_c, d_block_prod, n, v2_block_size);
        v2_phase_c<<<1, 1>>>(d_block_prod, d_G_left, d_u_entry, v2_K);
        v2_phase_d<<<v2_K, 1>>>(d_jac, d_z, d_c, d_G_left, d_u_entry,
                                 d_points, n, v2_block_size);
    } else if (g_window_w > 0) {
        /* v3: precompute (i*P) table, then run windowed-ladder kernel. */
        CUDA_CHECK(cudaMalloc(&d_table, (size_t)v3_W * 8 * sizeof(uint64_t)));
        v3_table_init<<<1, 1>>>(d_base, d_table, v3_W);
        if (g_window_w == 4) {
            secp_mul_batch_kernel_v3<16><<<blocks, threads_per_block>>>(
                d_base, d_table, d_scalars, d_points, n, 4u);
        } else {
            /* w == 8: 256-entry table, 16 KB __shared__ per block.    */
            /* Drop threads/block from 128 to 64 to keep occupancy.    */
            uint32_t tpb8 = 64;
            uint32_t blocks8 = (n + tpb8 - 1) / tpb8;
            secp_mul_batch_kernel_v3<256><<<blocks8, tpb8>>>(
                d_base, d_table, d_scalars, d_points, n, 8u);
        }
    } else if (g_use_batch_inv) {
        CUDA_CHECK(cudaMalloc(&d_jac,        (size_t)n * 12 * sizeof(uint64_t)));
        CUDA_CHECK(cudaMalloc(&d_z,          (size_t)n * 4  * sizeof(uint64_t)));
        CUDA_CHECK(cudaMalloc(&d_c,          (size_t)n * 4  * sizeof(uint64_t)));
        CUDA_CHECK(cudaMalloc(&d_block_prod, (size_t)v2_K * 4 * sizeof(uint64_t)));
        CUDA_CHECK(cudaMalloc(&d_G_left,     (size_t)v2_K * 4 * sizeof(uint64_t)));
        CUDA_CHECK(cudaMalloc(&d_u_entry,    (size_t)v2_K * 4 * sizeof(uint64_t)));

        /* Phase A: full N-thread parallel scalar mul. */
        v2_phase_a<<<blocks, threads_per_block>>>(d_base, d_scalars, d_jac, d_z, n);
        /* Phase B: K blocks, 1 thread each. (Use 1 thread/block for
         * cleanest semantics; could expand to 32-thread warps later.) */
        v2_phase_b<<<v2_K, 1>>>(d_z, d_c, d_block_prod, n, v2_block_size);
        /* Phase C: one thread total, walks blocks for the single inversion. */
        v2_phase_c<<<1, 1>>>(d_block_prod, d_G_left, d_u_entry, v2_K);
        /* Phase D: K blocks, 1 thread each. */
        v2_phase_d<<<v2_K, 1>>>(d_jac, d_z, d_c, d_G_left, d_u_entry,
                                 d_points, n, v2_block_size);
    } else {
        secp_mul_batch_kernel_v1<<<blocks, threads_per_block>>>(
            d_base, d_scalars, d_points, n);
    }
    CUDA_CHECK(cudaDeviceSynchronize());
    clock_gettime(CLOCK_MONOTONIC, &t1);
    double kernel_ms = (t1.tv_sec - t0.tv_sec) * 1e3
                     + (t1.tv_nsec - t0.tv_nsec) / 1e6;

    /* Assemble output */
    size_t out_size = 12 + (size_t)n * 64;
    uint8_t* out_buf = (uint8_t*)malloc(out_size);
    if (!out_buf) {
        cudaFree(d_base); cudaFree(d_scalars); cudaFree(d_points);
        free(in_buf); return -1;
    }
    memcpy(out_buf, "BSCR", 4);
    wr32_le(out_buf + 4, STATUS_OK);
    wr32_le(out_buf + 8, n);
    CUDA_CHECK(cudaMemcpy(out_buf + 12, d_points, (size_t)n * 64,
                          cudaMemcpyDeviceToHost));

    cudaFree(d_base); cudaFree(d_scalars); cudaFree(d_points);
    if (d_jac)        cudaFree(d_jac);
    if (d_z)          cudaFree(d_z);
    if (d_c)          cudaFree(d_c);
    if (d_block_prod) cudaFree(d_block_prod);
    if (d_G_left)     cudaFree(d_G_left);
    if (d_u_entry)    cudaFree(d_u_entry);
    if (d_table)      cudaFree(d_table);
    free(in_buf);

    int wr = write_buffer_to_file(out_path, out_buf, out_size);
    free(out_buf);
    const char* variant_label;
    char vN_label[40];
    if (g_window_w > 0 && g_use_batch_inv) {
        snprintf(vN_label, sizeof(vN_label), "v4-window-w%u+batch-inv", g_window_w);
        variant_label = vN_label;
    } else if (g_window_w > 0) {
        snprintf(vN_label, sizeof(vN_label), "v3-window-w%u", g_window_w);
        variant_label = vN_label;
    } else if (g_use_batch_inv) {
        variant_label = "v2-batch-inv";
    } else {
        variant_label = "v1-per-thread-inv";
    }
    fprintf(stderr, "secp op=0x%02x n=%u kernel=%.2fms %.3f Mkeys/s [%s]\n",
            op_id, n, kernel_ms,
            (kernel_ms > 0) ? ((double)n / 1000.0) / kernel_ms : 0.0,
            variant_label);
    return wr;
}

/* ------------------------------------------------------------------ */

int main(int argc, char** argv) {
    /* Pull --batch-inv / --no-batch-inv out of argv ahead of mode dispatch.
     * Default: ON (Day-2 batch inversion). --no-batch-inv selects v1.
     */
    int argi = 1;
    int new_argc = 1;
    char* new_argv[32];
    new_argv[0] = argv[0];
    while (argi < argc) {
        if (!strcmp(argv[argi], "--batch-inv")) {
            g_use_batch_inv = 1;
            argi++;
        } else if (!strcmp(argv[argi], "--no-batch-inv")) {
            g_use_batch_inv = 0;
            argi++;
        } else if (!strcmp(argv[argi], "--window-w") && argi + 1 < argc) {
            int w = atoi(argv[argi + 1]);
            if (w != 4 && w != 8) {
                fprintf(stderr, "--window-w must be 4 or 8 (got %d)\n", w);
                return 1;
            }
            if (w == 8) {
                fprintf(stderr, "--window-w 8: table init not yet implemented "
                                "(needs host-provided scratch for 256 entries). "
                                "Day-3 ships w=4 only.\n");
                return 1;
            }
            g_window_w = w;
            argi += 2;
        } else if (!strcmp(argv[argi], "--no-window")) {
            g_window_w = 0;
            argi++;
        } else {
            if (new_argc < 31) new_argv[new_argc++] = argv[argi];
            argi++;
        }
    }
    argc = new_argc;
    argv = new_argv;

    if (argc >= 2 && !strcmp(argv[1], "--daemon")) {
        int dev_count = 0;
        cudaGetDeviceCount(&dev_count);
        if (dev_count == 0) {
            fprintf(stderr, "no CUDA devices\n");
            return 1;
        }
        cudaFree(0);
        fprintf(stdout, "ready\n"); fflush(stdout);

        char line[4096];
        while (fgets(line, sizeof(line), stdin)) {
            char* nl = strchr(line, '\n');
            if (nl) *nl = 0;
            if (!strcmp(line, "quit") || !strcmp(line, "bye")) {
                fprintf(stdout, "bye\n"); fflush(stdout);
                break;
            }
            if (!strncmp(line, "process-bin ", 12)) {
                char* args = line + 12;
                char* sep = strchr(args, ' ');
                if (!sep) {
                    fprintf(stdout, "error bad-args\n"); fflush(stdout);
                    continue;
                }
                *sep = 0;
                char* in_path  = args;
                char* out_path = sep + 1;
                int r = process_one_bin(in_path, out_path);
                fprintf(stdout, r == 0 ? "done %s\n" : "error process-failed\n",
                        out_path);
                fflush(stdout);
            } else if (line[0]) {
                fprintf(stdout, "error unknown-command\n"); fflush(stdout);
            }
        }
        return 0;
    }
    if (argc >= 4 && !strcmp(argv[1], "--binary")) {
        return process_one_bin(argv[2], argv[3]);
    }
    fprintf(stderr,
        "usage: %s [variant] --binary <input.bin> <output.bin>\n"
        "       %s [variant] --daemon  (read commands on stdin)\n"
        "\n"
        "  --window-w 4|8   Day-3 v3: windowed-base ladder (precomputed table of i*P)\n"
        "  --batch-inv      Day-2 v2: Montgomery batch inversion across N points\n"
        "  --window-w 4 --batch-inv   Day-4 v4: v3 ladder + v2 batch inversion stacked\n"
        "  --no-batch-inv   Day-1 v1: per-thread _ModInv (default)\n"
        "  --no-window      disables window ladder (default)\n",
        argv[0], argv[0]);
    return 1;
}
