bend form A — Day-3 v3 windowed-base ladder, 1.73x at n=1M on 3090
Per Day-2 progress doc the kernel bottleneck was scalar_mul (256 doubles +
~128 adds per scalar), not _ModInv. v3 swaps the binary double-and-add for
a windowed-base ladder: precompute table[0..15] = i*P affine on-device,
walk scalar 4 bits at a time MSB->LSB, cutting per-scalar adds from ~128
to ~63. Table loaded into __shared__ (1024 B) per block.
Benched on 3090-ai (best-of-3, --no-batch-inv):
n=10k: v1 4.22 Mkeys/s v3 2.18 Mkeys/s 0.52x (init overhead dominates)
n=100k: v1 7.48 Mkeys/s v3 9.29 Mkeys/s 1.24x
n=1M: v1 7.87 Mkeys/s v3 13.60 Mkeys/s 1.73x
Byte-identical against coincurve at n in {32, 1000, 10000, 100000, 1000000}.
CLI: --window-w 4 selects v3 (Day-3 canonical). w=8 reserved but stub-
rejected since device-side table init is register-stack-bounded at W <= 16.
Default no-flag behaviour stays v1 (Day-1) so the gpu-worker.lsp daemon
inherits the safe baseline until fox routes traffic to v3.
Companion progress doc at plans/form-A-day3-progress.md.
This commit is contained in:
parent
b2209ac337
commit
ecfe27a46f
2 changed files with 837 additions and 15 deletions
83
examples/cuda-fanout/plans/form-A-day3-progress.md
Normal file
83
examples/cuda-fanout/plans/form-A-day3-progress.md
Normal file
|
|
@ -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<W>` (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<W>`, 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.
|
||||
|
|
@ -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 <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 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);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -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<<<blocks, threads_per_block>>>(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><<<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
|
||||
|
|
@ -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 <input.bin> <output.bin>\n"
|
||||
" %s --daemon (read commands on stdin)\n",
|
||||
"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"
|
||||
" --no-batch-inv Day-1 v1: per-thread _ModInv (default)\n"
|
||||
" --no-window disables window ladder (default)\n",
|
||||
argv[0], argv[0]);
|
||||
return 1;
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue