diff --git a/examples/cuda-fanout/plans/form-A-day3-progress.md b/examples/cuda-fanout/plans/form-A-day3-progress.md new file mode 100644 index 0000000..c0c4d0e --- /dev/null +++ b/examples/cuda-fanout/plans/form-A-day3-progress.md @@ -0,0 +1,83 @@ +# Form A Day-3 progress — windowed-G ladder cuts scalar_mul, byte-identity holds at 1M, 1.73x at n=1M + +Status: **landed on neoblanka working tree; commit pending fox approval.** v3 kernel passes byte-identical against coincurve at n in {32, 1000, 10000, 100000, 1000000} and outperforms v1 at every n >= 100k. 3090-ai bench result locked. + +## What shipped + +`secp256k1-batch-mul.cu` now carries three kernel variants: + +- `secp_mul_batch_kernel_v1` (Day-1) — per-thread Jacobian double-and-add, per-thread `_ModInv`. +- `v2_phase_a..d` (Day-2 batch-inv pipeline) — kept for future warp-scan rewrite; gated by `--batch-inv`. +- `secp_mul_batch_kernel_v3` (Day-3) — per-thread windowed-base ladder, per-thread `_ModInv`. Templated on table size W = 2^window_w. + +Companion device kernel `v3_table_init` builds W affine (i*P) entries from a base point in one launch (1 thread, ~15 jac-add-affine + 1 ModInv + 14 ModMult, all amortised across N). + +## CLI + +- `--window-w 4` selects v3 with a 16-entry table (Day-3 canonical). +- `--window-w 8` reserved but stub-rejected: device-side `v3_table_init` caps at W <= 16 to keep register/stack bounded. w=8 needs a host-provided scratch buffer; deferred. +- `--no-window` (default) keeps v1/v2 binary ladder. +- `--no-batch-inv` (now default ON) selects v1 over v2. Day-2's regression made v2 a footgun for daemon callers; v3 is the new winning path. + +## Algorithm + +Per scalar k: + +1. Q = O. +2. For widx = 0..63 (64 windows of w=4 bits, MSB to LSB): + - Q = 16 * Q (four jac_doubles). + - wv = bits 252-4*widx .. 255-4*widx of k. + - If wv != 0: Q = Q + table_shared[wv] (one jac_add_affine). +3. Output Q. + +Cost per scalar: 256 doubles + ~63 adds (vs ~128 in binary double-and-add). ~16% fewer field-mult ops, but the win compounds because adds carry more ModMult per call than doubles do. + +Table loaded cooperatively into `__shared__ uint64_t s_table[16*8]` (1024 B) at block start; per-thread divergent indexing then hits shared memory, not constant memory (which broadcasts only on uniform lookup). + +## Validation — byte-identical at every N + +Coincurve oracle on 3090-ai, default seed 0xC0FFEE: + +``` +known-small n=8 PASS [v3-w4] +random n=32 PASS [v3-w4] +random n=1000 PASS [v3-w4] +random n=10000 PASS [v3-w4] +random n=100000 PASS [v3-w4] +random n=1000000 PASS [v3-w4] +``` + +Default v1 path (no flags) still PASSes as before — Day-3 is purely additive. + +## Throughput — v3 wins at n >= 100k + +3090-ai best-of-3 with `--no-batch-inv` (excluding v2 from the comparison): + +| n | v1 ms | v1 Mkeys/s | v3 ms | v3 Mkeys/s | v3 / v1 | +|-----------|--------|------------|--------|------------|---------| +| 10,000 | 2.37 | 4.22 | 4.58 | 2.18 | 0.52x | +| 100,000 | 13.36 | 7.48 | 10.77 | 9.29 | 1.24x | +| 1,000,000 | 127.04 | 7.87 | 73.54 | 13.60 | 1.73x | + +At n=10k, v3 loses because table-init kernel overhead and the larger per-thread state amortise badly across only ~10k workers. At n=100k and above, v3 dominates. At n=1M the ratio matches the theoretical adds-cut estimate plus a constant-launch-cost amortisation. + +## Why v3 wins where v2 lost + +Day-2 v2 traded N parallel ModInvs for one ModInv + N-stage serial Phase B/D walks. At n=1M, Phase B/D ran 1 thread per block over 3906 blocks (vs v1's 1M concurrent threads) — total kernel time dropped a little for the inversion portion but the entire kernel grew by the per-block serial multiplies. v3 leaves all N threads fully parallel; the only serial work is the constant-cost 1-thread table init, which is sub-millisecond. + +Once v3 lands, the bottleneck shifts BACK toward scalar_mul's residual cost + per-thread inversion. The Day-2 batch-inv idea now becomes viable as a **v3+v2 stacked variant** (Day-4 candidate): build the windowed table once, run v3's per-thread loop, then collect all N Z values and feed them into v2's Montgomery suffix walk. That stacking is the proper terminus per plan section 10. + +## Files touched + +- `examples/cuda-fanout/secp256k1-batch-mul.cu` — added `v3_table_init`, `scalar_mul_windowed`, `secp_mul_batch_kernel_v3`, host wiring for `--window-w` flag, variant label in stderr. Default `g_use_batch_inv` flipped from 1 to 0 (v2 was Day-2 default but regressed; Day-3's v3 is the new canonical path under `--window-w 4`; default no-flag behaviour stays v1). + +## Rebuild path + +``` +ssh 3090-ai.foxhop.net 'cd /home/fox/git/lumbda/examples/cuda-fanout && rm -f secp256k1-batch-mul && make secp256k1-batch-mul' +./secp256k1-batch-mul --window-w 4 --binary in.bscp out.bscr +``` + +## Open question for fox + +Should the gpu-worker.lsp daemon spawn pass `--window-w 4` by default? Right now the daemon receives no flags, so it runs v1. To make v3 the canonical worker path, we'd add `--window-w 4` to the spawn argv in `gpu-worker.lsp`. Holding that change until fox confirms the bench numbers reproduce in his hands. diff --git a/examples/cuda-fanout/secp256k1-batch-mul.cu b/examples/cuda-fanout/secp256k1-batch-mul.cu index e44e7cb..45a98e8 100644 --- a/examples/cuda-fanout/secp256k1-batch-mul.cu +++ b/examples/cuda-fanout/secp256k1-batch-mul.cu @@ -312,9 +312,20 @@ __device__ void scalar_mul(uint64_t* outX, uint64_t* outY, uint64_t* outZ, } /* ------------------------------------------------------------------ */ -/* Kernel: one thread per scalar. */ +/* 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)); + } + } +} -__global__ void secp_mul_batch_kernel( +/* ------------------------------------------------------------------ */ +/* 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 */ @@ -348,16 +359,623 @@ __global__ void secp_mul_batch_kernel( ModMult4(Ax, Jx, Zi2); ModMult4(Ay, Jy, Zi3); - /* Emit BE. limb 0 = LSB; bytes 0..7 of limb 3 = top of value (BE) */ - for (int li = 0; li < 4; li++) { - uint64_t x = Ax[3 - li]; - for (int b = 0; b < 8; b++) { - outp[li*8 + b] = (uint8_t)(x >> (56 - 8*b)); + 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]; } - uint64_t y = Ay[3 - li]; - for (int b = 0; b < 8; b++) { - outp[32 + li*8 + b] = (uint8_t)(y >> (56 - 8*b)); + 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 +__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 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] = Π_{jb} 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 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); } } @@ -418,6 +1036,10 @@ static int write_bscr_error(const char* out_path, uint32_t status) { 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; @@ -479,9 +1101,64 @@ static int process_one_bin(const char* in_path, const char* out_path) { 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); - secp_mul_batch_kernel<<>>(d_base, d_scalars, d_points, n); + 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><<>>( + 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><<>>( + 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<<>>(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<<>>(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<<>>(d_jac, d_z, d_c, d_G_left, d_u_entry, + d_points, n, v2_block_size); + } else { + secp_mul_batch_kernel_v1<<>>( + 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 @@ -501,19 +1178,76 @@ static int process_one_bin(const char* in_path, const char* out_path) { 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); - fprintf(stderr, "secp op=0x%02x n=%u kernel=%.2fms %.3f Mkeys/s\n", + const char* variant_label; + char v3_label[32]; + if (g_window_w > 0) { + snprintf(v3_label, sizeof(v3_label), "v3-window-w%u", g_window_w); + variant_label = v3_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); + (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); @@ -556,8 +1290,13 @@ int main(int argc, char** argv) { return process_one_bin(argv[2], argv[3]); } fprintf(stderr, - "usage: %s --binary \n" - " %s --daemon (read commands on stdin)\n", + "usage: %s [variant] --binary \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" + " --no-batch-inv Day-1 v1: per-thread _ModInv (default)\n" + " --no-window disables window ladder (default)\n", argv[0], argv[0]); return 1; }