engine-cuda: batches of whole grids, early reject on element 0, nonce direction

Three scheduling changes measured against the #23 kernel, three interleaved
rounds each, parity clean on every host:

- Batches are rounded up to whole grids, two per batch by default
  (MINER_CUDA_BATCH_ALIGN, MINER_CUDA_BATCH_WAVES). The CLI's 1M-nonce
  batch left a half-empty tail wave and paid the launch gap every 1M
  hashes. beast +1.6% for one grid, +0.4% more for two; benjy +0.5% for
  the second grid; quadbrat neutral.
- LAIR_EARLY_REJECT: the second permutation stops before its last linear
  layer, computes output 0 alone (the hash's top 64 bits) and rejects on
  it; only the rare survivors pay for the other eleven outputs. +0.6%.
- LAIR_NONCE_DIR: the host precomputes the state after the first linear
  layer for the batch's first nonce (first_layer_after_absorb, unit test
  against the matrix column), and the kernel derives each nonce's state
  as that plus a scalar times column 7 of the external matrix; nonces
  past a carry out of the low limb take the old path. +0.25% on beast,
  +0.6% on the 4090 and 3060.

Also tried and rejected: 512-thread blocks (-0.4%), a smaller grid via
MINER_CUDA_THREADS_PER_SM=16384 (-0.6%), a larger one (65536, within
noise of two waves), and the whole S-box as one PTX block, which
compiles to byte-identical SASS. The threads-per-block knob
(MINER_CUDA_THREADS_PER_BLOCK with LAIR_TPB) stays for experiments.

beast 2x5090: 2146 -> ~2206 MH/s (+2.8%); benjy 4090: 711 -> ~718;
quadbrat 3060: 140.2 -> 141.2.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Ue5ZZm1Hiv5zPnucykKKuF

crates/engine-cuda/Cargo.toml

@@ -12,6 +12,7 @@ metrics = { path = "../metrics" }
 primitive-types = { workspace = true }
 log = { workspace = true }
 cudarc = { version = "0.19", default-features = false, features = ["driver", "nvrtc", "dynamic-loading", "cuda-13000"] }
+qp-poseidon-constants = "1.1"
 
 [build-dependencies]
 qp-poseidon-constants = "1.1"

crates/engine-cuda/src/kernels/mining.cu

@@ -239,6 +239,7 @@ __device__ __forceinline__ u64 gf_sqr(u64 a) {
         : "=l"(r) : "l"(a), "r"(e));
     return r;
 }
+
 #else
 __device__ __forceinline__ u64 gf_mul(u64 a, u64 b) {
     return gf_reduce(a * b, __umul64hi(a, b));
@@ -263,7 +264,9 @@ __device__ __forceinline__ u64 gf_sqr(u64 a) {
 }
 #endif // LAIR_FUSED_MUL
 
-// x^7
+// x^7. Writing the four steps as one PTX block in 32-bit halves was tried
+// (2026-09-13) and compiles to byte-identical SASS: the 64-bit pair moves
+// between blocks are free, so there is nothing to save here.
 __device__ __forceinline__ u64 gf_sbox(u64 x) {
     u64 x2 = gf_sqr(x);
     u64 x4 = gf_sqr(x2);
@@ -310,6 +313,26 @@ __device__ __forceinline__ u64 gf_canon(u64 a) {
 #ifndef LAIR_NOINLINE_PERMUTE
 #define LAIR_NOINLINE_PERMUTE 0
 #endif
+// Threads per block; the host's MINER_CUDA_THREADS_PER_BLOCK must agree.
+#ifndef LAIR_TPB
+#define LAIR_TPB 256
+#endif
+// Decide most nonces on element 0 of the final state before computing the
+// rest of the last linear layer: at mainnet difficulty the top ~47 bits of
+// the hash must be zero, and they live in that element.
+#ifndef LAIR_EARLY_REJECT
+#define LAIR_EARLY_REJECT 1
+#endif
+// Derive each nonce's state after the first linear layer from the host's
+// precomputed base state plus one scalar times column 7 of the external
+// matrix, instead of absorbing and recomputing the layer per nonce. Falls
+// back to the full path for the nonces past a carry out of the low limb.
+#ifndef LAIR_NONCE_DIR
+#define LAIR_NONCE_DIR 1
+#endif
+#if LAIR_NONCE_DIR && !(LAIR_LAZY_ADD && LAIR_FOLD_RC && LAIR_FUSED_MUL)
+#error "LAIR_NONCE_DIR needs LAIR_LAZY_ADD, LAIR_FOLD_RC and LAIR_FUSED_MUL"
+#endif
 
 // Deferred-carry accumulator: the value is lo + hi * 2^64 with hi small. Adds
 // cost an add-with-carry instead of two compare-and-fold steps; the folds are
@@ -513,6 +536,22 @@ __device__ __forceinline__ void ext_layer_rc(u64 st[12], const u64* rc) {
 
 __device__ __forceinline__ void ext_layer(u64 st[12]) { ext_layer_rc(st, nullptr); }
 
+// Output 0 of the external linear layer alone, inputs untouched: the chunk
+// sums it needs (about a tenth of the full layer).
+__device__ __forceinline__ u64 ext_layer_out0(const u64 st[12]) {
+    Acc a[3];
+#pragma unroll
+    for (int chunk = 0; chunk < 3; chunk++) {
+        int o = chunk * 4;
+        Acc t01 = acc_add(acc_of(st[o]), st[o + 1]);
+        Acc t23 = acc_add(acc_of(st[o + 2]), st[o + 3]);
+        Acc t01123 = acc_add(acc_add2(t01, t23), st[o + 1]);
+        a[chunk] = acc_add2(t01123, t01);
+    }
+    Acc sum0 = acc_add2(acc_add2(a[0], a[1]), a[2]);
+    return acc_reduce(acc_add2(a[0], sum0));
+}
+
 // Internal linear layer: one deferred sum, folded into each diagonal multiply.
 // `rc0` is the next internal round's constant for element 0, folded into that
 // element's multiply-add; `rc_row` a full row (the first terminal round's).
@@ -566,6 +605,18 @@ __device__ __forceinline__ void ext_layer(u64 st[12]) {
     }
 }
 
+__device__ __forceinline__ u64 ext_layer_out0(const u64 st[12]) {
+    u64 a[3];
+#pragma unroll
+    for (int chunk = 0; chunk < 3; chunk++) {
+        int o = chunk * 4;
+        u64 t01 = gf_add(st[o], st[o + 1]);
+        u64 t23 = gf_add(st[o + 2], st[o + 3]);
+        a[chunk] = gf_add(gf_add(gf_add(t01, t23), st[o + 1]), t01);
+    }
+    return gf_add(a[0], gf_add(gf_add(a[0], a[1]), a[2]));
+}
+
 // Internal linear layer: diagonal matrix plus full sum.
 __device__ __forceinline__ void int_layer(u64 st[12]) {
     u64 sum = st[0];
@@ -576,16 +627,19 @@ __device__ __forceinline__ void int_layer(u64 st[12]) {
 }
 #endif
 
+// `skip_last_ext`: leave the state before the final external linear layer,
+// for the caller to finish with ext_layer_out0 / ext_layer.
 #if LAIR_NOINLINE_PERMUTE
 // One copy of the (fully unrolled) permutation instead of three inlined ones:
 // 17k instructions instead of 51k, for instruction-cache pressure.
-__device__ __noinline__ void permute(u64 st[12]) {
+__device__ __noinline__ void permute(u64 st[12], bool skip_last_ext = false, bool skip_first_ext = false) {
 #else
-__device__ __forceinline__ void permute(u64 st[12]) {
+__device__ __forceinline__ void permute(u64 st[12], bool skip_last_ext = false, bool skip_first_ext = false) {
 #endif
 #if LAIR_LAZY_ADD && LAIR_FOLD_RC
     // Each linear layer folds the constant of the round that follows it.
-    ext_layer_rc(st, RC_INITIAL[0]);
+    // `skip_first_ext`: the caller supplies the state after this layer.
+    if (!skip_first_ext) ext_layer_rc(st, RC_INITIAL[0]);
 #pragma unroll
     for (int r = 0; r < N_EXTERNAL_HALF; r++) {
 #pragma unroll
@@ -612,7 +666,7 @@ __device__ __forceinline__ void permute(u64 st[12]) {
         for (int i = 0; i < 12; i++) st[i] = gf_sbox(st[i]);
         if (r + 1 < N_EXTERNAL_HALF) {
             ext_layer_rc(st, RC_TERMINAL[r + 1]);
-        } else {
+        } else if (!skip_last_ext) {
             ext_layer(st);
         }
     }
@@ -633,7 +687,7 @@ __device__ __forceinline__ void permute(u64 st[12]) {
     for (int r = 0; r < N_EXTERNAL_HALF; r++) {
 #pragma unroll
         for (int i = 0; i < 12; i++) st[i] = gf_sbox(gf_add(st[i], RC_TERMINAL[r][i]));
-        ext_layer(st);
+        if (r + 1 < N_EXTERNAL_HALF || !skip_last_ext) ext_layer(st);
     }
 #endif
 }
@@ -659,10 +713,20 @@ struct MiningUniforms {
     u64 midstate[12];
     u32 start_nonce[16];
     u32 target[16];
+    // State after the first external layer + its round constant for the
+    // batch's first nonce (host: first_layer_after_absorb).
+    u64 layer0_base[12];
 };
-static_assert(sizeof(MiningUniforms) == 224, "MiningUniforms must match the host layout");
+static_assert(sizeof(MiningUniforms) == 320, "MiningUniforms must match the host layout");
 
-extern "C" __global__ void __launch_bounds__(256)
+#if LAIR_NONCE_DIR
+// Column 7 of the external matrix circ(2*M4, M4, M4): what the layer adds per
+// unit of element 7, which is where the low nonce limb is absorbed. Checked by
+// engine-cuda's `column_seven_of_external_matrix` test.
+__device__ __constant__ u64 LAIR_DIR7[12] = {1, 1, 3, 2, 2, 2, 6, 4, 1, 1, 3, 2};
+#endif
+
+extern "C" __global__ void __launch_bounds__(LAIR_TPB)
 mining_main(u32* __restrict__ results,
             const MiningUniforms uni,
             u32 total_threads,
@@ -699,14 +763,49 @@ mining_main(u32* __restrict__ results,
         // pad, squeeze twice (3 permutations instead of 5; the second squeeze
         // only for candidates).
         u64 st[12];
+#if LAIR_NONCE_DIR
+        if (current_nonce[1] == uni.start_nonce[1]) {
+            // Only the low limb differs from the batch's first nonce, so only
+            // element 7 of the absorbed state does, by d; the layer is linear.
+            u64 a = (u64)bswap32(current_nonce[0]);
+            u64 b = (u64)bswap32(uni.start_nonce[0]);
+            u64 d = a >= b ? a - b : a + (P64 - b);
+#pragma unroll
+            for (int i = 0; i < 12; i++) st[i] = gf_mul_add(d, LAIR_DIR7[i], uni.layer0_base[i]);
+        } else {
+            // Past a carry out of the low limb: absorb and apply the layer.
+#pragma unroll
+            for (int i = 0; i < 12; i++) st[i] = uni.midstate[i];
+#pragma unroll
+            for (int i = 0; i < 8; i++) st[i] = gf_add(st[i], (u64)bswap32(current_nonce[7 - i]));
+            ext_layer_rc(st, RC_INITIAL[0]);
+        }
+        permute(st, false, true);
+#else
 #pragma unroll
         for (int i = 0; i < 12; i++) st[i] = uni.midstate[i];
 #pragma unroll
         for (int i = 0; i < 8; i++) st[i] = gf_add(st[i], (u64)bswap32(current_nonce[7 - i]));
         permute(st);
+#endif
         st[0] = gf_add(st[0], 1ull);
         st[1] = gf_add(st[1], 1ull);
+#if LAIR_EARLY_REJECT
+        // Element 0 of the final state is the hash's top 64 bits. Compute it
+        // alone, and only nonces that do not already exceed the target there
+        // pay for the other eleven outputs of the last linear layer.
+        permute(st, true);
+        {
+            u64 c0 = gf_canon(ext_layer_out0(st));
+            u32 h0 = bswap32((u32)(c0 & EPS));
+            u32 t0 = uni.target[15];
+            if (h0 > t0) continue;
+            if (h0 == t0 && bswap32((u32)(c0 >> 32)) > uni.target[14]) continue;
+        }
+        ext_layer(st);
+#else
         permute(st);
+#endif
 
         // First squeeze: most significant 256 bits of the hash.
         u32 first[8];

crates/engine-cuda/src/lib.rs

@@ -27,7 +27,44 @@ const FATBIN: &[u8] = include_bytes!(concat!(env!("OUT_DIR"), "/mining.fatbin"))
 /// SM list the fat binary carries, for `--version` and the config metric.
 pub const CUDA_ARCHS: &str = env!("MINER_CUDA_ARCHS");
 const KERNEL_ID: &str = "cuda";
-const THREADS_PER_BLOCK: u32 = 256;
+/// Threads per block. Must match the kernel's `__launch_bounds__` (LAIR_TPB,
+/// default 256); override for experiments with MINER_CUDA_THREADS_PER_BLOCK
+/// on a kernel compiled with the same value.
+const THREADS_PER_BLOCK_DEFAULT: u32 = 256;
+
+fn threads_per_block() -> u32 {
+    std::env::var("MINER_CUDA_THREADS_PER_BLOCK")
+        .ok()
+        .and_then(|v| v.parse().ok())
+        .filter(|v: &u32| *v > 0 && v.is_multiple_of(32))
+        .unwrap_or(THREADS_PER_BLOCK_DEFAULT)
+}
+
+/// lair: size every batch to whole grids. The CLI default batch (1M nonces)
+/// fills a 5090 about eleven waves deep with a half-empty tail wave and pays
+/// the per-launch gap every 1M hashes. Rounding the batch up to a multiple of
+/// the grid (SMs x threads per SM, one nonce per thread) measured +1.6% on
+/// beast, and two grids per batch (each thread loops over two nonces) a
+/// further +0.4% there and +0.5% on the 4090, neutral on the 3060
+/// (quantus/miner#3, 2026-09-13). A batch is then ~10 ms on every card, so a
+/// job switch discards about 0.1% of a block interval. Overrides:
+/// MINER_CUDA_BATCH_ALIGN=0 for the CLI batch as given, MINER_CUDA_BATCH_WAVES
+/// for the minimum number of grids per batch.
+const BATCH_WAVES_DEFAULT: u32 = 2;
+
+fn batch_align() -> bool {
+    std::env::var("MINER_CUDA_BATCH_ALIGN")
+        .map(|v| v != "0")
+        .unwrap_or(true)
+}
+
+fn batch_waves() -> u32 {
+    std::env::var("MINER_CUDA_BATCH_WAVES")
+        .ok()
+        .and_then(|v| v.parse().ok())
+        .filter(|v: &u32| *v > 0)
+        .unwrap_or(BATCH_WAVES_DEFAULT)
+}
 /// Threads per SM the grid is sized for (override: MINER_CUDA_THREADS_PER_SM).
 /// Sets how many nonces each thread loops over for a given batch; tuned with #2.
 // Measured on the 4090 (#3): 8192 caps the grid at ~1M threads, so batches
@@ -56,6 +93,7 @@ struct Device {
     func: CudaFunction,
     name: String,
     sm_count: u32,
+    threads_per_block: u32,
     threads_per_sm: u32,
     metrics: metrics::DeviceMetrics,
 }
@@ -86,9 +124,15 @@ struct MiningUniforms {
     midstate: [u64; 12],
     start_nonce: [u32; 16],
     target: [u32; 16],
+    /// lair: the state after the first external linear layer (and its round
+    /// constant) for the batch's first nonce. That layer is linear, and only
+    /// the least significant nonce limb changes within a batch, so the kernel
+    /// derives each nonce's state from this by one scalar times a fixed
+    /// small-integer column instead of recomputing the layer (quantus/miner#3).
+    layer0_base: [u64; 12],
 }
 
-// SAFETY: `#[repr(C)]` plain-old-data with no padding (12 x u64 then 32 x u32),
+// SAFETY: `#[repr(C)]` plain-old-data with no padding (12 x u64, 32 x u32, 12 x u64),
 // matching the kernel's parameter layout; there is nothing to validate beyond
 // the layout, which is what this marker asserts.
 unsafe impl DeviceRepr for MiningUniforms {}
@@ -96,7 +140,7 @@ unsafe impl DeviceRepr for MiningUniforms {}
 // A mismatch here would silently feed the kernel the wrong midstate or target,
 // so it is a build error rather than a mining bug. The kernel carries the same
 // assertion.
-const _: () = assert!(std::mem::size_of::<MiningUniforms>() == 224);
+const _: () = assert!(std::mem::size_of::<MiningUniforms>() == 320);
 
 pub struct CudaEngine {
     engine_id: usize,
@@ -184,6 +228,7 @@ impl Device {
             name,
             sm_count,
             threads_per_sm: threads_per_sm(),
+            threads_per_block: threads_per_block(),
             metrics: metrics::DeviceMetrics::new(ordinal, KERNEL_ID),
         })
     }
@@ -340,8 +385,15 @@ impl CudaEngine {
             let headroom =
                 (U512::one() << 256) - (current_start & ((U512::one() << 256) - U512::one()));
             let cap = remaining.min(headroom);
-            let this_batch: u32 = if cap > U512::from(self.batch_size) {
+            let want: u32 = if batch_align() {
+                let grid = (dev.sm_count * dev.threads_per_sm).max(1);
+                let waves = self.batch_size.div_ceil(grid).max(batch_waves());
+                waves.saturating_mul(grid)
+            } else {
                 self.batch_size
+            };
+            let this_batch: u32 = if cap > U512::from(want) {
+                want
             } else {
                 cap.low_u32()
             };
@@ -420,8 +472,9 @@ impl CudaEngine {
         // the remainder per thread.
         let max_threads = (dev.sm_count * dev.threads_per_sm) as u64;
         let logical_threads = (batch_size as u64).min(max_threads).max(1);
-        let num_blocks = ((logical_threads as u32).div_ceil(THREADS_PER_BLOCK)).max(1);
-        let total_threads = (num_blocks * THREADS_PER_BLOCK) as u64;
+        let tpb = dev.threads_per_block;
+        let num_blocks = ((logical_threads as u32).div_ceil(tpb)).max(1);
+        let total_threads = (num_blocks * tpb) as u64;
         let nonces_per_thread = ((batch_size as u64).div_ceil(total_threads)).max(1) as u32;
         let total_threads_u32 = total_threads as u32;
 
@@ -434,10 +487,12 @@ impl CudaEngine {
         let nonce_be = batch_start.to_big_endian();
         let midstate = pow_core::mining_midstate(ctx.header, nonce_be[..32].try_into().unwrap());
 
+        let layer0_base = first_layer_after_absorb(&midstate, &start_u32s[..8]);
         let uni = MiningUniforms {
             midstate,
             start_nonce: start_u32s,
             target: res.target_u32s,
+            layer0_base,
         };
 
         let stream = res.stream.clone();
@@ -445,7 +500,7 @@ impl CudaEngine {
             stream.memset_zeros(&mut res.results)?;
             let cfg = LaunchConfig {
                 grid_dim: (num_blocks, 1, 1),
-                block_dim: (THREADS_PER_BLOCK, 1, 1),
+                block_dim: (tpb, 1, 1),
                 shared_mem_bytes: 0,
             };
             let mut launch = stream.launch_builder(&dev.func);
@@ -507,3 +562,69 @@ impl CudaEngine {
 fn bytemuck_cast(words: &[u32]) -> Vec<u8> {
     words.iter().flat_map(|w| w.to_le_bytes()).collect()
 }
+
+// lair: host-side twin of the kernel's absorb + first external linear layer,
+// in canonical Goldilocks arithmetic. Bit-exact with `ext_layer_rc` in
+// kernels/mining.cu modulo p (the kernel keeps lazy representatives).
+const GOLDILOCKS_P: u128 = 0xFFFF_FFFF_0000_0001;
+
+fn gf_add(a: u64, b: u64) -> u64 {
+    ((a as u128 + b as u128) % GOLDILOCKS_P) as u64
+}
+
+/// State after absorbing the low nonce limbs into the midstate and applying
+/// the first external layer plus the first round constant, for the batch's
+/// first nonce. `nonce_le` are the low 8 little-endian u32 limbs.
+fn first_layer_after_absorb(midstate: &[u64; 12], nonce_le: &[u32]) -> [u64; 12] {
+    let mut st = *midstate;
+    for i in 0..8 {
+        st[i] = gf_add(st[i], nonce_le[7 - i].swap_bytes() as u64);
+    }
+    let mut out = [0u64; 12];
+    for chunk in 0..3 {
+        let o = chunk * 4;
+        let (x0, x1, x2, x3) = (st[o], st[o + 1], st[o + 2], st[o + 3]);
+        let t01 = gf_add(x0, x1);
+        let t23 = gf_add(x2, x3);
+        let t0123 = gf_add(t01, t23);
+        let t01123 = gf_add(t0123, x1);
+        let t01233 = gf_add(t0123, x3);
+        out[o + 3] = gf_add(t01233, gf_add(x0, x0));
+        out[o + 1] = gf_add(t01123, gf_add(x2, x2));
+        out[o] = gf_add(t01123, t01);
+        out[o + 2] = gf_add(t01233, t23);
+    }
+    let mut sums = [0u64; 4];
+    for k in 0..4 {
+        sums[k] = gf_add(gf_add(out[k], out[k + 4]), out[k + 8]);
+    }
+    let rc = &qp_poseidon_constants::POSEIDON2_INITIAL_EXTERNAL_CONSTANTS_RAW[0];
+    for i in 0..12 {
+        out[i] = gf_add(gf_add(out[i], sums[i & 3]), rc[i]);
+    }
+    out
+}
+
+#[cfg(test)]
+mod layer0_tests {
+    use super::*;
+
+    /// The kernel's fast path assumes column 7 of the external matrix is this
+    /// vector: applying the layer to the unit vector e7 must reproduce it.
+    #[test]
+    fn column_seven_of_external_matrix() {
+        let zero = [0u64; 12];
+        // absorb puts bswap(limb 0) into element 7: choose limb 0 so that it becomes 1
+        let mut nonce = [0u32; 8];
+        nonce[0] = 1u32.swap_bytes();
+        let with = first_layer_after_absorb(&zero, &nonce);
+        let without = first_layer_after_absorb(&zero, &[0u32; 8]);
+        let rc = &qp_poseidon_constants::POSEIDON2_INITIAL_EXTERNAL_CONSTANTS_RAW[0];
+        let expect: [u64; 12] = [1, 1, 3, 2, 2, 2, 6, 4, 1, 1, 3, 2];
+        for i in 0..12 {
+            assert_eq!(without[i], rc[i]);
+            let diff = (with[i] as u128 + GOLDILOCKS_P - without[i] as u128) % GOLDILOCKS_P;
+            assert_eq!(diff as u64, expect[i], "column entry {i}");
+        }
+    }
+}