From ea9507f1c1f97e453d02f3610a03558974de5240 Mon Sep 17 00:00:00 2001 From: eunhhu Date: Thu, 8 Oct 2026 13:39:18 +0100 Subject: 042 tiles (research 133): a register-tile interpreter in two forms, computed goto in one function (make.py: GCC spills every slot across the indirect branch) and threaded tail calls with the slots as vector arguments (make2.py: an operation its vector instructions, a load and a branch; clang's musttail, GCC's sibling calls); knapsack's inner loop, BFS's frontier and roots' Newton steps against tight passes and C: roots 5.8 ms (C 6.9, the compiled code 9.57), knapsack 29.4 (passes 2.9 with clang, 9.6 gcc; C 6.2; compiled 7.9), BFS 18.1 (passes 9.8; C 8.6; compiled 7.7); NEON gathers lane by lane cost the tiles what they gain in registers --- .../experiments/042-tiles/2026-10-08-tiles.jsonl | 54 +++ research/experiments/042-tiles/make.py | 355 +++++++++++++++++++ research/experiments/042-tiles/make2.py | 385 +++++++++++++++++++++ 3 files changed, 794 insertions(+) create mode 100644 research/experiments/042-tiles/2026-10-08-tiles.jsonl create mode 100644 research/experiments/042-tiles/make.py create mode 100644 research/experiments/042-tiles/make2.py 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 \n#include \n#include \n#include \n#include \n#include \n#include ') +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 +#include +#include +#include +#include +#include +#include +#include +#include +#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; +}''') -- cgit v1.2.3