040 forms.c: matmul's inner loop in the forms the tier may make it (index 20.7 ms, pointer with multiply-add 13.1, apart 10.0, with checks 10.2, a bit tested before each load 13.4): the position's multiply and the pre-load branch cost, the checks hardly; 023 section 5 says so and what comes first (pointer induction, then a full list's copy)

2 files changed, 93 insertions, 0 deletions

diff --git a/research/experiments/040-against-c/2026-10-05-forms.jsonl b/research/experiments/040-against-c/2026-10-05-forms.jsonl
new file mode 100644
index 0000000..8d7ad78
--- /dev/null
+++ b/research/experiments/040-against-c/2026-10-05-forms.jsonl
@@ -0,0 +1 @@
+{"index":20.69,"plain":20.11,"pointer":13.12,"apart":10.02,"checks":10.24,"bits":13.37,"c":3.76,"same":true}
diff --git a/research/experiments/040-against-c/forms.c b/research/experiments/040-against-c/forms.c
new file mode 100644
index 0000000..93c3e9b
--- /dev/null
+++ b/research/experiments/040-against-c/forms.c
@@ -0,0 +1,92 @@
+/* matmul's inner loop (200 x 200 cells of u32) in the forms the tier may
+ * make it, as AArch64 assembly, against natural.h's C on the same Pi
+ * (research 124, docs/design/023 section 5):
+ *   index    the cell at column + (i * n + k) and column + (k * n + j), each
+ *            found by a multiply and an add; the product and the sum each
+ *            checked below 2^53; a full list's bit tested before each load
+ *            (the tier's code at e4b4694)
+ *   plain    index, no checks, no bit tested, multiply-add into the sum
+ *   pointer  a pointer to each column stepped by its stride, multiply-add
+ *   apart    pointer, the multiply and the add apart (the sum's chain an add)
+ *   checks   apart, the product and the sum checked below 2^53
+ *   bits     apart, a full list's bit tested before each load
+ * Each sums every product (the answers must be equal); the median of 11. */
+#define _POSIX_C_SOURCE 200809L
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+enum { N = 200 };
+static uint32_t A[N * N], B[N * N];
+static double DA[N * N], DB[N * N];
+static double now(void) {
+    struct timespec t;
+    clock_gettime(CLOCK_MONOTONIC, &t);
+    return t.tv_sec * 1e3 + t.tv_nsec / 1e6;
+}
+#define LOOP(name, body, ...)                                                                                          \
+    static uint64_t name(void) {                                                                                       \
+        uint64_t total = 0, flag = 1;                                                                                  \
+        for (uint64_t i = 0; i < N; ++i)                                                                               \
+            for (uint64_t j = 0; j < N; ++j) {                                                                         \
+                uint64_t s = 0, k = 0, base = i * N;                                                                   \
+                const uint32_t *pa = A + i * N, *pb = B + j;                                                           \
+                __asm__ volatile(body : [s] "+r"(s), [k] "+r"(k), [pa] "+r"(pa), [pb] "+r"(pb)                         \
+                                 : [base] "r"(base), [a] "r"(A), [b] "r"(B), [n] "r"((uint64_t)N), [j] "r"(j),           \
+                                   [flag] "r"(flag)                                                                    \
+                                 : "x10", "x20", "x24", "x27", "cc", "memory");                                        \
+                total += s;                                                                                            \
+            }                                                                                                          \
+        return total;                                                                                                  \
+    }
+LOOP(index, "1: add x20, %[base], %[k]\n tbz %w[flag], #0, 9f\n ldr w24, [%[a], x20, lsl #2]\n mul x20, %[k], %[n]\n"
+            " add x20, x20, %[j]\n tbz %w[flag], #0, 9f\n ldr w27, [%[b], x20, lsl #2]\n mul x20, x24, x27\n"
+            " lsr x10, x20, #53\n cbnz x10, 9f\n add %[s], %[s], x20\n lsr x10, %[s], #53\n cbnz x10, 9f\n"
+            " add %[k], %[k], #1\n cmp %[k], %[n]\n b.lo 1b\n 9:\n")
+LOOP(plain, "1: add x20, %[base], %[k]\n ldr w24, [%[a], x20, lsl #2]\n mul x20, %[k], %[n]\n add x20, x20, %[j]\n"
+            " ldr w27, [%[b], x20, lsl #2]\n madd %[s], x24, x27, %[s]\n add %[k], %[k], #1\n cmp %[k], %[n]\n b.lo 1b\n")
+LOOP(pointer, "1: ldr w24, [%[pa]], #4\n ldr w27, [%[pb]]\n add %[pb], %[pb], #800\n madd %[s], x24, x27, %[s]\n"
+              " add %[k], %[k], #1\n cmp %[k], %[n]\n b.lo 1b\n")
+LOOP(apart, "1: ldr w24, [%[pa]], #4\n ldr w27, [%[pb]]\n add %[pb], %[pb], #800\n mul x20, x24, x27\n add %[s], %[s], x20\n"
+            " add %[k], %[k], #1\n cmp %[k], %[n]\n b.lo 1b\n")
+LOOP(checks, "1: ldr w24, [%[pa]], #4\n ldr w27, [%[pb]]\n add %[pb], %[pb], #800\n mul x20, x24, x27\n lsr x10, x20, #53\n"
+             " cbnz x10, 9f\n add %[s], %[s], x20\n lsr x10, %[s], #53\n cbnz x10, 9f\n add %[k], %[k], #1\n"
+             " cmp %[k], %[n]\n b.lo 1b\n 9:\n")
+LOOP(bits, "1: tbz %w[flag], #0, 9f\n ldr w24, [%[pa]], #4\n tbz %w[flag], #0, 9f\n ldr w27, [%[pb]]\n add %[pb], %[pb], #800\n"
+           " mul x20, x24, x27\n add %[s], %[s], x20\n add %[k], %[k], #1\n cmp %[k], %[n]\n b.lo 1b\n 9:\n")
+static uint64_t c(void) { /* natural.h's c_matmul, its sum */
+    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 += DA[i * N + k] * DB[k * N + j];
+            total += s;
+        }
+    return (uint64_t)total;
+}
+static int by(const void *a, const void *b) { return (*(const double *)a > *(const double *)b) - (*(const double *)a < *(const double *)b); }
+int main(void) {
+    uint64_t x = 1;
+    for (int i = 0; i < N * N; ++i) {
+        x = x * 48271 % 2147483647, A[i] = (uint32_t)(x % 100), DA[i] = A[i];
+        x = x * 48271 % 2147483647, B[i] = (uint32_t)(x % 100), DB[i] = B[i];
+    }
+    static uint64_t (*const loops[])(void) = {index, plain, pointer, apart, checks, bits, c};
+    static const char *const names[] = {"index", "plain", "pointer", "apart", "checks", "bits", "c"};
+    double ms[7][11];
+    uint64_t answer[7];
+    for (int r = 0; r < 11; ++r)
+        for (int l = 0; l < 7; ++l) {
+            double t = now();
+            answer[l] = loops[l]();
+            ms[l][r] = now() - t;
+        }
+    printf("{");
+    for (int l = 0; l < 7; ++l) {
+        qsort(ms[l], 11, sizeof(double), by);
+        printf("%s\"%s\":%.2f", l ? "," : "", names[l], ms[l][5]);
+    }
+    int same = 1;
+    for (int l = 1; l < 7; ++l) same = same && answer[l] == answer[0];
+    printf(",\"same\":%s}\n", same ? "true" : "false");
+}