041 batches: loops run a batch of rounds at a time, each operation a C loop made at build time (no code made at run time), on Pure's forms and checks, against C natural and C given the same restructuring: roots 6.79 ms (C 6.86), knapsack fused four-byte 2.70 (C 6.47, C in the batch's blocks 1.94), matmul lanes over j 6.67 (C 7.10, C swapped and vector 4.88), join 37 / 46 with its recurrence interpreted (C 30), BFS 10.3 (C 10.0, C read ahead 7.0); unfused interpreted operations 11.2 and 17.7-22.1

4 files changed, 707 insertions, 0 deletions

diff --git a/research/experiments/041-batches/2026-10-05-batches.jsonl b/research/experiments/041-batches/2026-10-05-batches.jsonl
new file mode 100644
index 0000000..6fe3fea
--- /dev/null
+++ b/research/experiments/041-batches/2026-10-05-batches.jsonl
@@ -0,0 +1,15 @@
+{"work":"matmul","c_natural_ms":7.59,"c_ms":7.09,"c_ikj_ms":8.51,"c_ikj_u32_ms":7.41,"c_ikj_vec_ms":4.89,"c_ikj_dvec_ms":7.87,"batch_k_ms":22.08,"batch_j_ms":17.75,"axpy_ms":6.68,"same":true}
+{"work":"knapsack","c_natural_ms":6.30,"c_ms":6.47,"c_u32_ms":7.62,"c_vec_ms":9.51,"c_block_ms":1.96,"batch_ms":11.21,"u32_ms":2.71,"same":true}
+{"work":"roots","c_ms":6.86,"batch_ms":6.79,"same":true}
+{"work":"join","c_ms":30.73,"batch_ms":38.10,"interpreted_ms":51.01,"reciprocal_ms":46.73,"same":true}
+{"work":"bfs","c_ms":9.99,"c_ahead_ms":7.10,"forms_ms":10.32,"batch_ms":10.29,"reached":99998,"same":true}
+{"work":"matmul","c_natural_ms":7.62,"c_ms":7.10,"c_ikj_ms":8.56,"c_ikj_u32_ms":7.41,"c_ikj_vec_ms":4.88,"c_ikj_dvec_ms":7.94,"batch_k_ms":22.06,"batch_j_ms":17.73,"axpy_ms":6.67,"same":true}
+{"work":"knapsack","c_natural_ms":6.24,"c_ms":6.47,"c_u32_ms":7.62,"c_vec_ms":9.50,"c_block_ms":1.94,"batch_ms":11.21,"u32_ms":2.70,"same":true}
+{"work":"roots","c_ms":6.86,"batch_ms":6.79,"same":true}
+{"work":"join","c_ms":30.34,"batch_ms":37.31,"interpreted_ms":50.63,"reciprocal_ms":45.01,"same":true}
+{"work":"bfs","c_ms":10.04,"c_ahead_ms":7.00,"forms_ms":10.36,"batch_ms":10.38,"reached":99998,"same":true}
+{"work":"matmul","c_natural_ms":7.61,"c_ms":7.11,"c_ikj_ms":8.52,"c_ikj_u32_ms":7.41,"c_ikj_vec_ms":4.88,"c_ikj_dvec_ms":7.91,"batch_k_ms":22.09,"batch_j_ms":17.74,"axpy_ms":6.68,"same":true}
+{"work":"knapsack","c_natural_ms":6.25,"c_ms":6.52,"c_u32_ms":7.62,"c_vec_ms":9.50,"c_block_ms":1.94,"batch_ms":11.21,"u32_ms":2.70,"same":true}
+{"work":"roots","c_ms":6.86,"batch_ms":6.79,"same":true}
+{"work":"join","c_ms":29.37,"batch_ms":35.69,"interpreted_ms":49.15,"reciprocal_ms":44.59,"same":true}
+{"work":"bfs","c_ms":10.04,"c_ahead_ms":7.01,"forms_ms":10.39,"batch_ms":10.44,"reached":99998,"same":true}
diff --git a/research/experiments/041-batches/batch.c b/research/experiments/041-batches/batch.c
new file mode 100644
index 0000000..d4da1a7
--- /dev/null
+++ b/research/experiments/041-batches/batch.c
@@ -0,0 +1,443 @@
+/* What a batch executor could do (docs/design/024): a loop's body as a short
+ * program of operations over a batch of rounds (a value one per round, a
+ * lane), each operation a C loop the compiler makes at build time (no code
+ * made at run time), one operation interpreted per batch; Pure's checks
+ * (whole numbers below 2^53, finite numbers, cells fitting their columns)
+ * done over the batch, a flag for it (a failed batch runs again round by
+ * round in the executor). Against natural.h's C on the same inputs, called
+ * where its sizes are not constants (c_natural: it makes its cells too), and
+ * against C given what the batches are given, the cells made before (c: the
+ * natural loop; c_ikj: the loops over j and k swapped; u32: Pure's four-byte
+ * cells; c_ikj_vec, c_ikj_dvec: those, and the numbers' eight bytes, at -O3,
+ * made vector; c_vec: the natural knapsack at -O3 with ivdep, which gcc
+ * leaves scalar; c_block: knapsack written as its batch runs it).
+ *   matmul   batch_k: the inner loop's rounds a batch (strided reads of b);
+ *            batch_j: the loops over j and k swapped (each sum still in k's
+ *            order: the same numbers), a row of b a batch; axpy: that with
+ *            the multiply and the add as one operation
+ *   knapsack each item's rounds a batch, lanes of 8 bytes; u32: lanes of
+ *            the column's 4 bytes, the add, max and store one operation
+ *   roots    a batch of i, each lane's Newton rounds while it goes on
+ * Built cc -O2 -std=c17 -I../040-against-c batch.c -lm (the batch operations
+ * at -O3: their loops made vector), as join.c and bfs.c; their runs in
+ * 2026-10-05-batches.jsonl, three of each (each the median of 11 turns). */
+#define _POSIX_C_SOURCE 200809L
+#include <math.h>
+#include <stdbool.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+#include "natural.h"
+
+static double now(void) {
+    struct timespec t;
+    clock_gettime(CLOCK_MONOTONIC, &t);
+    return t.tv_sec * 1e3 + t.tv_nsec / 1e6;
+}
+
+enum { B = 256, REGS = 16 };
+/* The operations: d = a op b over n lanes; ints in i[], numbers in f[]. */
+enum {
+    O_IOTA,    /* i[d] = k + lane * k2 */
+    O_ADDK,    /* i[d] = i[a] + k */
+    O_MULK,    /* i[d] = i[a] * k (checked below 2^53) */
+    O_SUBK_R,  /* i[d] = k - i[a] */
+    O_ADD,     /* i[d] = i[a] + i[b] (checked) */
+    O_MUL,     /* i[d] = i[a] * i[b] (checked) */
+    O_MAX,     /* i[d] = max(i[a], i[b]) */
+    O_LOAD32,  /* i[d] = col[i[a]] (u32 cells) */
+    O_STORE32, /* col[i[a]] = i[b] (fits 32 bits, checked) */
+    O_SUM,     /* acc += sum of i[a] (checked) */
+    O_FROMI,   /* f[d] = i[a] */
+    O_FDIV,    /* f[d] = f[a] / f[b] (b not 0, finite) */
+    O_FADD,    /* f[d] = f[a] + f[b] (finite) */
+    O_FMULK,   /* f[d] = f[a] * k (finite) */
+};
+typedef struct { uint8_t op, d, a, b; int64_t k, k2; uint32_t *col; } Op;
+typedef struct {
+    int64_t i[REGS][B];
+    double f[REGS][B];
+    int64_t acc;
+    bool bad; /* a check failed in the batch: the executor runs it again round by round */
+} Lanes;
+
+/* Each operation over n lanes; over collects bits whose any set means a check failed (a whole number past 2^53). */
+__attribute__((optimize("O3"))) static void op_add(int64_t *restrict d, const int64_t *restrict a, const int64_t *restrict b, unsigned n, int64_t *over) {
+    int64_t o = 0;
+    for (unsigned l = 0; l < n; ++l) { int64_t x = a[l] + b[l]; d[l] = x; o |= (x >> 53) ^ (x >> 63); }
+    *over |= o;
+}
+__attribute__((optimize("O3"))) static void op_addto(int64_t *restrict d, const int64_t *restrict b, unsigned n, int64_t *over) {
+    int64_t o = 0;
+    for (unsigned l = 0; l < n; ++l) { int64_t x = d[l] + b[l]; d[l] = x; o |= (x >> 53) ^ (x >> 63); }
+    *over |= o;
+}
+__attribute__((optimize("O3"))) static void op_mul(int64_t *restrict d, const int64_t *restrict a, const int64_t *restrict b, unsigned n, int64_t *over) {
+    int64_t o = 0;
+    for (unsigned l = 0; l < n; ++l) { int64_t x = a[l] * b[l]; d[l] = x; o |= (x >> 53) ^ (x >> 63); }
+    *over |= o;
+}
+__attribute__((optimize("O3"))) static void op_mulk(int64_t *restrict d, const int64_t *restrict a, int64_t k, unsigned n, int64_t *over) {
+    int64_t o = 0;
+    for (unsigned l = 0; l < n; ++l) { int64_t x = a[l] * k; d[l] = x; o |= (x >> 53) ^ (x >> 63); }
+    *over |= o;
+}
+__attribute__((optimize("O3"))) static void op_addk(int64_t *restrict d, const int64_t *restrict a, int64_t k, unsigned n) { for (unsigned l = 0; l < n; ++l) d[l] = a[l] + k; }
+__attribute__((optimize("O3"))) static void op_max(int64_t *restrict d, const int64_t *restrict a, const int64_t *restrict b, unsigned n) {
+    for (unsigned l = 0; l < n; ++l) d[l] = a[l] > b[l] ? a[l] : b[l];
+}
+/* A column's cells at an affine position base + lane * stride (as a key affine in the loop's round is): strided reads */
+__attribute__((optimize("O3"))) static void op_loads(int64_t *restrict d, const uint32_t *restrict col, int64_t base, int64_t stride, unsigned n) {
+    if (stride == 1) for (unsigned l = 0; l < n; ++l) d[l] = col[base + l];
+    else if (stride == -1) for (unsigned l = 0; l < n; ++l) d[l] = col[base - (int64_t)l];
+    else for (unsigned l = 0; l < n; ++l) d[l] = col[base + (int64_t)l * stride];
+}
+__attribute__((optimize("O3"))) static void op_stores(uint32_t *restrict col, int64_t base, int64_t stride, const int64_t *restrict v, unsigned n, int64_t *over) {
+    int64_t o = 0;
+    for (unsigned l = 0; l < n; ++l) o |= v[l] >> 32;
+    *over |= o;
+    if (o) return;
+    if (stride == -1) for (unsigned l = 0; l < n; ++l) col[base - (int64_t)l] = (uint32_t)v[l];
+    else for (unsigned l = 0; l < n; ++l) col[base + (int64_t)l * stride] = (uint32_t)v[l];
+}
+__attribute__((optimize("O3"))) static int64_t op_sum(const int64_t *restrict a, unsigned n) { int64_t s = 0; for (unsigned l = 0; l < n; ++l) s += a[l]; return s; }
+
+enum { P_LOADS, P_ADDK, P_MUL, P_MULK, P_ADD, P_MAX, P_STORES, P_SUM };
+typedef struct { uint8_t op, d, a, b; int64_t k, base, stride; uint32_t *col; } POp;
+static void run(const POp *prog, unsigned ops, Lanes *L, unsigned n) {
+    int64_t over = 0;
+    for (unsigned p = 0; p < ops; ++p) {
+        const POp *o = &prog[p];
+        switch (o->op) {
+        case P_LOADS: op_loads(L->i[o->d], o->col, o->base, o->stride, n); break;
+        case P_ADDK: op_addk(L->i[o->d], L->i[o->a], o->k, n); break;
+        case P_MUL: op_mul(L->i[o->d], L->i[o->a], L->i[o->b], n, &over); break;
+        case P_MULK: op_mulk(L->i[o->d], L->i[o->a], o->k, n, &over); break;
+        case P_ADD: op_add(L->i[o->d], L->i[o->a], L->i[o->b], n, &over); break;
+        case P_MAX: op_max(L->i[o->d], L->i[o->a], L->i[o->b], n); break;
+        case P_STORES: op_stores(o->col, o->base, o->stride, L->i[o->a], n, &over); break;
+        case P_SUM: L->acc += op_sum(L->i[o->a], n); over |= (L->acc >> 53) ^ (L->acc >> 63); break;
+        }
+    }
+    L->bad |= over != 0;
+}
+
+static uint64_t minstd_x;
+static uint32_t minstd_next(void) { return (uint32_t)(minstd_x = minstd_x * 48271 % 2147483647); }
+
+/* matmul, its inner loop (over k) a batch: s = sum of a[i n + k] * b[k n + j] (affine reads) */
+static double matmul_k(uint32_t n, uint32_t *a, uint32_t *b, double *last) {
+    static Lanes L;
+    double total = 0;
+    POp prog[] = {
+        {P_LOADS, 2, 0, 0, 0, 0, 1, a},
+        {P_LOADS, 4, 0, 0, 0, 0, 0, b},
+        {P_MUL, 5, 2, 4, 0, 0, 0, NULL},
+        {P_SUM, 0, 5, 0, 0, 0, 0, NULL},
+    };
+    for (uint32_t i = 0; i < n; ++i)
+        for (uint32_t j = 0; j < n; ++j) {
+            L.acc = 0, L.bad = false;
+            for (uint32_t k = 0; k < n; k += B) {
+                unsigned m = n - k < B ? n - k : B;
+                prog[0].base = (int64_t)i * n + k, prog[1].base = (int64_t)k * n + j, prog[1].stride = n;
+                run(prog, 4, &L, m);
+            }
+            total += (double)L.acc, *last = (double)L.acc;
+        }
+    return total;
+}
+
+/* matmul with the loops over j and k swapped (each s[j] still summed in k's order: the same numbers): a row of b a batch */
+static double matmul_j(uint32_t n, uint32_t *a, uint32_t *b, double *last) {
+    static Lanes L;
+    static int64_t s[1024];
+    double total = 0;
+    for (uint32_t i = 0; i < n; ++i) {
+        memset(s, 0, n * sizeof(*s));
+        int64_t over = 0;
+        for (uint32_t k = 0; k < n; ++k)
+            for (uint32_t j = 0; j < n; j += B) {
+                unsigned m = n - j < B ? n - j : B;
+                op_loads(L.i[1], b, (int64_t)k * n + j, 1, m);
+                op_mulk(L.i[2], L.i[1], a[i * n + k], m, &over);
+                op_addto(s + j, L.i[2], m, &over);
+            }
+        for (uint32_t j = 0; j < n; ++j) total += (double)s[j], *last = (double)s[j];
+    }
+    return total;
+}
+
+/* knapsack: each item's rounds a batch: best[c] = max(best[c], best[c - w] + v), c = capacity - k */
+static double knapsack(uint32_t n, uint32_t capacity, const uint32_t *weight, const uint32_t *worth, uint32_t *best) {
+    static Lanes L;
+    POp prog[] = {
+        {P_LOADS, 2, 0, 0, 0, 0, -1, best},   /* best[c - w] */
+        {P_ADDK, 2, 2, 0, 0, 0, 0, NULL},     /* + worth */
+        {P_LOADS, 3, 0, 0, 0, 0, -1, best},   /* best[c] */
+        {P_MAX, 4, 2, 3, 0, 0, 0, NULL},
+        {P_STORES, 0, 4, 0, 0, 0, -1, best},
+    };
+    for (uint32_t i = 0; i < n; ++i) {
+        uint32_t rounds = capacity + 1 - weight[i];
+        for (uint32_t k = 0; k < rounds; k += B) {
+            unsigned m = rounds - k < B ? rounds - k : B;
+            int64_t c = (int64_t)capacity - k;
+            prog[0].base = c - weight[i], prog[1].k = worth[i], prog[2].base = c, prog[4].base = c;
+            run(prog, 5, &L, m);
+        }
+    }
+    return best[capacity];
+}
+
+/* roots: a batch of i, each lane's Newton rounds while it goes on (lanes out as they stop) */
+static double roots(uint32_t n, double *steps) {
+    static double r[B], x[B], next[B];
+    static uint32_t live[B];
+    double sum = 0, count = 0;
+    for (uint32_t i0 = 1; i0 <= n; i0 += B) {
+        unsigned m = n + 1 - i0 < B ? n + 1 - i0 : B, alive = m;
+        for (unsigned l = 0; l < m; ++l) r[l] = x[l] = i0 + l, live[l] = l;
+        bool bad = false;
+        for (unsigned s = 0; s < 64 && alive; ++s) {
+            /* the body as three operations over the live lanes (gathered dense), then who goes on */
+            for (unsigned q = 0; q < alive; ++q) { unsigned l = live[q]; bad |= r[l] == 0; next[q] = x[l] / r[l]; bad |= !isfinite(next[q]); }
+            for (unsigned q = 0; q < alive; ++q) { unsigned l = live[q]; next[q] = r[l] + next[q]; bad |= !isfinite(next[q]); }
+            for (unsigned q = 0; q < alive; ++q) { next[q] = next[q] * 0.5; bad |= !isfinite(next[q]); }
+            unsigned kept = 0;
+            for (unsigned q = 0; q < alive; ++q) {
+                unsigned l = live[q];
+                if (next[q] >= r[l]) continue;
+                r[l] = next[q], ++count, live[kept++] = l;
+            }
+            alive = kept;
+        }
+        (void)bad;
+        for (unsigned l = 0; l < m; ++l) sum += r[l]; /* in i's order, as run adds them */
+    }
+    *steps = count;
+    return sum;
+}
+
+
+/* Narrow lanes: cells of four bytes as they are (u32), a sum past 32 bits a failed check (the column's width) */
+__attribute__((optimize("O3"))) static void u32_addk_max_store(uint32_t *restrict col, int64_t dst, int64_t src, uint32_t k, unsigned n, uint32_t *over) {
+    /* best[c - l] = max(best[c - l], best[c - w - l] + k), lanes descending (reads below every write of the batch) */
+    uint32_t o = 0;
+    const uint32_t *from = col + src - (n - 1);
+    uint32_t *to = col + dst - (n - 1);
+    uint32_t tmp[256];
+    for (unsigned l = 0; l < n; ++l) { uint32_t x = from[l] + k; o |= x < k; uint32_t y = to[l]; tmp[l] = x > y ? x : y; }
+    *over |= o;
+    if (!o) for (unsigned l = 0; l < n; ++l) to[l] = tmp[l];
+}
+static double knapsack_u32(uint32_t n, uint32_t capacity, const uint32_t *weight, const uint32_t *worth, uint32_t *best) {
+    uint32_t over = 0;
+    for (uint32_t i = 0; i < n; ++i) {
+        uint32_t rounds = capacity + 1 - weight[i];
+        for (uint32_t k = 0; k < rounds; k += B) {
+            unsigned m = rounds - k < B ? rounds - k : B;
+            int64_t c = (int64_t)capacity - k;
+            u32_addk_max_store(best, c, c - weight[i], worth[i], m, &over);
+        }
+    }
+    return over ? -1 : best[capacity];
+}
+/* s[j..] += k * b[k n + j..]: one operation (a multiply and an add of a kept column), s as whole numbers below 2^53 */
+__attribute__((optimize("O3"))) static void axpy(int64_t *restrict s, const uint32_t *restrict b, uint32_t k, unsigned n, int64_t *over) {
+    int64_t o = 0;
+    for (unsigned l = 0; l < n; ++l) { int64_t x = s[l] + (int64_t)((uint64_t)k * b[l]); s[l] = x; o |= (x >> 53) ^ (x >> 63); }
+    *over |= o;
+}
+static double matmul_axpy(uint32_t n, uint32_t *a, uint32_t *b, double *last) {
+    static int64_t s[1024];
+    double total = 0;
+    int64_t over = 0;
+    for (uint32_t i = 0; i < n; ++i) {
+        memset(s, 0, n * sizeof(*s));
+        for (uint32_t k = 0; k < n; ++k)
+            for (uint32_t j = 0; j < n; j += B) {
+                unsigned m = n - j < B ? n - j : B;
+                axpy(s + j, b + (size_t)k * n + j, a[i * n + k], m, &over);
+            }
+        for (uint32_t j = 0; j < n; ++j) total += (double)s[j], *last = (double)s[j];
+    }
+    return over ? -1 : total;
+}
+
+/* The C given what the batches were given, the cells made before as the batches' are:
+ * the loops over j and k swapped, Pure's four-byte cells, -O3 with ivdep (vector). */
+__attribute__((noinline)) static double c_ijk(uint32_t n, const double *a, const double *b, double *last) {
+    double total = 0;
+    for (uint32_t i = 0; i < n; ++i)
+        for (uint32_t j = 0; j < n; ++j) {
+            double s = 0;
+            for (uint32_t k = 0; k < n; ++k) s += a[i * n + k] * b[k * n + j];
+            total += s, *last = s;
+        }
+    return total;
+}
+__attribute__((noinline)) static double c_ikj(uint32_t n, const double *a, const double *b, double *last) {
+    static double s[1024];
+    double total = 0;
+    for (uint32_t i = 0; i < n; ++i) {
+        memset(s, 0, n * sizeof(*s));
+        for (uint32_t k = 0; k < n; ++k)
+            for (uint32_t j = 0; j < n; ++j) s[j] += a[i * n + k] * b[k * n + j];
+        for (uint32_t j = 0; j < n; ++j) total += s[j], *last = s[j];
+    }
+    return total;
+}
+__attribute__((noinline)) static double c_ikj_u32(uint32_t n, const uint32_t *a, const uint32_t *b, double *last) {
+    static int64_t s[1024];
+    double total = 0;
+    for (uint32_t i = 0; i < n; ++i) {
+        memset(s, 0, n * sizeof(*s));
+        for (uint32_t k = 0; k < n; ++k)
+            for (uint32_t j = 0; j < n; ++j) s[j] += (int64_t)a[i * n + k] * b[k * n + j];
+        for (uint32_t j = 0; j < n; ++j) total += (double)s[j], *last = (double)s[j];
+    }
+    return total;
+}
+__attribute__((noinline, optimize("O3"))) static double c_ikj_vec(uint32_t n, const uint32_t *restrict a, const uint32_t *restrict b, double *last) {
+    static int64_t s[1024];
+    double total = 0;
+    for (uint32_t i = 0; i < n; ++i) {
+        memset(s, 0, n * sizeof(*s));
+        for (uint32_t k = 0; k < n; ++k)
+            for (uint32_t j = 0; j < n; ++j) s[j] += (int64_t)a[i * n + k] * b[k * n + j];
+        for (uint32_t j = 0; j < n; ++j) total += (double)s[j], *last = (double)s[j];
+    }
+    return total;
+}
+__attribute__((noinline, optimize("O3"))) static double c_ikj_dvec(uint32_t n, const double *restrict a, const double *restrict b, double *last) {
+    static double s[1024];
+    double total = 0;
+    for (uint32_t i = 0; i < n; ++i) {
+        memset(s, 0, n * sizeof(*s));
+        for (uint32_t k = 0; k < n; ++k)
+            for (uint32_t j = 0; j < n; ++j) s[j] += a[i * n + k] * b[k * n + j];
+        for (uint32_t j = 0; j < n; ++j) total += s[j], *last = s[j];
+    }
+    return total;
+}
+__attribute__((noinline)) static double c_knap_i64(uint32_t n, uint32_t capacity, const uint32_t *weight, const uint32_t *worth, int64_t *best) {
+    for (uint32_t i = 0; i < n; ++i)
+        for (uint32_t c = capacity; c >= weight[i] && c <= capacity; --c) {
+            int64_t with = best[c - weight[i]] + worth[i];
+            if (with > best[c]) best[c] = with;
+        }
+    return (double)best[capacity];
+}
+__attribute__((noinline)) static double c_knap_u32(uint32_t n, uint32_t capacity, const uint32_t *weight, const uint32_t *worth, uint32_t *best) {
+    for (uint32_t i = 0; i < n; ++i)
+        for (uint32_t c = capacity; c >= weight[i] && c <= capacity; --c) {
+            uint32_t with = best[c - weight[i]] + worth[i];
+            if (with > best[c]) best[c] = with;
+        }
+    return best[capacity];
+}
+__attribute__((noinline, optimize("O3", "vect-cost-model=unlimited"))) static double c_knap_vec(uint32_t n, uint32_t capacity, const uint32_t *weight, const uint32_t *worth, uint32_t *best) {
+    for (uint32_t i = 0; i < n; ++i) {
+        uint32_t w = weight[i], v = worth[i], rounds = capacity + 1 - w;
+        /* best[c] for c descending reads best[c - w], written later in that order: no round reads another's write */
+#pragma GCC ivdep
+        for (uint32_t k = 0; k < rounds; ++k) {
+            uint32_t c = capacity - k, with = best[c - w] + v;
+            best[c] = with > best[c] ? with : best[c];
+        }
+    }
+    return best[capacity];
+}
+
+/* knapsack in C as the batch runs it: 256 rounds at a time, read, then written (ascending within the block) */
+__attribute__((noinline, optimize("O3"))) static double c_knap_block(uint32_t n, uint32_t capacity, const uint32_t *weight, const uint32_t *worth, uint32_t *restrict best) {
+    uint32_t tmp[256];
+    for (uint32_t i = 0; i < n; ++i) {
+        uint32_t w = weight[i], v = worth[i], rounds = capacity + 1 - w;
+        for (uint32_t k = 0; k < rounds; k += 256) {
+            uint32_t m = rounds - k < 256 ? rounds - k : 256, low = capacity - k - (m - 1);
+            const uint32_t *from = best + low - w;
+            uint32_t *to = best + low;
+            for (uint32_t l = 0; l < m; ++l) { uint32_t with = from[l] + v, was = to[l]; tmp[l] = with > was ? with : was; }
+            for (uint32_t l = 0; l < m; ++l) to[l] = tmp[l];
+        }
+    }
+    return best[capacity];
+}
+
+static int cmpd(const void *a, const void *b) { double x = *(const double *)a, y = *(const double *)b; return (x > y) - (x < y); }
+#define MEDIAN(arr) (qsort(arr, 11, sizeof(double), cmpd), arr[5])
+#define TIME(slot, call) do { double t_ = now(); call; (slot) = now() - t_; } while (0)
+static volatile uint32_t sizes[5] = {200, 100, 50000, 100000, 0}; /* read at run time: no C sees them as constants */
+static void (*volatile natural[3])(const uint32_t *, double *) = {c_matmul, c_knapsack, c_roots};
+
+int main(void) {
+    /* matmul, 200 x 200, cells as natural.h makes them */
+    uint32_t n = sizes[0];
+    uint32_t *a = malloc(4 * n * n), *b = malloc(4 * n * n);
+    double *da = malloc(8 * n * n), *db = malloc(8 * n * n);
+    minstd_x = 1;
+    for (uint32_t i = 0; i < n * n; ++i) a[i] = minstd_next() % 100, b[i] = minstd_next() % 100, da[i] = a[i], db[i] = b[i];
+    double tc[11], tn[11], ti[11], tu[11], tv[11], td[11], tk[11], tj[11], tx[11], out[4], last;
+    uint32_t in[2] = {n, 0};
+    double want = 0, got[8];
+    for (int r = 0; r < 11; ++r) {
+        TIME(tc[r], natural[0](in, out));
+        want = out[0];
+        TIME(tn[r], got[0] = c_ijk(n, da, db, &last));
+        TIME(ti[r], got[1] = c_ikj(n, da, db, &last));
+        TIME(tu[r], got[2] = c_ikj_u32(n, a, b, &last));
+        TIME(tv[r], got[6] = c_ikj_vec(n, a, b, &last));
+        TIME(td[r], got[7] = c_ikj_dvec(n, da, db, &last));
+        TIME(tk[r], got[3] = matmul_k(n, a, b, &last));
+        TIME(tj[r], got[4] = matmul_j(n, a, b, &last));
+        TIME(tx[r], got[5] = matmul_axpy(n, a, b, &last));
+    }
+    bool same = true;
+    for (int g = 0; g < 8; ++g) same &= got[g] == want;
+    printf("{\"work\":\"matmul\",\"c_natural_ms\":%.2f,\"c_ms\":%.2f,\"c_ikj_ms\":%.2f,\"c_ikj_u32_ms\":%.2f,\"c_ikj_vec_ms\":%.2f,\"c_ikj_dvec_ms\":%.2f,\"batch_k_ms\":%.2f,\"batch_j_ms\":%.2f,\"axpy_ms\":%.2f,\"same\":%s}\n",
+           MEDIAN(tc), MEDIAN(tn), MEDIAN(ti), MEDIAN(tu), MEDIAN(tv), MEDIAN(td), MEDIAN(tk), MEDIAN(tj), MEDIAN(tx), same ? "true" : "false");
+    /* knapsack, 100 items, capacity 50,000 */
+    uint32_t items = sizes[1], capacity = sizes[2], *weight = malloc(4 * items), *worth = malloc(4 * items), *best = malloc(4 * (capacity + 1));
+    int64_t *best64 = malloc(8 * (capacity + 1));
+    minstd_x = 1;
+    for (uint32_t i = 0; i < items; ++i) {
+        weight[i] = minstd_next() % 1000 + 1;
+        worth[i] = minstd_next() % 1000;
+    }
+    double kc[11], k6[11], k3[11], kv[11], kk[11], kb[11], ku[11], kgot[6];
+    uint32_t kin[2] = {items, capacity};
+    for (int r = 0; r < 11; ++r) {
+        TIME(kc[r], natural[1](kin, out));
+        want = out[0];
+        memset(best64, 0, 8 * (capacity + 1));
+        TIME(k6[r], kgot[0] = c_knap_i64(items, capacity, weight, worth, best64));
+        memset(best, 0, 4 * (capacity + 1));
+        TIME(k3[r], kgot[1] = c_knap_u32(items, capacity, weight, worth, best));
+        memset(best, 0, 4 * (capacity + 1));
+        TIME(kv[r], kgot[2] = c_knap_vec(items, capacity, weight, worth, best));
+        memset(best, 0, 4 * (capacity + 1));
+        memset(best, 0, 4 * (capacity + 1));
+        TIME(kk[r], kgot[5] = c_knap_block(items, capacity, weight, worth, best));
+        memset(best, 0, 4 * (capacity + 1));
+        TIME(kb[r], kgot[3] = knapsack(items, capacity, weight, worth, best));
+        memset(best, 0, 4 * (capacity + 1));
+        TIME(ku[r], kgot[4] = knapsack_u32(items, capacity, weight, worth, best));
+    }
+    same = true;
+    for (int g = 0; g < 6; ++g) same &= kgot[g] == want;
+    printf("{\"work\":\"knapsack\",\"c_natural_ms\":%.2f,\"c_ms\":%.2f,\"c_u32_ms\":%.2f,\"c_vec_ms\":%.2f,\"c_block_ms\":%.2f,\"batch_ms\":%.2f,\"u32_ms\":%.2f,\"same\":%s}\n",
+           MEDIAN(kc), MEDIAN(k6), MEDIAN(k3), MEDIAN(kv), MEDIAN(kk), MEDIAN(kb), MEDIAN(ku), same ? "true" : "false");
+    /* roots, 100,000 */
+    double rc[11], rb[11], rgot = 0, steps;
+    uint32_t rin[1] = {sizes[3]};
+    for (int r = 0; r < 11; ++r) {
+        TIME(rc[r], natural[2](rin, out));
+        want = out[0];
+        TIME(rb[r], rgot = roots(rin[0], &steps));
+    }
+    printf("{\"work\":\"roots\",\"c_ms\":%.2f,\"batch_ms\":%.2f,\"same\":%s}\n", MEDIAN(rc), MEDIAN(rb), rgot == want && steps == out[1] ? "true" : "false");
+}
diff --git a/research/experiments/041-batches/bfs.c b/research/experiments/041-batches/bfs.c
new file mode 100644
index 0000000..3c25a98
--- /dev/null
+++ b/research/experiments/041-batches/bfs.c
@@ -0,0 +1,114 @@
+/* 033's breadth first search (graph.aup's route) as a batch executor could
+ * run it, on the forms Pure gives its lists here (docs/design/021, 022): edges
+ * packed (a column per field in key order, from * 2^20 + to; an index of
+ * buckets of 2^20 keys: a node's edges from index[node] to index[node + 1]),
+ * seen DENSE (keyed by node: a bit per key, the depth at the key, four
+ * bytes), queue FLAT (appended to).
+ *   forms  round by round, as the executor runs the loop
+ *   batch  the rounds whose test (i below queue.count) holds when the batch
+ *          opens (count only grows): each step for every lane (queue read,
+ *          the depth read, the edges walked into lanes of their own, seen's
+ *          bit read for each), then the adds and appends in order, each
+ *          lane's bit read again first (an earlier lane may have added it)
+ * Against natural.h's c_bfs (offsets per node: CSR) and c_bfs_ahead (that
+ * read ahead), on the suite's graph and routes (100,000 nodes, 1,000,000
+ * edges drawn, seed 121). Built cc -O2. */
+#define _POSIX_C_SOURCE 200809L
+#include <math.h>
+#include <stdbool.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+#include "natural.h"
+
+static double now(void) {
+    struct timespec t;
+    clock_gettime(CLOCK_MONOTONIC, &t);
+    return t.tv_sec * 1e3 + t.tv_nsec / 1e6;
+}
+
+enum { B = 256, LANES = 8192 };
+
+/* seen: list of Seen by node, DENSE */
+typedef struct { uint64_t *bits; uint32_t *depth; size_t count; } Seen;
+static void seen_start(Seen *s, uint32_t nodes) {
+    s->bits = calloc((nodes + 63) / 64, 8), s->depth = malloc(4ull * nodes), s->count = 0;
+}
+static bool seen_has(const Seen *s, uint32_t node) { return s->bits[node / 64] >> (node % 64) & 1; }
+static void seen_add(Seen *s, uint32_t node, uint32_t depth) {
+    s->bits[node / 64] |= UINT64_C(1) << (node % 64), s->depth[node] = depth, ++s->count;
+}
+
+/* edges: list of Edge by link, packed */
+static uint32_t *to_col, *index_col;
+
+/* round by round, as the executor runs the loop */
+static void forms_bfs(uint32_t from, uint32_t target, uint32_t nodes, double *out) {
+    Seen seen;
+    seen_start(&seen, nodes);
+    uint32_t *queue = malloc(4ull * nodes), count = 0;
+    queue[count++] = from;
+    seen_add(&seen, from, 0);
+    for (uint32_t i = 0; i < nodes; ++i) {
+        if (i >= count) break;
+        uint32_t node = queue[i];
+        for (uint32_t k = index_col[node]; k < index_col[node + 1]; ++k) {
+            uint32_t to = to_col[k];
+            if (!seen_has(&seen, to)) seen_add(&seen, to, seen.depth[node] + 1), queue[count++] = to;
+        }
+    }
+    out[0] = seen_has(&seen, target) ? seen.depth[target] : -1, out[1] = (double)seen.count;
+    free(seen.bits), free(seen.depth), free(queue);
+}
+
+/* a batch of rounds at a time */
+static void batch_bfs(uint32_t from, uint32_t target, uint32_t nodes, double *out) {
+    static uint32_t node[B], depth[B], lo[B], hi[B], edge_to[LANES], edge_round[LANES];
+    static bool has[LANES];
+    Seen seen;
+    seen_start(&seen, nodes);
+    uint32_t *queue = malloc(4ull * nodes), count = 0;
+    queue[count++] = from;
+    seen_add(&seen, from, 0);
+    for (uint32_t i = 0; i < nodes && i < count;) {
+        unsigned m = count - i < B ? count - i : B; /* the rounds whose test holds now (count only grows) */
+        if (m > nodes - i) m = nodes - i;
+        for (unsigned l = 0; l < m; ++l) node[l] = queue[i + l];
+        for (unsigned l = 0; l < m; ++l) lo[l] = index_col[node[l]], hi[l] = index_col[node[l] + 1];
+        for (unsigned l = 0; l < m; ++l) depth[l] = seen.depth[node[l]] + 1; /* seen[q.node].depth + 1 */
+        unsigned e = 0, rounds = 0; /* the edges into lanes, as many rounds as fit */
+        for (; rounds < m && e + (hi[rounds] - lo[rounds]) <= LANES; ++rounds)
+            for (uint32_t k = lo[rounds]; k < hi[rounds]; ++k) edge_to[e] = to_col[k], edge_round[e] = rounds, ++e;
+        for (unsigned q = 0; q < e; ++q) has[q] = seen_has(&seen, edge_to[q]); /* seen.has(e.to), each lane */
+        for (unsigned q = 0; q < e; ++q) /* the adds and appends, in order */
+            if (!has[q] && !seen_has(&seen, edge_to[q])) seen_add(&seen, edge_to[q], depth[edge_round[q]]), queue[count++] = edge_to[q];
+        i += rounds;
+    }
+    out[0] = seen_has(&seen, target) ? seen.depth[target] : -1, out[1] = (double)seen.count;
+    free(seen.bits), free(seen.depth), free(queue);
+}
+
+static int cmpd(const void *a, const void *b) { double x = *(const double *)a, y = *(const double *)b; return (x > y) - (x < y); }
+static void (*volatile natural[2])(const uint32_t *, double *) = {c_bfs, c_bfs_ahead};
+int main(void) {
+    uint64_t *keys;
+    size_t edge_count = graph_build(100000, 1000000, &keys);
+    to_col = malloc(4 * edge_count), index_col = calloc(100001, 4);
+    for (size_t k = 0; k < edge_count; ++k) to_col[k] = (uint32_t)(keys[k] & 0xFFFFF), ++index_col[(keys[k] >> 20) + 1];
+    for (uint32_t v = 0; v < 100000; ++v) index_col[v + 1] += index_col[v];
+    double tc[11], ta[11], tf[11], tb[11], oc[2], oa[2], of[2], ob[2];
+    bool same = true;
+    for (int r = 0; r < 11; ++r) {
+        uint32_t in[2] = {(uint32_t)(graph_next() % 100000), (uint32_t)(graph_next() % 100000)};
+        double t = now(); natural[0](in, oc); tc[r] = now() - t;
+        t = now(); natural[1](in, oa); ta[r] = now() - t;
+        t = now(); forms_bfs(in[0], in[1], 100000, of); tf[r] = now() - t;
+        t = now(); batch_bfs(in[0], in[1], 100000, ob); tb[r] = now() - t;
+        for (int k = 0; k < 2; ++k) same &= oc[k] == oa[k] && oc[k] == of[k] && oc[k] == ob[k];
+    }
+    qsort(tc, 11, 8, cmpd), qsort(ta, 11, 8, cmpd), qsort(tf, 11, 8, cmpd), qsort(tb, 11, 8, cmpd);
+    printf("{\"work\":\"bfs\",\"c_ms\":%.2f,\"c_ahead_ms\":%.2f,\"forms_ms\":%.2f,\"batch_ms\":%.2f,\"reached\":%.0f,\"same\":%s}\n", tc[5], ta[5], tf[5], tb[5], ob[1],
+           same ? "true" : "false");
+}
diff --git a/research/experiments/041-batches/join.c b/research/experiments/041-batches/join.c
new file mode 100644
index 0000000..35b8462
--- /dev/null
+++ b/research/experiments/041-batches/join.c
@@ -0,0 +1,135 @@
+/* join as a batch executor could run it, on Pure's own forms: all (FLAT, ids as
+ * u32 cells), byid (HASHED: key bits in a table of slots 4/5 full at most, the
+ * region at its key's slot), totals (DENSE, 16 sums). Both loops a batch of
+ * rounds at a time: the MINSTD recurrence in order (in C, or interpreted:
+ * the executor's part), then each step of the round for every lane (the
+ * setup: the appends, the adds in order; the orders: ids gathered, slots
+ * probed, regions gathered, the sums added in order). Against natural.h's
+ * c_join (its own hash table), both timed from nothing to the answer, the
+ * sizes read at run time. */
+#define _POSIX_C_SOURCE 200809L
+#include <math.h>
+#include <stdbool.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+#include "natural.h"
+static double now(void) { struct timespec t; clock_gettime(CLOCK_MONOTONIC, &t); return t.tv_sec * 1e3 + t.tv_nsec / 1e6; }
+static const uint64_t VACANT = 0x7FF8000000000001ull;
+static uint64_t key_bits(double k) { if (k == 0) k = 0; uint64_t b; memcpy(&b, &k, 8); return b; }
+static size_t slot_of(uint64_t bits, size_t mask) { bits ^= bits >> 33; bits *= 0xff51afd7ed558ccdull; bits ^= bits >> 33; return bits & mask; }
+enum { B = 256 };
+/* A section that cannot run as a batch (each round needs the last's x) run
+ * round by round by an interpreter of the same operations (computed goto,
+ * Pure's checks): x = x * 48271 - floor(x * 48271 / 2147483647) * 2147483647,
+ * each round's x into its lane. */
+enum { S_MULK, S_FLOORDIVK, S_SUB, S_PUT, S_END };
+typedef struct { uint8_t op, d, a, b; int64_t k; uint64_t magic; unsigned shift; } SOp; /* magic, shift: k's reciprocal, made with the program */
+static SOp minstd_ops[] = {
+    {S_MULK, 1, 0, 0, 48271},          /* t = x * 48271 */
+    {S_FLOORDIVK, 2, 1, 0, 2147483647}, /* q = floor(t / 2147483647) */
+    {S_MULK, 2, 2, 0, 2147483647},      /* q * 2147483647 */
+    {S_SUB, 0, 1, 2, 0},                /* x = t - that */
+    {S_PUT, 0, 0, 0, 0},                /* the lane's x */
+    {S_END, 0, 0, 0, 0},
+};
+static void rounds_prepare(SOp *prog, bool magic) { /* a divisor's reciprocal, k not a power of two: q = (x m) >> (64 + s), 2^s < k, exact for x below 2^53 */
+    for (; prog->op != S_END; ++prog)
+        if (prog->op == S_FLOORDIVK && prog->k > 2 && (prog->k & (prog->k - 1)) && magic) {
+            unsigned s = 0;
+            while ((2ull << s) < (uint64_t)prog->k) ++s;
+            prog->magic = (uint64_t)((((unsigned __int128)1 << (64 + s)) / (uint64_t)prog->k) + 1), prog->shift = s;
+        } else prog->magic = 0;
+}
+__attribute__((noinline)) static bool rounds_run(const SOp *prog, int64_t *r, uint64_t *lanes, unsigned m) {
+    static void *const at[] = {&&mulk, &&floordivk, &&sub, &&put, &&end};
+    int64_t over = 0;
+    for (unsigned l = 0; l < m; ++l) {
+        const SOp *o = prog;
+        goto *at[o->op];
+    mulk: { int64_t x = r[o->a] * o->k; over |= (x >> 53) ^ (x >> 63); r[o->d] = x; ++o; goto *at[o->op]; }
+    floordivk: {
+        int64_t x = r[o->a], q;
+        if (x >= 0 && o->magic) q = (int64_t)(((unsigned __int128)(uint64_t)x * o->magic) >> 64 >> o->shift); /* below 2^53 */
+        else q = x / o->k, q -= (x % o->k != 0) & ((x < 0) != (o->k < 0));
+        r[o->d] = q; ++o; goto *at[o->op];
+    }
+    sub: { int64_t x = r[o->a] - r[o->b]; over |= (x >> 53) ^ (x >> 63); r[o->d] = x; ++o; goto *at[o->op]; }
+    put: lanes[l] = (uint64_t)r[o->a]; ++o; goto *at[o->op];
+    end: ;
+    }
+    return over != 0;
+}
+static void pure_join(uint32_t customers, uint32_t orders, int interpreted, double *out) {
+    rounds_prepare(minstd_ops, interpreted == 2);
+    uint32_t *ids = malloc(4 * (size_t)customers);
+    size_t cap = 16;
+    while (5 * (size_t)customers > 4 * cap) cap *= 2;
+    uint64_t *members = malloc(8 * cap);
+    uint32_t *region = malloc(4 * cap);
+    for (size_t i = 0; i < cap; ++i) members[i] = VACANT;
+    uint64_t x = 1;
+    int64_t totals[16] = {0};
+    static uint64_t X[B], C[B], BITS[B], SLOT[B];
+    static uint32_t R[B];
+    bool bad = false;
+    for (uint32_t i0 = 0; i0 < customers; i0 += B) { /* the setup: the recurrence, then the batch's appends, then its adds in order */
+        unsigned m = customers - i0 < B ? customers - i0 : B;
+        if (interpreted) {
+            int64_t r[3] = {(int64_t)x};
+            bad |= rounds_run(minstd_ops, r, X, m);
+            x = (uint64_t)r[0];
+        } else
+            for (unsigned l = 0; l < m; ++l) x = x * 48271 % 2147483647, X[l] = x;
+        for (unsigned l = 0; l < m; ++l) ids[i0 + l] = (uint32_t)X[l];
+        for (unsigned l = 0; l < m; ++l) {
+            uint64_t bits = key_bits((double)X[l]);
+            size_t at = slot_of(bits, cap - 1);
+            while (members[at] != VACANT) at = (at + 1) & (cap - 1);
+            members[at] = bits, region[at] = (i0 + l) % 16;
+        }
+    }
+    for (uint32_t o0 = 0; o0 < orders; o0 += B) {
+        unsigned m = orders - o0 < B ? orders - o0 : B;
+        if (interpreted) {                                                                 /* the recurrence, in order */
+            int64_t r[3] = {(int64_t)x};
+            bad |= rounds_run(minstd_ops, r, X, m);
+            x = (uint64_t)r[0];
+        } else
+            for (unsigned l = 0; l < m; ++l) x = x * 48271 % 2147483647, X[l] = x;
+        for (unsigned l = 0; l < m; ++l) C[l] = X[l] % customers;                          /* first + (x - floor(x / n) n) */
+        for (unsigned l = 0; l < m; ++l) BITS[l] = key_bits((double)ids[C[l]]);           /* c.id, gathered */
+        for (unsigned l = 0; l < m; ++l) {                                                 /* byid[c.id]: each lane's slot */
+            size_t at = slot_of(BITS[l], cap - 1);
+            while (members[at] != BITS[l]) {
+                if (members[at] == VACANT) { bad = true; break; }
+                at = (at + 1) & (cap - 1);
+            }
+            SLOT[l] = at;
+        }
+        for (unsigned l = 0; l < m; ++l) R[l] = region[SLOT[l]];                           /* .region, gathered */
+        for (unsigned l = 0; l < m; ++l) totals[R[l]] += (int64_t)(X[l] % 1000);           /* the sums, in order */
+    }
+    int64_t most = 0, sum = 0;
+    for (int r = 0; r < 16; ++r) most = totals[r] > most ? totals[r] : most, sum += totals[r];
+    out[0] = bad ? -1 : (double)most, out[1] = (double)sum;
+    free(ids), free(members), free(region);
+}
+static int cmpd(const void *a, const void *b) { double x = *(const double *)a, y = *(const double *)b; return (x > y) - (x < y); }
+static volatile uint32_t sizes[2] = {100000, 1000000}; /* read at run time */
+static void (*volatile natural)(const uint32_t *, double *) = c_join;
+int main(void) {
+    double tc[11], tp[11], ti[11], tm[11], oc[2], op[2], oi[2], om[2];
+    uint32_t in[2] = {sizes[0], sizes[1]};
+    for (int r = 0; r < 11; ++r) {
+        double t = now(); natural(in, oc); tc[r] = now() - t;
+        t = now(); pure_join(in[0], in[1], 0, op); tp[r] = now() - t;
+        t = now(); pure_join(in[0], in[1], 1, oi); ti[r] = now() - t;
+        t = now(); pure_join(in[0], in[1], 2, om); tm[r] = now() - t;
+    }
+    qsort(tc, 11, 8, cmpd), qsort(tp, 11, 8, cmpd), qsort(ti, 11, 8, cmpd), qsort(tm, 11, 8, cmpd);
+    printf("{\"work\":\"join\",\"c_ms\":%.2f,\"batch_ms\":%.2f,\"interpreted_ms\":%.2f,\"reciprocal_ms\":%.2f,\"same\":%s}\n", tc[5], tp[5], ti[5], tm[5],
+           oc[0] == op[0] && oc[1] == op[1] && oc[0] == oi[0] && oc[1] == oi[1] && oc[0] == om[0] && oc[1] == om[1] ? "true" : "false");
+}