aboutsummaryrefslogtreecommitdiffstats
path: root/research/experiments
diff options
context:
space:
mode:
Diffstat (limited to 'research/experiments')
-rw-r--r--research/experiments/042-tiles/2026-10-08-tiles.jsonl54
-rw-r--r--research/experiments/042-tiles/make.py355
-rw-r--r--research/experiments/042-tiles/make2.py385
3 files changed, 794 insertions, 0 deletions
diff --git a/research/experiments/042-tiles/2026-10-08-tiles.jsonl b/research/experiments/042-tiles/2026-10-08-tiles.jsonl
new file mode 100644
index 0000000..cb32e72
--- /dev/null
+++ b/research/experiments/042-tiles/2026-10-08-tiles.jsonl
@@ -0,0 +1,54 @@
+{"binary":"tail-2-4-clang","body":"knapsack","K":2,"S":4,"T":4,"tiles_ms":29.523,"passes_ms":2.931,"c_ms":8.262,"same":true,"answer":49720}
+{"binary":"tail-2-4-clang","body":"knapsack","K":2,"S":4,"T":4,"tiles_ms":29.431,"passes_ms":2.923,"c_ms":8.252,"same":true,"answer":49720}
+{"binary":"tail-2-4-clang","body":"knapsack","K":2,"S":4,"T":4,"tiles_ms":29.435,"passes_ms":2.920,"c_ms":8.255,"same":true,"answer":49720}
+{"binary":"tail-2-4-clang","body":"roots","K":2,"S":4,"T":4,"tiles_ms":5.782,"c_ms":7.750,"same":true,"sum":21082008.97391776,"steps":1238280}
+{"binary":"tail-2-4-clang","body":"roots","K":2,"S":4,"T":4,"tiles_ms":5.780,"c_ms":7.779,"same":true,"sum":21082008.97391776,"steps":1238280}
+{"binary":"tail-2-4-clang","body":"roots","K":2,"S":4,"T":4,"tiles_ms":5.778,"c_ms":7.804,"same":true,"sum":21082008.97391776,"steps":1238280}
+{"binary":"tail-2-4-clang","body":"bfs","K":2,"S":4,"T":4,"tiles_ms":18.683,"passes_ms":9.881,"c_ms":8.618,"same":true,"reached":99998}
+{"binary":"tail-2-4-clang","body":"bfs","K":2,"S":4,"T":4,"tiles_ms":18.693,"passes_ms":9.842,"c_ms":8.650,"same":true,"reached":99998}
+{"binary":"tail-2-4-clang","body":"bfs","K":2,"S":4,"T":4,"tiles_ms":18.068,"passes_ms":9.951,"c_ms":8.874,"same":true,"reached":99998}
+{"binary":"tail-2-4-gcc","body":"knapsack","K":2,"S":4,"T":4,"tiles_ms":41.693,"passes_ms":9.591,"c_ms":6.169,"same":true,"answer":49720}
+{"binary":"tail-2-4-gcc","body":"knapsack","K":2,"S":4,"T":4,"tiles_ms":41.597,"passes_ms":9.587,"c_ms":6.167,"same":true,"answer":49720}
+{"binary":"tail-2-4-gcc","body":"knapsack","K":2,"S":4,"T":4,"tiles_ms":41.594,"passes_ms":9.590,"c_ms":6.166,"same":true,"answer":49720}
+{"binary":"tail-2-4-gcc","body":"roots","K":2,"S":4,"T":4,"tiles_ms":6.879,"c_ms":6.877,"same":true,"sum":21082008.97391776,"steps":1238280}
+{"binary":"tail-2-4-gcc","body":"roots","K":2,"S":4,"T":4,"tiles_ms":6.804,"c_ms":6.881,"same":true,"sum":21082008.97391776,"steps":1238280}
+{"binary":"tail-2-4-gcc","body":"roots","K":2,"S":4,"T":4,"tiles_ms":6.670,"c_ms":6.878,"same":true,"sum":21082008.97391776,"steps":1238280}
+{"binary":"tail-2-4-gcc","body":"bfs","K":2,"S":4,"T":4,"tiles_ms":24.460,"passes_ms":11.581,"c_ms":11.448,"same":true,"reached":99998}
+{"binary":"tail-2-4-gcc","body":"bfs","K":2,"S":4,"T":4,"tiles_ms":23.762,"passes_ms":12.632,"c_ms":11.434,"same":true,"reached":99998}
+{"binary":"tail-2-4-gcc","body":"bfs","K":2,"S":4,"T":4,"tiles_ms":23.743,"passes_ms":11.977,"c_ms":11.477,"same":true,"reached":99998}
+{"binary":"tail-1-8-clang","body":"knapsack","K":1,"S":8,"T":2,"tiles_ms":35.823,"c_ms":8.266,"same":true,"answer":49720}
+{"binary":"tail-1-8-clang","body":"knapsack","K":1,"S":8,"T":2,"tiles_ms":35.750,"c_ms":8.272,"same":true,"answer":49720}
+{"binary":"tail-1-8-clang","body":"knapsack","K":1,"S":8,"T":2,"tiles_ms":35.748,"c_ms":8.265,"same":true,"answer":49720}
+{"binary":"tail-1-8-clang","body":"roots","K":1,"S":8,"T":2,"tiles_ms":9.746,"c_ms":7.801,"same":true,"sum":21082008.973917805,"steps":1238280}
+{"binary":"tail-1-8-clang","body":"roots","K":1,"S":8,"T":2,"tiles_ms":9.735,"c_ms":7.811,"same":true,"sum":21082008.973917805,"steps":1238280}
+{"binary":"tail-1-8-clang","body":"roots","K":1,"S":8,"T":2,"tiles_ms":9.745,"c_ms":7.775,"same":true,"sum":21082008.973917805,"steps":1238280}
+{"binary":"tail-1-8-clang","body":"bfs","K":1,"S":8,"T":2,"tiles_ms":19.944,"passes_ms":10.015,"c_ms":9.460,"same":true,"reached":99998}
+{"binary":"tail-1-8-clang","body":"bfs","K":1,"S":8,"T":2,"tiles_ms":20.277,"passes_ms":9.950,"c_ms":9.441,"same":true,"reached":99998}
+{"binary":"tail-1-8-clang","body":"bfs","K":1,"S":8,"T":2,"tiles_ms":19.555,"passes_ms":9.825,"c_ms":9.677,"same":true,"reached":99998}
+{"binary":"tail-1-8-gcc","body":"knapsack","K":1,"S":8,"T":2,"tiles_ms":40.733,"c_ms":6.166,"same":true,"answer":49720}
+{"binary":"tail-1-8-gcc","body":"knapsack","K":1,"S":8,"T":2,"tiles_ms":40.717,"c_ms":6.168,"same":true,"answer":49720}
+{"binary":"tail-1-8-gcc","body":"knapsack","K":1,"S":8,"T":2,"tiles_ms":40.692,"c_ms":6.171,"same":true,"answer":49720}
+{"binary":"tail-1-8-gcc","body":"roots","K":1,"S":8,"T":2,"tiles_ms":9.942,"c_ms":6.889,"same":true,"sum":21082008.973917805,"steps":1238280}
+{"binary":"tail-1-8-gcc","body":"roots","K":1,"S":8,"T":2,"tiles_ms":9.946,"c_ms":6.875,"same":true,"sum":21082008.973917805,"steps":1238280}
+{"binary":"tail-1-8-gcc","body":"roots","K":1,"S":8,"T":2,"tiles_ms":9.944,"c_ms":6.872,"same":true,"sum":21082008.973917805,"steps":1238280}
+{"binary":"tail-1-8-gcc","body":"bfs","K":1,"S":8,"T":2,"tiles_ms":23.364,"passes_ms":11.229,"c_ms":11.616,"same":true,"reached":99998}
+{"binary":"tail-1-8-gcc","body":"bfs","K":1,"S":8,"T":2,"tiles_ms":20.582,"passes_ms":11.138,"c_ms":11.718,"same":true,"reached":99998}
+{"binary":"tail-1-8-gcc","body":"bfs","K":1,"S":8,"T":2,"tiles_ms":21.081,"passes_ms":11.145,"c_ms":11.474,"same":true,"reached":99998}
+{"binary":"tiles-8-O3 (computed goto, spilled)","body":"knapsack","T":8,"tiles_ms":46.762,"passes_ms":27.043,"c_ms":6.170,"same":true,"answer":49720}
+{"binary":"tiles-8-O3 (computed goto, spilled)","body":"knapsack","T":8,"tiles_ms":46.750,"passes_ms":27.014,"c_ms":6.168,"same":true,"answer":49720}
+{"binary":"tiles-8-O3 (computed goto, spilled)","body":"knapsack","T":8,"tiles_ms":46.650,"passes_ms":27.008,"c_ms":6.171,"same":true,"answer":49720}
+{"binary":"tiles-8-O3 (computed goto, spilled)","body":"roots","T":8,"tiles_ms":12.846,"c_ms":6.882,"same":false,"sum":21082008.973917745,"steps":1238280}
+{"binary":"tiles-8-O3 (computed goto, spilled)","body":"roots","T":8,"tiles_ms":12.900,"c_ms":6.877,"same":false,"sum":21082008.973917745,"steps":1238280}
+{"binary":"tiles-8-O3 (computed goto, spilled)","body":"roots","T":8,"tiles_ms":12.817,"c_ms":6.876,"same":false,"sum":21082008.973917745,"steps":1238280}
+{"binary":"tiles-8-O3 (computed goto, spilled)","body":"bfs","T":8,"tiles_ms":18.845,"passes_ms":11.196,"c_ms":11.413,"same":true,"reached":99998}
+{"binary":"tiles-8-O3 (computed goto, spilled)","body":"bfs","T":8,"tiles_ms":18.868,"passes_ms":11.011,"c_ms":11.387,"same":true,"reached":99998}
+{"binary":"tiles-8-O3 (computed goto, spilled)","body":"bfs","T":8,"tiles_ms":19.576,"passes_ms":11.237,"c_ms":11.400,"same":true,"reached":99998}
+{"binary":"tiles-16-O3 (computed goto, spilled)","body":"knapsack","T":16,"tiles_ms":43.138,"passes_ms":27.044,"c_ms":6.315,"same":true,"answer":49720}
+{"binary":"tiles-16-O3 (computed goto, spilled)","body":"knapsack","T":16,"tiles_ms":43.138,"passes_ms":27.025,"c_ms":6.313,"same":true,"answer":49720}
+{"binary":"tiles-16-O3 (computed goto, spilled)","body":"knapsack","T":16,"tiles_ms":43.179,"passes_ms":27.040,"c_ms":6.316,"same":true,"answer":49720}
+{"binary":"tiles-16-O3 (computed goto, spilled)","body":"roots","T":16,"tiles_ms":11.597,"c_ms":6.887,"same":false,"sum":21082008.973917749,"steps":1238280}
+{"binary":"tiles-16-O3 (computed goto, spilled)","body":"roots","T":16,"tiles_ms":11.565,"c_ms":6.877,"same":false,"sum":21082008.973917749,"steps":1238280}
+{"binary":"tiles-16-O3 (computed goto, spilled)","body":"roots","T":16,"tiles_ms":11.645,"c_ms":6.882,"same":false,"sum":21082008.973917749,"steps":1238280}
+{"binary":"tiles-16-O3 (computed goto, spilled)","body":"bfs","T":16,"tiles_ms":16.199,"passes_ms":11.197,"c_ms":11.365,"same":true,"reached":99998}
+{"binary":"tiles-16-O3 (computed goto, spilled)","body":"bfs","T":16,"tiles_ms":15.623,"passes_ms":11.305,"c_ms":11.387,"same":true,"reached":99998}
+{"binary":"tiles-16-O3 (computed goto, spilled)","body":"bfs","T":16,"tiles_ms":15.889,"passes_ms":11.203,"c_ms":11.597,"same":true,"reached":99998}
diff --git a/research/experiments/042-tiles/make.py b/research/experiments/042-tiles/make.py
new file mode 100644
index 0000000..060c5d8
--- /dev/null
+++ b/research/experiments/042-tiles/make.py
@@ -0,0 +1,355 @@
+#!/usr/bin/env python3
+"""Research 132's deciding experiment: a register-tile interpreter (tiles.c).
+
+A body's instructions run over a tile of T lanes held in vector registers
+(GCC vector types), on a stack of DEPTH vector slots whose depth is known
+per instruction (each operation has one label a depth: the C compiler sees
+fixed variables, never an indexed array), dispatched as threaded code
+(computed goto). No code is made at run time. Three bodies: knapsack's
+inner loop, BFS's frontier (an edge a lane), roots' Newton steps with lanes
+leaving. Against the same bodies as one pass per operation over 1024
+lanes (the batch executor's shape) and as plain C.
+
+ python3 make.py > tiles.c
+ cc -O2 -DT=8 -o tiles-8 tiles.c -lm (T 4, 8, 16)
+ ./tiles-8 BODY ROUNDS (knapsack | bfs | roots | all)
+
+The bar (132): cycles an element within 1.3x the compiled code's on all
+three, measured with 041's icount over ROUNDS 1 and 11.
+"""
+DEPTH = 6
+ops = [] # (name, pops, pushes, body as a format with {d}, {a}, {b}, {c} slot names)
+
+
+def op(name, pops, pushes, body, needs=None):
+ ops.append((name, pops, pushes, body, needs))
+
+
+# stack: s0 .. s{DEPTH-1}; at depth d the top is s{d-1}
+op('PUSHK', 0, 1, 's{t} = (vd){{}} + ip->k;')
+op('PUSH_AFFINE', 0, 1, 's{t} = (vd){{}} + ip->k + (lanes_v + base) * ip->k2;') # k + step * (tile base + lane)
+op('PUSH_ALL', 0, 1, 's{t} = (vd)all_v;')
+op('PUSH_ROWS', 0, 1, 'for (int i = 0; i < T; ++i) s{t}[i] = (double)((const uint32_t *)ip->p)[tile * T + i];')
+op('PICK', 0, 1, 's{t} = s{pick};', needs='pick') # push a copy of slot (d-1-a)
+op('SET', 1, 0, 's{set} = s{top};', needs='set') # pop the top into slot (d-2-a)
+op('POP', 1, 0, '')
+op('ADD', 2, 1, 's{b} = s{b} + s{top};')
+op('SUB', 2, 1, 's{b} = s{b} - s{top};')
+op('MUL', 2, 1, 's{b} = s{b} * s{top};')
+op('DIV', 2, 1, 's{b} = s{b} / s{top};')
+op('ADDK', 1, 1, 's{top} = s{top} + ip->k;')
+op('SUBK', 1, 1, 's{top} = s{top} - ip->k;')
+op('MULK', 1, 1, 's{top} = s{top} * ip->k;')
+op('DIVK', 1, 1, 's{top} = s{top} / ip->k;')
+op('LT', 2, 1, 's{b} = (vd)(s{b} < s{top});')
+op('GT', 2, 1, 's{b} = (vd)(s{b} > s{top});')
+op('AND', 2, 1, 's{b} = (vd)((vm)s{b} & (vm)s{top});')
+op('NOT', 1, 1, 's{top} = (vd)(~(vm)s{top});')
+op('GATHER4', 1, 1, '{{ const uint32_t *col = ip->p; for (int i = 0; i < T; ++i) s{top}[i] = col[(size_t)s{top}[i]]; }}')
+op('GATHER8', 1, 1, '{{ const double *col = ip->p; for (int i = 0; i < T; ++i) s{top}[i] = col[(size_t)s{top}[i]]; }}')
+op('BITS', 1, 1, '{{ const uint64_t *bits = ip->p; vm m; for (int i = 0; i < T; ++i) {{ size_t k = (size_t)s{top}[i]; m[i] = -(int64_t)(bits[k >> 6] >> (k & 63) & 1); }} s{top} = (vd)m; }}')
+op('STOREM4', 3, 0, '{{ uint32_t *col = ip->p; vm m = (vm)s{top} & valid; for (int i = 0; i < T; ++i) if (m[i]) col[(size_t)s{c}[i]] = (uint32_t)s{b}[i]; }}')
+op('SEL', 3, 1, 's{c} = (vd)(((vm)s{b} & (vm)s{top}) | ((vm)s{c} & ~(vm)s{top}));') # [a, b, m] -> [m ? b : a]
+op('ACCM', 1, 1, 'acc[ip->a] += (vd)((vm)s{top} & valid & (vm)one_v);') # acc += top mask ? 1 : 0 (peek)
+op('ACC_ADD', 1, 0, 'acc[ip->a] += (vd)((vm)s{top} & valid);')
+op('COMPACT', 1, 0, '{{ vm m = (vm)s{top} & valid; for (int i = 0; i < T; ++i) if (m[i]) ctx->out[ctx->count++] = (uint32_t)(tile * T + i); }}')
+op('JANYS', 0, 0, '{{ vm m = (vm)s{pick} & valid; int64_t any = 0; for (int i = 0; i < T; ++i) any |= m[i]; if (any && ++steps < 64) {{ ip = prog + ip->a; goto *ip->label; }} steps = 0; }}', needs='pickj')
+op('TILE', 0, 0, 'if (++tile == ctx->tiles) goto halt; base += T; for (int i = 0; i < T; ++i) valid[i] = -(int64_t)(tile * T + i < ctx->count_lanes); ip = prog; goto *ip->label;')
+
+print('/* Made by make.py (research 132): a register-tile interpreter. See make.py. */')
+print('#define _POSIX_C_SOURCE 200809L')
+print('#include <stdint.h>\n#include <stdbool.h>\n#include <stdio.h>\n#include <stdlib.h>\n#include <string.h>\n#include <math.h>\n#include <time.h>')
+print('#ifndef T\n#define T 8\n#endif')
+print(f'enum {{ DEPTH = {DEPTH} }};')
+print('typedef double vd __attribute__((vector_size(8 * T)));')
+print('typedef int64_t vm __attribute__((vector_size(8 * T)));')
+print('typedef struct { const void *label; int op, a; double k, k2; void *p; } Insn;')
+print('typedef struct { size_t tiles, count_lanes; uint32_t *out; size_t count; vd acc[2]; } Ctx;')
+print('static vd lanes_v, one_v; static vm all_v;')
+print('static void init_vectors(void) { for (int i = 0; i < T; ++i) lanes_v[i] = i, one_v[i] = 1, all_v[i] = -1; }')
+# op numbering: (name, depth)
+names = {}
+labels = []
+for name, pops, pushes, body, needs in ops:
+ for d in range(DEPTH + 1):
+ if d < pops or d - pops + pushes > DEPTH:
+ continue
+ if needs == 'pick' or needs == 'pickj':
+ for a in range(d):
+ labels.append((name, d, a))
+ elif needs == 'set':
+ for a in range(d - 1):
+ labels.append((name, d, a))
+ else:
+ labels.append((name, d, -1))
+print('enum {')
+for i, (name, d, a) in enumerate(labels):
+ print(f' OP_{name}_{d}' + (f'_{a}' if a >= 0 else '') + f' = {i},')
+print(' OPS\n};')
+print('static const void *table[OPS];')
+print('''
+/* The program run tile by tile from ctx (tiles, the valid lanes of the last), until TILE finds none left. */
+static __attribute__((noinline)) void run(const Insn *prog, Ctx *ctx) {
+ if (!prog) {''')
+for i, (name, d, a) in enumerate(labels):
+ print(f' table[{i}] = &&L_{name}_{d}' + (f'_{a}' if a >= 0 else '') + ';')
+print(''' return;
+ }
+ register const Insn *ip = prog;
+ vd ''' + ', '.join(f's{i} = {{0}}' for i in range(DEPTH)) + ''';
+ vd acc[2] = {ctx->acc[0], ctx->acc[1]};
+ vm valid;
+ double base = 0;
+ size_t tile = 0;
+ unsigned steps = 0;
+ for (int i = 0; i < T; ++i) valid[i] = -(int64_t)((size_t)i < ctx->count_lanes);
+ goto *ip->label;''')
+for name, d, a in labels:
+ spec = next(o for o in ops if o[0] == name)
+ _, pops, pushes, body, needs = spec
+ top = f'{d-1}'
+ b = f'{d-2}'
+ c = f'{d-3}'
+ t = f'{d}'
+ pick = f'{d-1-a}' if a >= 0 and needs in ('pick', 'pickj') else ''
+ sett = f'{d-2-a}' if a >= 0 and needs == 'set' else ''
+ code = body.format(t=t, top=top, b=b, c=c, pick=pick, set=sett)
+ suffix = f'_{a}' if a >= 0 else ''
+ print(f' L_{name}_{d}{suffix}: {code} ++ip; goto *ip->label;')
+print(''' halt:
+ ctx->acc[0] = acc[0], ctx->acc[1] = acc[1];
+}
+
+static Insn *at;
+static Insn *emit(int opc, int a, double k, double k2, void *p) {
+ Insn *i = at++;
+ i->label = table[opc], i->op = opc, i->a = a, i->k = k, i->k2 = k2, i->p = p;
+ return i;
+}
+#define I(name, d) emit(OP_##name##_##d, 0, 0, 0, NULL)
+#define IK(name, d, k) emit(OP_##name##_##d, 0, k, 0, NULL)
+#define IP(name, d, p) emit(OP_##name##_##d, 0, 0, 0, p)
+#define IA(name, d, a) emit(OP_##name##_##d##_##a, a, 0, 0, NULL)
+
+static double now_ms(void) { struct timespec t; clock_gettime(CLOCK_MONOTONIC, &t); return t.tv_sec * 1e3 + t.tv_nsec / 1e6; }
+static uint64_t minstd(uint64_t x) { return x * 48271 % 2147483647; }
+
+/* ---- knapsack: 100 items, capacity 50,000 (the suite's), the table in u32; the inner loop the body */
+static double knapsack_tiles(unsigned n, unsigned capacity, const uint32_t *weight, const uint32_t *worth, uint32_t *best, double *ms) {
+ Insn prog[16];
+ double a = now_ms();
+ for (unsigned i = 0; i < n; ++i) {
+ uint32_t w = weight[i];
+ if (w > capacity) continue;
+ at = prog;
+ emit(OP_PUSH_AFFINE_0, 0, (double)capacity, -1, NULL); /* c = capacity - k */
+ IA(PICK, 1, 0); IK(SUBK, 2, (double)w); IP(GATHER4, 2, best); IK(ADDK, 2, (double)worth[i]); /* [c, with] */
+ IA(PICK, 2, 0); IA(PICK, 3, 2); IP(GATHER4, 4, best); I(GT, 4); IP(STOREM4, 3, best); I(TILE, 0); /* [c, with, with > best[c]] */
+ size_t count = capacity + 1 - w;
+ Ctx ctx = {.tiles = (count + T - 1) / T, .count_lanes = count};
+ run(prog, &ctx);
+ }
+ *ms = now_ms() - a;
+ return best[capacity];
+}
+
+static double knapsack_passes(unsigned n, unsigned capacity, const uint32_t *weight, const uint32_t *worth, uint32_t *best, double *ms) {
+ enum { LANES = 1024 };
+ static double c_[LANES], with[LANES], old[LANES];
+ static uint32_t mask[LANES];
+ double a = now_ms();
+ for (unsigned i = 0; i < n; ++i) {
+ uint32_t w = weight[i];
+ if (w > capacity) continue;
+ size_t count = capacity + 1 - w;
+ for (size_t k0 = 0; k0 < count; k0 += LANES) {
+ size_t m = count - k0 < LANES ? count - k0 : LANES;
+ for (size_t l = 0; l < m; ++l) c_[l] = (double)capacity - (double)(k0 + l);
+ for (size_t l = 0; l < m; ++l) with[l] = best[(size_t)(c_[l] - w)];
+ for (size_t l = 0; l < m; ++l) with[l] += worth[i];
+ for (size_t l = 0; l < m; ++l) old[l] = best[(size_t)c_[l]];
+ for (size_t l = 0; l < m; ++l) mask[l] = with[l] > old[l];
+ for (size_t l = 0; l < m; ++l) if (mask[l]) best[(size_t)c_[l]] = (uint32_t)with[l];
+ }
+ }
+ *ms = now_ms() - a;
+ return best[capacity];
+}
+
+static double knapsack_c(unsigned n, unsigned capacity, const uint32_t *weight, const uint32_t *worth, uint32_t *best, double *ms) {
+ double a = now_ms();
+ 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;
+ }
+ *ms = now_ms() - a;
+ return best[capacity];
+}
+
+/* ---- roots: the suite's 100,000 Newton square roots, a tile of lanes stepping together, each leaving as it converges */
+static void roots_tiles(unsigned n, double *sum, double *steps, double *ms) {
+ Insn prog[32];
+ at = prog;
+ emit(OP_PUSH_AFFINE_0, 0, 1, 1, NULL); /* [i] */
+ I(PUSH_ALL, 1); /* [i, live] */
+ IA(PICK, 2, 1); /* [i, live, r] */
+ Insn *loop = at;
+ IA(PICK, 3, 2); IA(PICK, 4, 1); I(DIV, 5); IA(PICK, 4, 1); I(ADD, 5); IK(DIVK, 4, 2); /* [i, live, r, next] */
+ IA(PICK, 4, 0); IA(PICK, 5, 2); I(LT, 6); IA(PICK, 5, 3); I(AND, 6); /* [i, live, r, next, live'] */
+ IA(SET, 5, 2); /* [i, live', r, next] */
+ IA(PICK, 4, 2); emit(OP_ACCM_5, 1, 0, 0, NULL); I(SEL, 5); /* count += live'; [i, live', r'] */
+ emit(OP_JANYS_3_1, 1, 0, 0, NULL)->a = (int)(loop - prog);
+ emit(OP_ACC_ADD_3, 0, 0, 0, NULL); /* sum += r */
+ I(POP, 2); I(POP, 1); I(TILE, 0);
+ double a = now_ms();
+ Ctx ctx = {.tiles = (n + T - 1) / T, .count_lanes = n};
+ run(prog, &ctx);
+ *ms = now_ms() - a;
+ *sum = *steps = 0;
+ for (int i = 0; i < T; ++i) *sum += ctx.acc[0][i], *steps += ctx.acc[1][i];
+}
+
+static void roots_c(unsigned n, double *sum, double *steps, double *ms) {
+ double a = now_ms(), s = 0, c = 0;
+ for (uint32_t i = 1; i <= n; ++i) {
+ double r = i;
+ for (int k = 0; k < 64; ++k) {
+ double next = (r + i / r) / 2;
+ if (next >= r) break;
+ r = next;
+ ++c;
+ }
+ s += r;
+ }
+ *sum = s, *steps = c, *ms = now_ms() - a;
+}
+
+/* ---- bfs: the suite's graph (100,000 nodes, 1,000,000 edges drawn by xorshift from 121), a search from a node; the
+ * frontier's edges spread a lane each: to gathered, seen's bits asked, the candidates compacted; then the tail in C. */
+static uint64_t gstate = 121;
+static uint64_t gnext(void) { gstate ^= gstate << 13, gstate ^= gstate >> 7, gstate ^= gstate << 17; return gstate; }
+static int cmp_u64(const void *a, const void *b) { uint64_t x = *(const uint64_t *)a, y = *(const uint64_t *)b; return (x > y) - (x < y); }
+typedef struct { uint32_t nodes; uint32_t *start, *to; } Graph;
+static Graph graph_build(uint32_t nodes, size_t edges) {
+ uint64_t *keys = malloc(edges * 8);
+ size_t c = 0;
+ gstate = 121;
+ while (c < edges) { uint64_t a = gnext() % nodes, b = gnext() % nodes; if (a != b) keys[c++] = a << 20 | b; }
+ qsort(keys, c, 8, cmp_u64);
+ size_t u = 0;
+ for (size_t i = 0; i < c; ++i) if (!u || keys[i] != keys[u - 1]) keys[u++] = keys[i];
+ Graph g = {nodes, calloc(nodes + 1, 4), malloc(4 * u)};
+ for (size_t i = 0; i < u; ++i) ++g.start[(keys[i] >> 20) + 1], g.to[i] = keys[i] & 0xFFFFF;
+ for (uint32_t v = 0; v < nodes; ++v) g.start[v + 1] += g.start[v];
+ free(keys);
+ return g;
+}
+
+enum { BATCH = 1024, LANES = 65536 };
+static uint32_t rows[LANES], owner[LANES], cand[LANES];
+
+/* One search: the queue walked a batch of nodes at a time, the batch's edges spread into lanes (owner: the node). tiles: the
+ * frontier by the tile program; else by one pass an operation over the lanes (the batch executor's shape). */
+static uint32_t bfs_search(const Graph *g, uint32_t from, bool tiles, uint64_t *seen, uint32_t *depth, uint32_t *queue, double *ms) {
+ Insn prog[8];
+ at = prog;
+ IP(PUSH_ROWS, 0, rows); IP(GATHER4, 1, g->to); IA(PICK, 1, 0); IP(BITS, 2, seen); I(NOT, 2); I(COMPACT, 2); I(TILE, 0);
+ double a = now_ms();
+ memset(seen, 0, (g->nodes + 63) / 64 * 8);
+ size_t count = 1, i = 0;
+ queue[0] = from, seen[from >> 6] |= UINT64_C(1) << (from & 63), depth[from] = 0;
+ while (i < count) {
+ size_t end = count < i + BATCH ? count : i + BATCH, n = 0;
+ for (size_t r = i; r < end; ++r) { /* spread: an edge a lane (a batch's edges fit: 1024 nodes x 10 on average) */
+ uint32_t v = queue[r];
+ for (uint32_t e = g->start[v]; e < g->start[v + 1]; ++e) {
+ if (n == LANES) { fputs("too many edges in a batch\\n", stderr); exit(1); }
+ rows[n] = e, owner[n++] = v;
+ }
+ }
+ Ctx ctx = {.tiles = (n + T - 1) / T, .count_lanes = n, .out = cand, .count = 0};
+ if (tiles) run(prog, &ctx);
+ else {
+ static uint32_t to_[LANES], has[LANES];
+ for (size_t l = 0; l < n; ++l) to_[l] = g->to[rows[l]];
+ for (size_t l = 0; l < n; ++l) has[l] = seen[to_[l] >> 6] >> (to_[l] & 63) & 1;
+ for (size_t l = 0; l < n; ++l) if (!has[l]) cand[ctx.count++] = (uint32_t)l;
+ }
+ for (size_t k = 0; k < ctx.count; ++k) { /* the tail: each candidate asked again in order, added, appended */
+ uint32_t l = cand[k], t = g->to[rows[l]];
+ if (seen[t >> 6] >> (t & 63) & 1) continue;
+ seen[t >> 6] |= UINT64_C(1) << (t & 63);
+ depth[t] = depth[owner[l]] + 1;
+ queue[count++] = t;
+ }
+ i = end;
+ }
+ *ms = now_ms() - a;
+ return (uint32_t)count;
+}
+
+static uint32_t bfs_c(const Graph *g, uint32_t from, uint64_t *seen, uint32_t *depth, uint32_t *queue, double *ms) {
+ double a = now_ms();
+ memset(seen, 0, (g->nodes + 63) / 64 * 8);
+ size_t count = 1;
+ queue[0] = from, seen[from >> 6] |= UINT64_C(1) << (from & 63), depth[from] = 0;
+ for (size_t i = 0; i < count; ++i) {
+ uint32_t v = queue[i];
+ for (uint32_t e = g->start[v]; e < g->start[v + 1]; ++e) {
+ uint32_t t = g->to[e];
+ if (seen[t >> 6] >> (t & 63) & 1) continue;
+ seen[t >> 6] |= UINT64_C(1) << (t & 63);
+ depth[t] = depth[v] + 1;
+ queue[count++] = t;
+ }
+ }
+ *ms = now_ms() - a;
+ return (uint32_t)count;
+}
+
+int main(int argc, char **argv) {
+ const char *body = argc > 1 ? argv[1] : "all";
+ int rounds = argc > 2 ? atoi(argv[2]) : 1;
+ init_vectors();
+ run(NULL, NULL);
+ bool all = !strcmp(body, "all");
+ if (all || !strcmp(body, "knapsack")) {
+ unsigned n = 100, capacity = 50000;
+ uint32_t *weight = malloc(4 * n), *worth = malloc(4 * n), *best = calloc(capacity + 1, 4);
+ uint64_t x = 1;
+ for (unsigned i = 0; i < n; ++i) x = minstd(x), weight[i] = (uint32_t)(x % 1000) + 1, x = minstd(x), worth[i] = (uint32_t)(x % 1000);
+ for (int r = 0; r < rounds; ++r) {
+ double ms[3], answers[3];
+ memset(best, 0, 4 * (capacity + 1)), answers[0] = knapsack_tiles(n, capacity, weight, worth, best, &ms[0]);
+ memset(best, 0, 4 * (capacity + 1)), answers[1] = knapsack_passes(n, capacity, weight, worth, best, &ms[1]);
+ memset(best, 0, 4 * (capacity + 1)), answers[2] = knapsack_c(n, capacity, weight, worth, best, &ms[2]);
+ printf("{\\"body\\":\\"knapsack\\",\\"T\\":%d,\\"tiles_ms\\":%.3f,\\"passes_ms\\":%.3f,\\"c_ms\\":%.3f,\\"same\\":%s,\\"answer\\":%.17g}\\n", T, ms[0], ms[1], ms[2],
+ answers[0] == answers[2] && answers[1] == answers[2] ? "true" : "false", answers[0]);
+ }
+ }
+ if (all || !strcmp(body, "roots")) {
+ for (int r = 0; r < rounds; ++r) {
+ double sum[2], steps[2], ms[2];
+ roots_tiles(100000, &sum[0], &steps[0], &ms[0]);
+ roots_c(100000, &sum[1], &steps[1], &ms[1]);
+ printf("{\\"body\\":\\"roots\\",\\"T\\":%d,\\"tiles_ms\\":%.3f,\\"c_ms\\":%.3f,\\"same\\":%s,\\"sum\\":%.17g,\\"steps\\":%.0f}\\n", T, ms[0], ms[1],
+ sum[0] == sum[1] && steps[0] == steps[1] ? "true" : "false", sum[0], steps[0]);
+ }
+ }
+ if (all || !strcmp(body, "bfs")) {
+ Graph g = graph_build(100000, 1000000);
+ uint64_t *seen = calloc((g.nodes + 63) / 64, 8);
+ uint32_t *depth = calloc(g.nodes, 4), *queue = malloc(4 * g.nodes);
+ for (int r = 0; r < rounds; ++r) {
+ uint32_t from = (uint32_t)(gnext() % g.nodes), reached[3];
+ double ms[3];
+ reached[0] = bfs_search(&g, from, true, seen, depth, queue, &ms[0]);
+ reached[1] = bfs_search(&g, from, false, seen, depth, queue, &ms[1]);
+ reached[2] = bfs_c(&g, from, seen, depth, queue, &ms[2]);
+ printf("{\\"body\\":\\"bfs\\",\\"T\\":%d,\\"tiles_ms\\":%.3f,\\"passes_ms\\":%.3f,\\"c_ms\\":%.3f,\\"same\\":%s,\\"reached\\":%u}\\n", T, ms[0], ms[1], ms[2],
+ reached[0] == reached[2] && reached[1] == reached[2] ? "true" : "false", reached[0]);
+ }
+ }
+ return 0;
+}''')
diff --git a/research/experiments/042-tiles/make2.py b/research/experiments/042-tiles/make2.py
new file mode 100644
index 0000000..5278d1b
--- /dev/null
+++ b/research/experiments/042-tiles/make2.py
@@ -0,0 +1,385 @@
+#!/usr/bin/env python3
+"""Research 132's deciding experiment, second form: a register-tile
+interpreter as threaded tail calls (tail.c).
+
+make.py's form (computed goto, one function) spilled every slot to the
+stack: GCC keeps nothing in registers across an indirect branch. Here each
+operation is a function whose arguments are the slots (homogeneous vector
+aggregates: AAPCS64 passes them in v0-v7) and which tail-calls the next
+(clang's musttail; GCC's sibling calls), so the slots stay in registers:
+an operation is its vector instructions, a load of the next address and a
+branch. K registers a slot (a tile of 2K doubles), S slots (K * S <= 8).
+Operations name their slots (three-address), so each (op, slots) is a
+function; only the ones the three programs use are made.
+
+ python3 make2.py K S > tail-K-S.c
+ clang -O2 -o tail-K-S tail-K-S.c -lm (or gcc)
+ ./tail-K-S BODY ROUNDS (knapsack | bfs | roots | all)
+"""
+import sys
+
+K = int(sys.argv[1]) if len(sys.argv) > 1 else 2
+S = int(sys.argv[2]) if len(sys.argv) > 2 else 4
+T = 2 * K
+assert K * S <= 8
+F = 'abcd'[:K] # a slot's fields, a register each
+
+# ---- programs: (op, dst, a, b, k, k2, p, target); slots by number; None where not used
+programs = {}
+# knapsack: slots 0 c, 1 with, 2 best then the mask
+programs['knapsack'] = [
+ ('AFFINE', 0, None, None, 'capacity', -1.0, None, None),
+ ('SUBK', 1, 0, None, 'w', None, None, None),
+ ('GATHER4', 1, 1, None, None, None, 'best', None),
+ ('ADDK', 1, 1, None, 'worth', None, None, None),
+ ('GATHER4', 2, 0, None, None, None, 'best', None),
+ ('GT', 2, 1, 2, None, None, None, None),
+ ('STOREM4', None, 0, 1, None, None, 'best', 2), # rows a, value b, mask target
+ ('TILE', None, None, None, None, None, None, None),
+]
+# bfs: slots 0 row then to, 1 has then its not
+programs['bfs'] = [
+ ('ROWS', 0, None, None, None, None, 'rows', None),
+ ('GATHER4', 0, 0, None, None, None, 'to', None),
+ ('BITS', 1, 0, None, None, None, 'seen', None),
+ ('NOT', 1, 1, None, None, None, None, None),
+ ('COMPACT', None, 1, None, None, None, 'cand', None),
+ ('TILE', None, None, None, None, None, None, None),
+]
+if S >= 5: # roots: 0 i, 1 live, 2 r, 3 next, 4 a mask
+ programs['roots'] = [
+ ('AFFINE', 0, None, None, 1.0, 1.0, None, None),
+ ('ALL', 1, None, None, None, None, None, None),
+ ('MOV', 2, 0, None, None, None, None, None),
+ ('DIV', 3, 0, 2, None, None, None, None), # 3: loop
+ ('ADD', 3, 3, 2, None, None, None, None),
+ ('DIVK', 3, 3, None, 2.0, None, None, None),
+ ('LT', 4, 3, 2, None, None, None, None),
+ ('AND', 1, 1, 4, None, None, None, None),
+ ('ACCM', None, 1, None, None, None, None, 1), # acc[1] += live
+ ('SEL', 2, 3, 1, None, None, None, None), # r = live ? next : r
+ ('JANY', None, 1, None, None, None, None, 3),
+ ('ACC_ADD', None, 2, None, None, None, None, 0),
+ ('TILE', None, None, None, None, None, None, None),
+ ]
+else: # roots in four slots: the compare and the and as one (CMPLT_AND: dst &= a < b), the only fused operation here
+ programs['roots'] = [
+ ('AFFINE', 0, None, None, 1.0, 1.0, None, None),
+ ('ALL', 1, None, None, None, None, None, None),
+ ('MOV', 2, 0, None, None, None, None, None),
+ ('DIV', 3, 0, 2, None, None, None, None), # 3: loop
+ ('ADD', 3, 3, 2, None, None, None, None),
+ ('DIVK', 3, 3, None, 2.0, None, None, None),
+ ('CMPLT_AND', 1, 3, 2, None, None, None, None),
+ ('ACCM', None, 1, None, None, None, None, 1),
+ ('SEL', 2, 3, 1, None, None, None, None),
+ ('JANY', None, 1, None, None, None, None, 3),
+ ('ACC_ADD', None, 2, None, None, None, None, 0),
+ ('TILE', None, None, None, None, None, None, None),
+ ]
+
+# ---- operations: C for the fields, with {d} {a} {b} slot names and {f} the field
+def per_field(fmt):
+ return ' '.join(fmt.replace('{f}', f) for f in F)
+
+bodies = {
+ 'AFFINE': lambda d, a, b: per_field('s{d}.{f} = vfmaq_f64(vdupq_n_f64(ip->k), vaddq_f64(lanes_{f}, vdupq_n_f64(c->base)), vdupq_n_f64(ip->k2));'.replace('{d}', str(d))),
+ 'ALL': lambda d, a, b: per_field('s{d}.{f} = all_f64;'.replace('{d}', str(d))),
+ 'MOV': lambda d, a, b: per_field('s{d}.{f} = s{a}.{f};'.replace('{d}', str(d)).replace('{a}', str(a))),
+ 'ADD': lambda d, a, b: per_field('s{d}.{f} = vaddq_f64(s{a}.{f}, s{b}.{f});'.replace('{d}', str(d)).replace('{a}', str(a)).replace('{b}', str(b))),
+ 'SUB': lambda d, a, b: per_field('s{d}.{f} = vsubq_f64(s{a}.{f}, s{b}.{f});'.replace('{d}', str(d)).replace('{a}', str(a)).replace('{b}', str(b))),
+ 'MUL': lambda d, a, b: per_field('s{d}.{f} = vmulq_f64(s{a}.{f}, s{b}.{f});'.replace('{d}', str(d)).replace('{a}', str(a)).replace('{b}', str(b))),
+ 'DIV': lambda d, a, b: per_field('s{d}.{f} = vdivq_f64(s{a}.{f}, s{b}.{f});'.replace('{d}', str(d)).replace('{a}', str(a)).replace('{b}', str(b))),
+ 'ADDK': lambda d, a, b: per_field('s{d}.{f} = vaddq_f64(s{a}.{f}, vdupq_n_f64(ip->k));'.replace('{d}', str(d)).replace('{a}', str(a))),
+ 'SUBK': lambda d, a, b: per_field('s{d}.{f} = vsubq_f64(s{a}.{f}, vdupq_n_f64(ip->k));'.replace('{d}', str(d)).replace('{a}', str(a))),
+ 'DIVK': lambda d, a, b: per_field('s{d}.{f} = vdivq_f64(s{a}.{f}, vdupq_n_f64(ip->k));'.replace('{d}', str(d)).replace('{a}', str(a))),
+ 'LT': lambda d, a, b: per_field('s{d}.{f} = vreinterpretq_f64_u64(vcltq_f64(s{a}.{f}, s{b}.{f}));'.replace('{d}', str(d)).replace('{a}', str(a)).replace('{b}', str(b))),
+ 'GT': lambda d, a, b: per_field('s{d}.{f} = vreinterpretq_f64_u64(vcgtq_f64(s{a}.{f}, s{b}.{f}));'.replace('{d}', str(d)).replace('{a}', str(a)).replace('{b}', str(b))),
+ 'CMPLT_AND': lambda d, a, b: per_field('s{d}.{f} = vreinterpretq_f64_u64(vandq_u64(vreinterpretq_u64_f64(s{d}.{f}), vcltq_f64(s{a}.{f}, s{b}.{f})));'.replace('{d}', str(d)).replace('{a}', str(a)).replace('{b}', str(b))),
+ 'AND': lambda d, a, b: per_field('s{d}.{f} = vreinterpretq_f64_u64(vandq_u64(vreinterpretq_u64_f64(s{a}.{f}), vreinterpretq_u64_f64(s{b}.{f})));'.replace('{d}', str(d)).replace('{a}', str(a)).replace('{b}', str(b))),
+ 'NOT': lambda d, a, b: per_field('s{d}.{f} = vreinterpretq_f64_u64(veorq_u64(vreinterpretq_u64_f64(s{a}.{f}), all_u64));'.replace('{d}', str(d)).replace('{a}', str(a))),
+ 'SEL': lambda d, a, b: per_field('s{d}.{f} = vbslq_f64(vreinterpretq_u64_f64(s{b}.{f}), s{a}.{f}, s{d}.{f});'.replace('{d}', str(d)).replace('{a}', str(a)).replace('{b}', str(b))),
+ 'GATHER4': lambda d, a, b: '{ const uint32_t *col = ip->p; ' + per_field('s{d}.{f} = vsetq_lane_f64(col[(size_t)vgetq_lane_f64(s{a}.{f}, 0)], s{d}.{f}, 0); s{d}.{f} = vsetq_lane_f64(col[(size_t)vgetq_lane_f64(s{a}.{f}, 1)], s{d}.{f}, 1);'.replace('{d}', str(d)).replace('{a}', str(a))) + ' }',
+ 'BITS': lambda d, a, b: '{ const uint64_t *bits = ip->p; uint64x2_t m = vdupq_n_u64(0); size_t k; ' + per_field('k = (size_t)vgetq_lane_f64(s{a}.{f}, 0); m = vsetq_lane_u64(-(bits[k >> 6] >> (k & 63) & 1), m, 0); k = (size_t)vgetq_lane_f64(s{a}.{f}, 1); m = vsetq_lane_u64(-(bits[k >> 6] >> (k & 63) & 1), m, 1); s{d}.{f} = vreinterpretq_f64_u64(m);'.replace('{d}', str(d)).replace('{a}', str(a))) + ' }',
+ 'ROWS': lambda d, a, b: '{ const uint32_t *rows = (const uint32_t *)ip->p + c->tile * T; ' + ' '.join(f's{d}.{f} = vcvtq_f64_u64(vmovl_u32(vld1_u32(rows + {2 * i})));' for i, f in enumerate(F)) + ' }',
+ 'STOREM4': lambda d, a, b: '{ uint32_t *col = ip->p; uint64x2_t m; ' + ' '.join(f'm = vandq_u64(vreinterpretq_u64_f64(s{{m}}.{f}), c->valid.{f}); if (vgetq_lane_u64(m, 0)) col[(size_t)vgetq_lane_f64(s{a}.{f}, 0)] = (uint32_t)vgetq_lane_f64(s{b}.{f}, 0); if (vgetq_lane_u64(m, 1)) col[(size_t)vgetq_lane_f64(s{a}.{f}, 1)] = (uint32_t)vgetq_lane_f64(s{b}.{f}, 1);' for f in F) + ' }',
+ 'COMPACT': lambda d, a, b: '{ uint32_t *out = ip->p; uint64x2_t m; ' + ' '.join(f'm = vandq_u64(vreinterpretq_u64_f64(s{a}.{f}), c->valid.{f}); if (vgetq_lane_u64(m, 0)) out[c->count++] = (uint32_t)(c->tile * T + {2 * i}); if (vgetq_lane_u64(m, 1)) out[c->count++] = (uint32_t)(c->tile * T + {2 * i + 1});' for i, f in enumerate(F)) + ' }',
+ 'ACCM': lambda d, a, b: per_field('c->acc[{n}].{f} = vaddq_f64(c->acc[{n}].{f}, vbslq_f64(vandq_u64(vreinterpretq_u64_f64(s{a}.{f}), c->valid.{f}), one_f64, zero_f64));'.replace('{a}', str(a)).replace('{n}', '{n}')),
+ 'ACC_ADD': lambda d, a, b: per_field('c->acc[{n}].{f} = vaddq_f64(c->acc[{n}].{f}, vbslq_f64(c->valid.{f}, s{a}.{f}, zero_f64));'.replace('{a}', str(a)).replace('{n}', '{n}')),
+ 'JANY': lambda d, a, b: '{ uint64x2_t m = ' + ' | '.join(f'vandq_u64(vreinterpretq_u64_f64(s{a}.{f}), c->valid.{f})' for f in F) + '; if ((vgetq_lane_u64(m, 0) | vgetq_lane_u64(m, 1)) && ++c->steps < 64) { ip = c->prog + ip->a - 1; } else c->steps = 0; }',
+ 'TILE': lambda d, a, b: '{ if (++c->tile == c->tiles) return; c->base += T; ' + ' '.join(f'c->valid.{f} = vsetq_lane_u64(-(uint64_t)(c->tile * T + {2 * i} < c->count_lanes), c->valid.{f}, 0); c->valid.{f} = vsetq_lane_u64(-(uint64_t)(c->tile * T + {2 * i + 1} < c->count_lanes), c->valid.{f}, 1);' for i, f in enumerate(F)) + ' ip = c->prog - 1; }',
+}
+
+used = {}
+for name, prog in programs.items():
+ for op, d, a, b, k, k2, p, target in prog:
+ n = target if op in ('ACCM', 'ACC_ADD', 'STOREM4') else None
+ used[(op, d, a, b, n)] = True
+
+def fname(op, d, a, b, n):
+ return 'op_' + op.lower() + ''.join('_' + str(x) for x in (d, a, b, n) if x is not None)
+
+args = ', '.join(f'slot s{i}' for i in range(S))
+pass_args = ', '.join(f's{i}' for i in range(S))
+print(f'/* Made by make2.py {K} {S} (research 132): a register-tile interpreter as threaded tail calls. See make2.py. */')
+print('''#define _POSIX_C_SOURCE 200809L
+#include <arm_neon.h>
+#include <stddef.h>
+#include <stdint.h>
+#include <stdbool.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <time.h>
+#ifdef __clang__
+#define MUSTTAIL __attribute__((musttail))
+#else
+#define MUSTTAIL
+#endif''')
+print(f'enum {{ K = {K}, S = {S}, T = {T} }};')
+print('typedef struct { ' + ' '.join(f'float64x2_t {f};' for f in F) + ' } slot;')
+print('typedef struct { ' + ' '.join(f'uint64x2_t {f};' for f in F) + ' } vmask;')
+print('typedef struct insn insn;')
+print('typedef struct { double base; size_t tile, tiles, count_lanes, count; unsigned steps; vmask valid; slot acc[2]; const insn *prog; uint32_t *out; } Ctx;')
+print(f'typedef void op_f(const insn *ip, Ctx *c, {args});')
+print('struct insn { op_f *f; double k, k2; void *p; int a; };')
+print('static float64x2_t ' + ', '.join(f'lanes_{f}' for f in F) + ', all_f64, one_f64, zero_f64; static uint64x2_t all_u64;')
+print('static void init_vectors(void) { ' + ' '.join(f'lanes_{f} = (float64x2_t){{{2 * i}, {2 * i + 1}}};' for i, f in enumerate(F)) + ' all_u64 = vdupq_n_u64(~UINT64_C(0)); all_f64 = vreinterpretq_f64_u64(all_u64); one_f64 = vdupq_n_f64(1); zero_f64 = vdupq_n_f64(0); }')
+print(f'#define NEXT do {{ const insn *n = ip + 1; MUSTTAIL return n->f(n, c, {pass_args}); }} while (0)')
+for (op, d, a, b, n) in used:
+ body = bodies[op](d, a, b)
+ if n is not None:
+ body = body.replace('{n}', str(n)).replace('{m}', str(n))
+ print(f'static void {fname(op, d, a, b, n)}(const insn *ip, Ctx *c, {args}) {{ {body} NEXT; }}')
+print('''
+static insn *at;
+static insn *emit(op_f *f, double k, double k2, void *p, int a) {
+ insn *i = at++;
+ i->f = f, i->k = k, i->k2 = k2, i->p = p, i->a = a;
+ return i;
+}
+static void start(const insn *prog, Ctx *c) {
+ slot z = {''' + ', '.join('vdupq_n_f64(0)' for _ in F) + '''};
+ c->prog = prog, c->tile = 0, c->base = 0, c->steps = 0;
+ ''' + ' '.join(f'c->valid.{f} = (uint64x2_t){{-(uint64_t)({2 * i} < c->count_lanes), -(uint64_t)({2 * i + 1} < c->count_lanes)}};' for i, f in enumerate(F)) + '''
+ prog->f(prog, c, ''' + ', '.join('z' for _ in range(S)) + ''');
+}
+static double now_ms(void) { struct timespec t; clock_gettime(CLOCK_MONOTONIC, &t); return t.tv_sec * 1e3 + t.tv_nsec / 1e6; }
+static uint64_t minstd(uint64_t x) { return x * 48271 % 2147483647; }
+''')
+
+def emit_program(name, prog, vars):
+ lines = []
+ for op, d, a, b, k, k2, p, target in prog:
+ n = target if op in ('ACCM', 'ACC_ADD', 'STOREM4') else None
+ kk = vars.get(k, k) if isinstance(k, str) else (k if k is not None else 0)
+ kk2 = k2 if k2 is not None else 0
+ pp = vars.get(p, 'NULL') if p else 'NULL'
+ aa = target if op in ('JANY', 'STOREM4') else 0
+ lines.append(f' emit({fname(op, d, a, b, n)}, {kk}, {kk2}, {pp}, {aa});')
+ return '\n'.join(lines)
+
+print('''/* ---- knapsack: 100 items, capacity 50,000 (the suite's), the table in u32; the inner loop the body */
+static double knapsack_tiles(unsigned n, unsigned capacity, const uint32_t *weight, const uint32_t *worth, uint32_t *best, double *ms) {
+ insn prog[16];
+ double a = now_ms();
+ for (unsigned i = 0; i < n; ++i) {
+ uint32_t w = weight[i];
+ if (w > capacity) continue;
+ at = prog;
+ double capacity_ = capacity, w_ = w, worth_ = worth[i];
+''' + emit_program('knapsack', programs['knapsack'], {'capacity': 'capacity_', 'w': 'w_', 'worth': 'worth_', 'best': 'best'}) + '''
+ size_t count = capacity + 1 - w;
+ Ctx ctx = {.tiles = (count + T - 1) / T, .count_lanes = count};
+ start(prog, &ctx);
+ }
+ *ms = now_ms() - a;
+ return best[capacity];
+}
+
+/* The same body as the batch executor's shape done tightly: four-byte lanes, affine rows read as strides, two passes a batch
+ * of 1024 rounds (with = best[c - w] + worth; best[c] = max(best[c], with), the store where it is greater as a select). */
+static double knapsack_passes(unsigned n, unsigned capacity, const uint32_t *weight, const uint32_t *worth, uint32_t *best, double *ms) {
+ enum { LANES = 1024 };
+ static uint32_t with[LANES];
+ double a = now_ms();
+ for (unsigned i = 0; i < n; ++i) {
+ uint32_t w = weight[i], v = worth[i];
+ if (w > capacity) continue;
+ size_t count = capacity + 1 - w;
+ for (size_t k0 = 0; k0 < count; k0 += LANES) {
+ size_t m = count - k0 < LANES ? count - k0 : LANES;
+ const uint32_t *from = best + (capacity - k0 - w);
+ uint32_t *to = best + (capacity - k0);
+ for (size_t l = 0; l < m; ++l) with[l] = from[-(ptrdiff_t)l] + v;
+ for (size_t l = 0; l < m; ++l) { uint32_t cur = to[-(ptrdiff_t)l]; to[-(ptrdiff_t)l] = with[l] > cur ? with[l] : cur; }
+ }
+ }
+ *ms = now_ms() - a;
+ return best[capacity];
+}
+
+static double knapsack_c(unsigned n, unsigned capacity, const uint32_t *weight, const uint32_t *worth, uint32_t *best, double *ms) {
+ double a = now_ms();
+ 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;
+ }
+ *ms = now_ms() - a;
+ return best[capacity];
+}
+
+/* ---- roots: the suite's 100,000 Newton square roots, a tile stepping together, each lane leaving as it converges */
+static void roots_tiles(unsigned n, double *sum, double *steps, double *ms) {
+ insn prog[16];
+ at = prog;
+''' + emit_program('roots', programs['roots'], {}) + '''
+ Ctx ctx = {.tiles = (n + T - 1) / T, .count_lanes = n};
+ double a = now_ms();
+ start(prog, &ctx);
+ *ms = now_ms() - a;
+ *sum = *steps = 0;
+ ''' + ' '.join(f'*sum += vgetq_lane_f64(ctx.acc[0].{f}, 0) + vgetq_lane_f64(ctx.acc[0].{f}, 1); *steps += vgetq_lane_f64(ctx.acc[1].{f}, 0) + vgetq_lane_f64(ctx.acc[1].{f}, 1);' for f in F) + '''
+}
+
+static void roots_c(unsigned n, double *sum, double *steps, double *ms) {
+ double a = now_ms(), s = 0, c = 0;
+ for (uint32_t i = 1; i <= n; ++i) {
+ double r = i;
+ for (int k = 0; k < 64; ++k) {
+ double next = (r + i / r) / 2;
+ if (next >= r) break;
+ r = next;
+ ++c;
+ }
+ s += r;
+ }
+ *sum = s, *steps = c, *ms = now_ms() - a;
+}
+
+/* ---- bfs: the suite's graph, a search; the frontier's edges a lane each: to gathered, seen's bits asked, candidates compacted; the tail in C */
+static uint64_t gstate = 121;
+static uint64_t gnext(void) { gstate ^= gstate << 13, gstate ^= gstate >> 7, gstate ^= gstate << 17; return gstate; }
+static int cmp_u64(const void *a, const void *b) { uint64_t x = *(const uint64_t *)a, y = *(const uint64_t *)b; return (x > y) - (x < y); }
+typedef struct { uint32_t nodes; uint32_t *start, *to; } Graph;
+static Graph graph_build(uint32_t nodes, size_t edges) {
+ uint64_t *keys = malloc(edges * 8);
+ size_t c = 0;
+ gstate = 121;
+ while (c < edges) { uint64_t a = gnext() % nodes, b = gnext() % nodes; if (a != b) keys[c++] = a << 20 | b; }
+ qsort(keys, c, 8, cmp_u64);
+ size_t u = 0;
+ for (size_t i = 0; i < c; ++i) if (!u || keys[i] != keys[u - 1]) keys[u++] = keys[i];
+ Graph g = {nodes, calloc(nodes + 1, 4), malloc(4 * u)};
+ for (size_t i = 0; i < u; ++i) ++g.start[(keys[i] >> 20) + 1], g.to[i] = keys[i] & 0xFFFFF;
+ for (uint32_t v = 0; v < nodes; ++v) g.start[v + 1] += g.start[v];
+ free(keys);
+ return g;
+}
+enum { BATCH = 1024, LANES = 65536 };
+static uint32_t rows[LANES + 16], owner[LANES], cand[LANES];
+
+static uint32_t bfs_search(const Graph *g, uint32_t from, bool tiles, uint64_t *seen, uint32_t *depth, uint32_t *queue, double *ms) {
+ insn prog[8];
+ at = prog;
+ uint32_t *to = g->to;
+''' + emit_program('bfs', programs['bfs'], {'rows': 'rows', 'to': 'to', 'seen': 'seen', 'cand': 'cand'}) + '''
+ double a = now_ms();
+ memset(seen, 0, (g->nodes + 63) / 64 * 8);
+ size_t count = 1, i = 0;
+ queue[0] = from, seen[from >> 6] |= UINT64_C(1) << (from & 63), depth[from] = 0;
+ while (i < count) {
+ size_t end = count < i + BATCH ? count : i + BATCH, n = 0;
+ for (size_t r = i; r < end; ++r) {
+ uint32_t v = queue[r];
+ for (uint32_t e = g->start[v]; e < g->start[v + 1]; ++e) {
+ if (n == LANES) { fputs("too many edges in a batch\\n", stderr); exit(1); }
+ rows[n] = e, owner[n++] = v;
+ }
+ }
+ Ctx ctx = {.tiles = (n + T - 1) / T, .count_lanes = n, .out = cand, .count = 0};
+ if (tiles) start(prog, &ctx);
+ else {
+ static uint32_t to_[LANES], has[LANES];
+ for (size_t l = 0; l < n; ++l) to_[l] = g->to[rows[l]];
+ for (size_t l = 0; l < n; ++l) has[l] = seen[to_[l] >> 6] >> (to_[l] & 63) & 1;
+ for (size_t l = 0; l < n; ++l) if (!has[l]) cand[ctx.count++] = (uint32_t)l;
+ }
+ for (size_t k = 0; k < ctx.count; ++k) {
+ uint32_t l = cand[k], t = g->to[rows[l]];
+ if (seen[t >> 6] >> (t & 63) & 1) continue;
+ seen[t >> 6] |= UINT64_C(1) << (t & 63);
+ depth[t] = depth[owner[l]] + 1;
+ queue[count++] = t;
+ }
+ i = end;
+ }
+ *ms = now_ms() - a;
+ return (uint32_t)count;
+}
+
+static uint32_t bfs_c(const Graph *g, uint32_t from, uint64_t *seen, uint32_t *depth, uint32_t *queue, double *ms) {
+ double a = now_ms();
+ memset(seen, 0, (g->nodes + 63) / 64 * 8);
+ size_t count = 1;
+ queue[0] = from, seen[from >> 6] |= UINT64_C(1) << (from & 63), depth[from] = 0;
+ for (size_t i = 0; i < count; ++i) {
+ uint32_t v = queue[i];
+ for (uint32_t e = g->start[v]; e < g->start[v + 1]; ++e) {
+ uint32_t t = g->to[e];
+ if (seen[t >> 6] >> (t & 63) & 1) continue;
+ seen[t >> 6] |= UINT64_C(1) << (t & 63);
+ depth[t] = depth[v] + 1;
+ queue[count++] = t;
+ }
+ }
+ *ms = now_ms() - a;
+ return (uint32_t)count;
+}
+
+int main(int argc, char **argv) {
+ const char *body = argc > 1 ? argv[1] : "all";
+ int rounds = argc > 2 ? atoi(argv[2]) : 1;
+ init_vectors();
+ bool all = !strcmp(body, "all");
+ if (all || !strcmp(body, "knapsack")) {
+ unsigned n = 100, capacity = 50000;
+ uint32_t *weight = malloc(4 * n), *worth = malloc(4 * n), *best = calloc(capacity + 1, 4);
+ uint64_t x = 1;
+ for (unsigned i = 0; i < n; ++i) x = minstd(x), weight[i] = (uint32_t)(x % 1000) + 1, x = minstd(x), worth[i] = (uint32_t)(x % 1000);
+ for (int r = 0; r < rounds; ++r) {
+ double ms[3], answers[3];
+ memset(best, 0, 4 * (capacity + 1)), answers[0] = knapsack_tiles(n, capacity, weight, worth, best, &ms[0]);
+ memset(best, 0, 4 * (capacity + 1)), answers[1] = knapsack_passes(n, capacity, weight, worth, best, &ms[1]);
+ memset(best, 0, 4 * (capacity + 1)), answers[2] = knapsack_c(n, capacity, weight, worth, best, &ms[2]);
+ printf("{\\"body\\":\\"knapsack\\",\\"K\\":%d,\\"S\\":%d,\\"T\\":%d,\\"tiles_ms\\":%.3f,\\"passes_ms\\":%.3f,\\"c_ms\\":%.3f,\\"same\\":%s,\\"answer\\":%.17g}\\n", K, S, T, ms[0], ms[1], ms[2],
+ answers[0] == answers[2] && answers[1] == answers[2] ? "true" : "false", answers[0]);
+ }
+ }
+ if (all || !strcmp(body, "roots")) {
+ for (int r = 0; r < rounds; ++r) {
+ double sum[2], steps[2], ms[2];
+ roots_tiles(100000, &sum[0], &steps[0], &ms[0]);
+ roots_c(100000, &sum[1], &steps[1], &ms[1]);
+ printf("{\\"body\\":\\"roots\\",\\"K\\":%d,\\"S\\":%d,\\"T\\":%d,\\"tiles_ms\\":%.3f,\\"c_ms\\":%.3f,\\"same\\":%s,\\"sum\\":%.17g,\\"steps\\":%.0f}\\n", K, S, T, ms[0], ms[1],
+ steps[0] == steps[1] && fabs(sum[0] - sum[1]) < 1e-6 ? "true" : "false", sum[0], steps[0]);
+ }
+ }
+ if (all || !strcmp(body, "bfs")) {
+ Graph g = graph_build(100000, 1000000);
+ uint64_t *seen = calloc((g.nodes + 63) / 64, 8);
+ uint32_t *depth = calloc(g.nodes, 4), *queue = malloc(4 * g.nodes);
+ for (int r = 0; r < rounds; ++r) {
+ uint32_t from = (uint32_t)(gnext() % g.nodes), reached[3];
+ double ms[3];
+ reached[0] = bfs_search(&g, from, true, seen, depth, queue, &ms[0]);
+ reached[1] = bfs_search(&g, from, false, seen, depth, queue, &ms[1]);
+ reached[2] = bfs_c(&g, from, seen, depth, queue, &ms[2]);
+ printf("{\\"body\\":\\"bfs\\",\\"K\\":%d,\\"S\\":%d,\\"T\\":%d,\\"tiles_ms\\":%.3f,\\"passes_ms\\":%.3f,\\"c_ms\\":%.3f,\\"same\\":%s,\\"reached\\":%u}\\n", K, S, T, ms[0], ms[1], ms[2],
+ reached[0] == reached[2] && reached[1] == reached[2] ? "true" : "false", reached[0]);
+ }
+ }
+ return 0;
+}''')
generated by cgit v1.2.3 (git 2.46.0) at 2026-10-20 14:48:14 +0000