form E day-1 — Montgomery batch-inv 0x0A regressed 300x vs 0x05, parked

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

@@ -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<<<blocks, threads_per_block>>>(d_r_wide, d_a, d_b, n);
             break;
+        case OP_MOD_INV_BATCH:
+            mod_inv_batch_256_kernel<<<batch_inv_blocks, batch_inv_threads>>>(
+                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);
 

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.

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-add",       driver_mod_add,       True),
-    ("mod-sub",    driver_mod_sub,    True),
+    ("mod-sub",       driver_mod_sub,       True),
-    ("mod-mul",    driver_mod_mul,    True),
+    ("mod-mul",       driver_mod_mul,       True),
-    ("mod-sqr",    driver_mod_sqr,    True),
+    ("mod-sqr",       driver_mod_sqr,       True),
-    ("mod-inv",    driver_mod_inv,    True),
+    ("mod-inv",       driver_mod_inv,       True),
-    ("mod-exp",    driver_mod_exp,    True),
+    ("mod-exp",       driver_mod_exp,       True),
-    ("mod-reduce", driver_mod_reduce, True),
+    ("mod-reduce",    driver_mod_reduce,    True),
-    ("add-no-mod", driver_add_no_mod, False),
+    ("add-no-mod",    driver_add_no_mod,    False),
-    ("mul-no-mod", driver_mul_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"),
     ]