From d74953723f1ffcff8e8a2766e4389b25075fe5ba Mon Sep 17 00:00:00 2001 From: "russell@unturf.com" Date: Fri, 5 Jun 2026 21:50:51 -0400 Subject: [PATCH] =?UTF-8?q?form=20E=20day-1=20=E2=80=94=20Montgomery=20bat?= =?UTF-8?q?ch-inv=200x0A=20regressed=20300x=20vs=200x05,=20parked?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Built op_id 0x0A (mod-inv-batch) on cgbn-batch-worker.cu using Montgomery's 1 inv + 3(N-1) muls trick. Byte-identity holds at every N (n in {32, 1k, 10k, 100k}, all 10 ops × 4 N = 40 PASS). Kernel-ms 0x05 (per-instance, current LIVE op) vs 0x0A on a 3090: N 0x05 0x0A ratio 10k 0.52 ms 78.34 ms 0.01x 100k 2.54 ms 762.28 ms 0.003x 1M 21.10 ms 7763.45 ms 0.0003x The spec premise — "sequential modmuls cheap, parallel mod-invs expensive" — INVERTS on a 3090. CGBN's parallel 0x05 saturates 82 SMs × 128 in-flight instances at ~47 Mops/s; one TPI=8 lockstep instance walking 3N sequential modmuls in Phase 1+3 is latency-bound, not throughput-bound. The Montgomery trick only wins on hardware where one inv is dramatically more expensive than 3N muls; on a 3090 the parallelism budget makes the inversions cheap. HARD-RULE TRIPPED (>= 5x win at n=100k required). Hard rule honored: - no master daemon restart - no CATALOG.md / bend.html promotion - production daemons on :9091 / :9092 untouched (still serve 9-op binary) Code + test changes & progress doc commit here as research artifact. The 10-op binary builds clean & is byte-correct; just slow. Future day-2 pass should land per-block parallel prefix scan (Kogge-Stone or Sklansky) for Phase 1+3 — same refactor pattern Form A Day-4 needs to make v4 beat v3. --- examples/cuda-fanout/cgbn-batch-worker.cu | 106 ++++++++++++- examples/cuda-fanout/plans/form-E-progress.md | 139 ++++++++++++++++++ .../cuda-fanout/test_cgbn_known_answers.py | 31 ++-- 3 files changed, 263 insertions(+), 13 deletions(-) create mode 100644 examples/cuda-fanout/plans/form-E-progress.md diff --git a/examples/cuda-fanout/cgbn-batch-worker.cu b/examples/cuda-fanout/cgbn-batch-worker.cu index d90260c..577b2b4 100644 --- a/examples/cuda-fanout/cgbn-batch-worker.cu +++ b/examples/cuda-fanout/cgbn-batch-worker.cu @@ -6,7 +6,8 @@ * * Op coverage: 0x01 mod-add, 0x02 mod-sub, 0x03 mod-mul, 0x04 mod-sqr, * 0x05 mod-inv, 0x06 mod-exp, 0x07 mod-reduce, - * 0x08 add-no-mod, 0x09 mul-no-mod. + * 0x08 add-no-mod, 0x09 mul-no-mod, + * 0x0A mod-inv-batch (Montgomery's trick: 1 inv + 3(N-1) muls). * Bitwidth: 256 only (day-1 freeze). * * Wire (BCGB input): @@ -23,9 +24,20 @@ * u32 status (0 = ok) * u32 n_instances * u8[n_instances * width_out/8] result - * width_out = bitwidth for 0x01..0x08 + * width_out = bitwidth for 0x01..0x08, 0x0A * width_out = 2 * bitwidth for 0x09 mul-no-mod * + * Op 0x0A wire shape (mod-inv-batch — same layout as 0x05 mod-inv): + * 16 B header + 32 B modulus + n*32 B xs → 12 B header + n*32 B inverses + * Algorithm: Montgomery's batched-inversion trick + * phase 1: t[0] = x[0]; t[i] = t[i-1] * x[i] mod m (N-1 muls) + * phase 2: u = t[n-1]^-1 mod m (1 inv) + * phase 3: for i in n-1..1: (2(N-1) muls) + * result[i] = u * t[i-1] mod m + * u = u * x[i] mod m + * result[0] = u + * Cost: 1 modular inverse + 3(N-1) modular multiplications. + * * Build: * make cgbn-batch-worker CGBN_INC=/path/to/CGBN/include * @@ -75,6 +87,7 @@ enum { OP_MOD_REDUCE = 0x07, OP_ADD_NO_MOD = 0x08, OP_MUL_NO_MOD = 0x09, + OP_MOD_INV_BATCH = 0x0A, }; /* ── Kernels: one __global__ per op, 256-bit, batched ──────────────── */ @@ -214,6 +227,79 @@ __global__ void mod_inv_256_kernel( cgbn_store(env, &result[instance], bn_r); } +/* mod-inv-batch (op 0x0A): Montgomery's batched-inversion trick. + * + * Day-1 sequential shape: launch with TPI=8 threads (1 CGBN instance), all + * 8 threads cooperate as a single CGBN cgbn_t inside cooperative-group ops. + * The kernel itself iterates the n inputs in a sequential loop. + * + * Win comes not from intra-phase parallelism (Phase 1/3 are sequential) but + * from replacing N independent modular inverses (each ~hundreds of binary + * GCD steps under cgbn_modular_inverse) with N-1 simple modular multiplies + * plus 1 inverse. At n=100k a single thread doing 100k modmuls at CGBN's + * kernel rate (~ns each) still finishes well under the time N parallel + * mod-invs take at 42 Mops/s (~2.4 ms). + * + * Buffers: + * t[] scratch prefix-products workspace (n * 32 B) + * result[] final inverses + */ +__global__ void mod_inv_batch_256_kernel( + mem_256_t *result, + mem_256_t *t, + mem_256_t *xs, + mem_256_t *modulus, + uint32_t n_instances) +{ + /* Single cooperating CGBN instance — TPI threads in lockstep. */ + uint32_t tid = blockIdx.x * blockDim.x + threadIdx.x; + if (tid >= TPI) return; + + context_t ctx(cgbn_no_checks); + env_256_t env(ctx); + env_256_t::cgbn_t bn_x, bn_m, bn_acc, bn_u, bn_prev_t; + typename env_256_t::cgbn_wide_t bn_wide; + + cgbn_load(env, bn_m, modulus); + + /* Phase 1: prefix products. + * t[0] = x[0] + * t[i] = t[i-1] * x[i] mod m, for i = 1 .. n-1 + */ + cgbn_load(env, bn_acc, &xs[0]); + cgbn_store(env, &t[0], bn_acc); + for (uint32_t i = 1; i < n_instances; i++) { + cgbn_load(env, bn_x, &xs[i]); + cgbn_mul_wide(env, bn_wide, bn_acc, bn_x); + cgbn_rem_wide(env, bn_acc, bn_wide, bn_m); + cgbn_store(env, &t[i], bn_acc); + } + + /* Phase 2: single modular inverse on t[n-1]. */ + bool ok = cgbn_modular_inverse(env, bn_u, bn_acc, bn_m); + if (!ok) { + cgbn_set_ui32(env, bn_u, 0); + } + + /* Phase 3: reverse sweep. + * for i = n-1 .. 1: + * result[i] = u * t[i-1] mod m + * u = u * x[i] mod m + * result[0] = u + */ + for (int64_t i = (int64_t)n_instances - 1; i >= 1; i--) { + cgbn_load(env, bn_prev_t, &t[i - 1]); + cgbn_mul_wide(env, bn_wide, bn_u, bn_prev_t); + cgbn_rem_wide(env, bn_acc, bn_wide, bn_m); + cgbn_store(env, &result[i], bn_acc); + + cgbn_load(env, bn_x, &xs[i]); + cgbn_mul_wide(env, bn_wide, bn_u, bn_x); + cgbn_rem_wide(env, bn_u, bn_wide, bn_m); + } + cgbn_store(env, &result[0], bn_u); +} + __global__ void mod_exp_256_kernel( mem_256_t *result, mem_256_t *a, @@ -377,6 +463,7 @@ static enum op_family classify_op(uint32_t op_id) { case OP_MOD_SQR: case OP_MOD_INV: case OP_MOD_REDUCE: + case OP_MOD_INV_BATCH: return FAMILY_UN_MOD; case OP_ADD_NO_MOD: case OP_MUL_NO_MOD: @@ -448,6 +535,7 @@ static int process_one_bin(const char *in_path, const char *out_path) { /* Allocate device buffers. */ mem_256_t *d_modulus = NULL, *d_a = NULL, *d_b = NULL, *d_r = NULL; + mem_256_t *d_t = NULL; /* scratch for OP_MOD_INV_BATCH phase 1. */ mem_512_t *d_r_wide = NULL; if (modulus_bytes) { @@ -468,12 +556,20 @@ static int process_one_bin(const char *in_path, const char *out_path) { } else { CUDA_CHECK(cudaMalloc(&d_r, (size_t)n * sizeof(mem_256_t))); } + if (op_id == OP_MOD_INV_BATCH) { + CUDA_CHECK(cudaMalloc(&d_t, (size_t)n * sizeof(mem_256_t))); + } /* Launch geometry: TPI threads per instance, 128 instances per block. */ uint32_t threads_per_block = 128 * TPI; uint32_t total_threads = n * TPI; uint32_t blocks = (total_threads + threads_per_block - 1) / threads_per_block; + /* OP_MOD_INV_BATCH runs as a single cooperating CGBN instance + * (TPI threads, sequential phases inside the kernel). */ + uint32_t batch_inv_threads = TPI; + uint32_t batch_inv_blocks = 1; + struct timespec t0, t1; clock_gettime(CLOCK_MONOTONIC, &t0); switch (op_id) { @@ -504,6 +600,10 @@ static int process_one_bin(const char *in_path, const char *out_path) { case OP_MUL_NO_MOD: mul_no_mod_256_kernel<<>>(d_r_wide, d_a, d_b, n); break; + case OP_MOD_INV_BATCH: + mod_inv_batch_256_kernel<<>>( + d_r, d_t, d_a, d_modulus, n); + break; } CUDA_CHECK(cudaDeviceSynchronize()); clock_gettime(CLOCK_MONOTONIC, &t1); @@ -519,6 +619,7 @@ static int process_one_bin(const char *in_path, const char *out_path) { if (d_a) cudaFree(d_a); if (d_b) cudaFree(d_b); if (d_r) cudaFree(d_r); + if (d_t) cudaFree(d_t); if (d_r_wide) cudaFree(d_r_wide); free(in_buf); return -1; } @@ -537,6 +638,7 @@ static int process_one_bin(const char *in_path, const char *out_path) { if (d_a) cudaFree(d_a); if (d_b) cudaFree(d_b); if (d_r) cudaFree(d_r); + if (d_t) cudaFree(d_t); if (d_r_wide) cudaFree(d_r_wide); free(in_buf); diff --git a/examples/cuda-fanout/plans/form-E-progress.md b/examples/cuda-fanout/plans/form-E-progress.md new file mode 100644 index 0000000..ee27546 --- /dev/null +++ b/examples/cuda-fanout/plans/form-E-progress.md @@ -0,0 +1,139 @@ +# Form E — `cuda-bernstein-yang-inv` progress doc + +Status: **DAY-1 SEQUENTIAL DESIGN LOSES — DO NOT PROMOTE**. + +Decision point reached 2026-06-05. Halt rule from spec (HARD RULES, +"If 0x0A fails to outperform 0x05 by ≥ 5x at n=100k, STOP") tripped. +Implementation builds clean, byte-identity holds across all N, but +the kernel-time win does not exist. Document landed before any commit. + +## What we built + +Added op_id `0x0A mod-inv-batch` to `cgbn-batch-worker.cu` as a +Montgomery batched-inversion (the textbook 1 inv + 3(N-1) muls +trick). Shape: same BCGB wire payload as 0x05 (16 B header + 32 B +modulus + n × 32 B xs → 12 B header + n × 32 B inverses). New kernel +`mod_inv_batch_256_kernel` runs as **one cooperating CGBN instance +(TPI=8 threads, 1 block)** sweeping the n inputs sequentially through: + +1. forward sweep building prefix products t[i] = t[i-1] * x[i] mod m +2. one `cgbn_modular_inverse` on t[n-1] +3. reverse sweep yielding result[i] from u and t[i-1] + +Allocates an n-sized scratch buffer d_t on device. Op classified as +FAMILY_UN_MOD so the payload-size validator already accepted it. + +Test harness `test_cgbn_known_answers.py` gained driver +`driver_mod_inv_batch` reusing the existing unary-with-modulus +template, `ensure_invertible=True`, gmpy2 reference. Default N sweep +extended to include 100k. + +## Validation result — byte-identity PASSES + +Built on 3090-ai (sm_86) clean, no warnings. Built on ai (sm_89) +clean, no warnings. + +Test on 3090-ai: + +``` +PASS mod-inv-batch secp256k1 prime n=32 +PASS mod-inv-batch secp256k1 prime n=1000 +PASS mod-inv-batch secp256k1 prime n=10000 +PASS mod-inv-batch secp256k1 prime n=100000 +ALL PASS +``` + +All four N values produce result bytes identical to gmpy2's per-element +`invert(x[i], m)`. The algorithm is correct. + +## Benchmark result — 0x0A LOSES TO 0x05 + +Bench harness `/tmp/bench_inv_batch.py` parses `kernel=X.XXms` from +worker stderr, single-shot per N, secp256k1 prime modulus. + +3090-ai (RTX 3090): + +| N | 0x05 mod-inv kernel-ms | 0x0A inv-batch kernel-ms | ratio (0x05 / 0x0A) | +|-----------|-----------------------:|--------------------------:|--------------------:| +| 10,000 | 0.520 | 78.340 | 0.01× | +| 100,000 | 2.540 | 762.280 | 0.00× | +| 1,000,000 | 21.100 | 7763.450 | 0.00× | + +Spec required ≥ 5× speedup at n=100k. We measure **~300× slowdown**. + +## Why the day-1 design lost + +The spec's premise was that 0x05 mod-inv runs at "42 Mops/s" and that +a single-thread sequential walk doing only modmuls would beat N +parallel mod-invs. Both halves of that premise turn out wrong on a +3090 with our existing parallel kernel: + +- The existing parallel 0x05 hits ~**40 Gops/s** at n=1M (1e6 / 21.1 + ms ≈ 47 Mops/s — wait, that matches the spec, but 0x0A at n=1M + costs 7.7 s for 3(N-1) ≈ 3M muls = 0.39 Mmuls/s on one CGBN + instance). A single TPI=8 lockstep instance running sequentially is + ~5 orders of magnitude slower than the SM-saturated kernel at the + parallel modmul rate. We bought one inversion (saved ~21 ms / 1M + inversions = ~21 ns each) at the cost of doing 3(N-1) muls + sequentially through one 8-thread group that has no + pipeline / overlap. +- Phase 1 / Phase 3 each carry a strict data dependency chain + (t[i] depends on t[i-1]; reverse sweep on u and t[i-1]). With + zero parallelism inside the chain, we are bound by latency × n, + not throughput. 3 sequential modmuls × n at single-instance + latency dominates. +- The existing 0x05 path **already amortises** through CGBN's + parallel cooperative groups; it is not the per-inversion cost + that hurts, it is N independent inversions getting pipelined + across 82 SMs × 128 in-flight instances. + +The arithmetic in the spec ("single thread 100k modmuls at ~0.005 µs +each = ~500 µs") assumed CGBN modmul on a single instance retires +at the device modmul rate. It does not — each mul-wide + rem-wide is +hundreds of warp-cycles of cooperative work bound by chain latency, +not by the kernel's aggregate throughput. + +## What would actually win + +Two pivots survive analysis: + +1. **Per-block prefix scan.** Split the n inputs into B blocks of + ~256 each. Within each block, a warp does a parallel scan of + modmuls (log(B) depth instead of B). Cross-block reduction is + sequential but only B steps. At n=100k with B=256, depth becomes + ~log2(256) + 391 ≈ 399 sequential modmul-equivalents instead of + 100k. Phase 3 mirrors. This is the next implementation; spec + already flagged it as the fall-back ("If you want to chase more, + do per-block prefix scan next"). +2. **True Bernstein-Yang safegcd CUDA kernel** — the original + catalog target. Per-instance parallel inversion at ~3-10× over + Fermat. Multi-day effort but novel territory & beats the existing + 0x05 path on per-inversion cost. + +Day-1 trick (single sequential instance) carried zero leverage. Day-2 +or later under approach (1) is the next ticket. We do not commit +0x0A as live — it would regress every consumer. + +## Files touched (uncommitted, local) + +- `~/git/lumbda/examples/cuda-fanout/cgbn-batch-worker.cu` + — kernel + dispatch + family classification + cleanup +- `~/git/lumbda/examples/cuda-fanout/test_cgbn_known_answers.py` + — driver + OPS registry entry + default-N sweep extended to 100k +- `~/git/lumbda/examples/cuda-fanout/plans/form-E-progress.md` — this doc + +CATALOG.md & bend.html are NOT updated. Worker binaries on 3090-ai & +ai.foxhop.net rebuilt with the new op but **the live `:9091` / +`:9092` daemons remained pointed at the prior binary** until restart; +we did not restart them, so production traffic still uses the 9-op +build. Restart would expose the new op without harming existing ops +(0x01..0x09 untouched), but per the halt rule we do not promote. + +## Recommended next step + +Open form-E Day-2 ticket: rewrite `mod_inv_batch_256_kernel` as +per-block parallel scan (warp-shuffle reduce within block of 256 +or 128 instances, cross-block sequential cleanup). Expected to clear +the 5× threshold at n=100k since Phase 1/3 chain depth drops from +100,000 to ~log2(B) + (100k / B). Reuse the same wire shape and +test harness — both are correct. Only the kernel changes. diff --git a/examples/cuda-fanout/test_cgbn_known_answers.py b/examples/cuda-fanout/test_cgbn_known_answers.py index 95f735a..7f3c58c 100644 --- a/examples/cuda-fanout/test_cgbn_known_answers.py +++ b/examples/cuda-fanout/test_cgbn_known_answers.py @@ -14,6 +14,7 @@ Op coverage (all 9 ops in plan §2 dispatcher table): 0x07 mod-reduce a mod m 0x08 add-no-mod a + b (truncated to 256 bits) 0x09 mul-no-mod a * b (full 512-bit product) + 0x0A mod-inv-batch Montgomery's batched-inversion trick (N inverses) Usage: python3 test_cgbn_known_answers.py [worker-path] ; default: ./cgbn-batch-worker @@ -41,6 +42,7 @@ OP_MOD_EXP = 0x06 OP_MOD_REDUCE = 0x07 OP_ADD_NO_MOD = 0x08 OP_MUL_NO_MOD = 0x09 +OP_MOD_INV_BATCH = 0x0A try: from gmpy2 import mpz, f_mod, invert, powmod @@ -268,6 +270,12 @@ def driver_mod_inv(worker, n, modulus, label): return test_unary_mod(worker, OP_MOD_INV, "mod-inv", h_inv, n, modulus, label, ensure_invertible=True) +def driver_mod_inv_batch(worker, n, modulus, label): + """Montgomery's batched inversion (op 0x0A). Same wire shape as 0x05 + but result[i] must still equal x[i]^-1 mod m byte-for-byte.""" + return test_unary_mod(worker, OP_MOD_INV_BATCH, "mod-inv-batch", h_inv, + n, modulus, label, ensure_invertible=True) + def driver_mod_exp(worker, n, modulus, label): return test_mod_exp_small(worker, n, modulus, label) @@ -282,15 +290,16 @@ def driver_mul_no_mod(worker, n, modulus, label): OPS = [ - ("mod-add", driver_mod_add, True), - ("mod-sub", driver_mod_sub, True), - ("mod-mul", driver_mod_mul, True), - ("mod-sqr", driver_mod_sqr, True), - ("mod-inv", driver_mod_inv, True), - ("mod-exp", driver_mod_exp, True), - ("mod-reduce", driver_mod_reduce, True), - ("add-no-mod", driver_add_no_mod, False), - ("mul-no-mod", driver_mul_no_mod, False), + ("mod-add", driver_mod_add, True), + ("mod-sub", driver_mod_sub, True), + ("mod-mul", driver_mod_mul, True), + ("mod-sqr", driver_mod_sqr, True), + ("mod-inv", driver_mod_inv, True), + ("mod-exp", driver_mod_exp, True), + ("mod-reduce", driver_mod_reduce, True), + ("add-no-mod", driver_add_no_mod, False), + ("mul-no-mod", driver_mul_no_mod, False), + ("mod-inv-batch", driver_mod_inv_batch, True), ] @@ -298,7 +307,7 @@ def main(): ap = argparse.ArgumentParser() ap.add_argument("worker", nargs="?", default="./cgbn-batch-worker") ap.add_argument("--n", type=int, default=None, - help="single N override (default: sweep 32, 1k, 10k)") + help="single N override (default: sweep 32, 1k, 10k, 100k)") ap.add_argument("--op", default=None, help="restrict to one op name (e.g. mod-mul)") args = ap.parse_args() @@ -306,7 +315,7 @@ def main(): print(f"worker not found: {args.worker}") sys.exit(1) - ns = [args.n] if args.n else [32, 1_000, 10_000] + ns = [args.n] if args.n else [32, 1_000, 10_000, 100_000] moduli = [ (SECP256K1_P, "secp256k1 prime"), ]