/* Research 139's deciding experiment: a batch loop's body as the program a * load would make of it (not built yet: docs/research/139), interpreted a * step at a time by one switch. Scalar steps make the uniform and affine * values as the executor's rules do (a number's arithmetic finite, an affine * value's parts whole within 2^52, a DENSE key's row whole and in range); * kernel steps run over four-byte lanes in registers fixed at load (no * buffer sought at run time); a list's cells are read once a batch; a row * range is checked at its ends; sums and products past 2^32 fail the batch; * a set waits for the batch's end, after its reads' distance is checked. * The bodies are the suite's as their bytecode has them: * matmul the j loop a batch (200 lanes), k's rounds together: * s = s + a[i * n + k].v * b[k * n + j].v * knapsack each item's rounds a batch (1024 lanes): c = capacity - k, * with = table[c - w].best + v, if with > table[c].best { set } * Each also with the same kernels called straight from C, the scalars made * in C (direct: no interpretation), against natural.h's C. * cc -O3 -falign-loops=64 -std=c17 -I../040-against-c programs.c -lm -o programs * ./programs [matmul|knapsack] [program|direct|c|all] [rounds] (one JSON line each) */ #define _POSIX_C_SOURCE 200809L #include #include #include #include #include #include #include #include "natural.h" static double now(void) { struct timespec t; clock_gettime(CLOCK_MONOTONIC, &t); return t.tv_sec * 1e3 + t.tv_nsec / 1e6; } enum { LANES = 1024, SREGS = 16, LREGS = 8 }; /* ---- Kernels (as module-batch.c's words_op, cells_at, put_words), each its own function: inlined in run()'s switch, * GCC made them scalar (their lanes' pointers from one struct) */ __attribute__((noinline)) static bool k_mulk(uint32_t *restrict out, const uint32_t *restrict xs, uint32_t k, unsigned n) { uint32_t over = 0; for (unsigned l = 0; l < n; ++l) { uint64_t p = (uint64_t)xs[l] * k; out[l] = (uint32_t)p, over |= (uint32_t)(p >> 32); } return !over; } __attribute__((noinline)) static bool k_addk(uint32_t *restrict out, const uint32_t *restrict xs, uint32_t k, unsigned n) { uint32_t over = 0; for (unsigned l = 0; l < n; ++l) { uint32_t s = xs[l] + k; out[l] = s, over |= s < k; } return !over; } __attribute__((noinline)) static bool k_add(uint32_t *restrict out, const uint32_t *restrict xs, const uint32_t *restrict ys, unsigned n) { uint32_t over = 0; for (unsigned l = 0; l < n; ++l) { uint32_t s = xs[l] + ys[l]; out[l] = s, over |= s < xs[l]; } return !over; } __attribute__((noinline)) static unsigned k_gt(uint32_t *restrict out, const uint32_t *restrict xs, const uint32_t *restrict ys, unsigned n) { uint32_t count = 0; for (unsigned l = 0; l < n; ++l) { uint32_t t = xs[l] > ys[l]; out[l] = t, count += t; } return count; } /* a one-byte column's cells at rows first + lane: four-byte lanes */ __attribute__((noinline)) static void k_widen(uint32_t *restrict out, const uint8_t *restrict from, unsigned n) { for (unsigned l = 0; l < n; ++l) out[l] = from[l]; } __attribute__((noinline)) static void k_put_masked(uint32_t *restrict to, const uint32_t *restrict from, const uint32_t *restrict mask, unsigned n) { for (unsigned l = 0; l < n; ++l) to[l] = mask[l] ? from[l] : to[l]; } /* ---- Lists as a batch reads them: a DENSE list's column (each key its row) */ typedef struct { size_t count; void *column; uint8_t width; } Cells; /* ---- The program */ enum { S_FIELD, /* s[d] = the message's field arg */ S_LOCAL, /* s[d] = local arg (a local the block does not set) */ S_ROUND, /* s[d] = the nested loop's round (uniform) */ S_LANE, /* s[d] = the batch loop's value (affine: first, lane_step) */ S_MUL, /* s[d] = s[a] * s[b] (uniform by uniform, or uniform by affine) */ S_ADD, /* s[d] = s[a] + s[b] */ S_SUB, /* s[d] = s[a] - s[b] */ S_GET, /* s[d] = list arg's field at the row of key s[a] (uniform; DENSE: the key its row) */ K_READ, /* L[d] = list arg's field at affine rows s[a] (four bytes: the column itself; one byte: made four) */ K_MULK, /* L[d] = L[a] * s[b] */ K_ADDK, /* L[d] = L[a] + s[b] */ K_ADD, /* L[d] = L[a] + L[b] */ K_GT, /* L[d] = L[a] > L[b], its 1s counted (a branch's lanes) */ K_IF, /* none of L[a]'s lanes: on past the next b steps */ K_SET, /* list arg's field at affine rows s[a] = L[b] in L[d]'s lanes, at the batch's end */ }; typedef struct { uint8_t op, d, a, b; uint32_t arg; } Step; typedef struct { double value, step; bool affine; } Scalar; typedef struct { uint32_t list; Scalar rows; const uint32_t *value, *mask; } Set; typedef struct { Scalar s[SREGS]; const uint32_t *lane[LREGS]; /* each lane register's lanes: its own, or a column's */ uint32_t own[LREGS][2][LANES + 24]; /* each two (laid out at load, not 4 KB apart as the executor's buffers): a result never its operand's */ unsigned ones[LREGS]; const uint32_t *field; /* the message */ const double *local; double round, first, lane_step; unsigned n; Cells cells[4]; /* each list's, once a batch */ Set sets[2]; unsigned set_count; bool failed; } Machine; static bool whole(double x) { return x == floor(x) && fabs(x) <= 0x1p52 && !(x == 0 && signbit(x)); } /* value + step × lane, as affine_of: whole parts at both ends (else the executor makes lanes; here the batch fails) */ static bool affine(Machine *M, double value, double step, Scalar *out) { if (!whole(value) || !whole(step) || !whole(value + step * (M->n - 1))) return false; *out = (Scalar){value, step, step != 0}; return true; } /* a op b of scalars as the executor's rules: uniforms as numbers (finite), affine values as affine_op */ static bool scalar_op(Machine *M, uint8_t op, const Scalar *a, const Scalar *b, Scalar *out) { if (!a->affine && !b->affine) { double v = op == S_MUL ? a->value * b->value : op == S_ADD ? a->value + b->value : a->value - b->value; if (!isfinite(v)) return false; *out = (Scalar){v, 0, false}; return true; } if ((!a->affine && !whole(a->value)) || (!b->affine && !whole(b->value))) return false; double as = a->affine ? a->step : 0, bs = b->affine ? b->step : 0; switch (op) { case S_ADD: return affine(M, a->value + b->value, as + bs, out); case S_SUB: return affine(M, a->value - b->value, as - bs, out); default: if (!a->affine && a->value > 0) return affine(M, a->value * b->value, a->value * bs, out); if (!b->affine && b->value > 0) return affine(M, a->value * b->value, as * b->value, out); return false; } } /* Where register d's result goes: of its two buffers the one neither operand holds (s = s + e: the other). */ static uint32_t *target(Machine *M, unsigned d, unsigned a, unsigned b) { uint32_t *out = M->own[d][0]; return out == M->lane[a] || out == M->lane[b] ? M->own[d][1] : out; } /* The program's steps once for the batch's lanes (a nested loop's round once). false: the batch fails (the executor's rounds). */ static bool run(const Step *prog, unsigned count, Machine *M) { unsigned n = M->n; for (unsigned i = 0; i < count; ++i) { const Step *p = &prog[i]; switch (p->op) { case S_FIELD: M->s[p->d] = (Scalar){M->field[p->arg], 0, false}; break; case S_LOCAL: M->s[p->d] = (Scalar){M->local[p->arg], 0, false}; break; case S_ROUND: M->s[p->d] = (Scalar){M->round, 0, false}; break; case S_LANE: if (!affine(M, M->first, M->lane_step, &M->s[p->d])) return false; break; case S_MUL: case S_ADD: case S_SUB: if (!scalar_op(M, p->op, &M->s[p->a], &M->s[p->b], &M->s[p->d])) return false; break; case S_GET: { const Cells *c = &M->cells[p->arg]; double k = M->s[p->a].value; if (M->s[p->a].affine || !(k >= 0 && k < (double)c->count && k == floor(k))) return false; size_t row = (size_t)k; M->s[p->d] = (Scalar){c->width == 1 ? ((const uint8_t *)c->column)[row] : ((const uint32_t *)c->column)[row], 0, false}; break; } case K_READ: { const Cells *c = &M->cells[p->arg]; const Scalar *r = &M->s[p->a]; double end = r->value + r->step * (n - 1); if (!r->affine || r->step != 1 || !(r->value >= 0 && end < (double)c->count)) return false; /* (the executor's other ways) */ size_t first = (size_t)r->value; if (c->width == 4) M->lane[p->d] = (const uint32_t *)c->column + first; /* a view */ else k_widen(M->own[p->d][0], (const uint8_t *)c->column + first, n), M->lane[p->d] = M->own[p->d][0]; break; } case K_MULK: { uint32_t *out = target(M, p->d, p->a, p->a); if (M->s[p->b].affine || !(M->s[p->b].value >= 0 && M->s[p->b].value < 4294967296.0) || !k_mulk(out, M->lane[p->a], (uint32_t)M->s[p->b].value, n)) return false; M->lane[p->d] = out; break; } case K_ADDK: { uint32_t *out = target(M, p->d, p->a, p->a); if (M->s[p->b].affine || !(M->s[p->b].value >= 0 && M->s[p->b].value < 4294967296.0) || !k_addk(out, M->lane[p->a], (uint32_t)M->s[p->b].value, n)) return false; M->lane[p->d] = out; break; } case K_ADD: { uint32_t *out = target(M, p->d, p->a, p->b); if (!k_add(out, M->lane[p->a], M->lane[p->b], n)) return false; M->lane[p->d] = out; break; } case K_GT: { uint32_t *out = target(M, p->d, p->a, p->b); M->ones[p->d] = k_gt(out, M->lane[p->a], M->lane[p->b], n); M->lane[p->d] = out; break; } case K_IF: if (!M->ones[p->a]) i += p->b; break; case K_SET: if (M->set_count == 2) return false; M->sets[M->set_count++] = (Set){p->arg, M->s[p->a], M->lane[p->b], M->ones[p->d] == n ? NULL : M->lane[p->d]}; break; } } return true; } /* ---- matmul */ static uint64_t minstd_x; static uint32_t minstd_next(void) { return (uint32_t)(minstd_x = minstd_x * 48271 % 2147483647); } /* The j loop a batch (as the executor runs it: the lanes j, k's rounds together), s in lane register 0. */ static double matmul_program(uint32_t n, const uint8_t *a, const uint8_t *b, double *last, uint64_t *work) { static Machine M; enum { R_I, R_N, R_IN, R_K, R_KEYA, R_AV, R_KN, R_J, R_ROWS }; const Step body[] = { {S_LOCAL, R_I, 0, 0, 0}, {S_FIELD, R_N, 0, 0, 0}, {S_MUL, R_IN, R_I, R_N, 0}, {S_ROUND, R_K, 0, 0, 0}, {S_ADD, R_KEYA, R_IN, R_K, 0}, {S_GET, R_AV, R_KEYA, 0, 0}, /* a[i * n + k].v */ {S_MUL, R_KN, R_K, R_N, 0}, {S_LANE, R_J, 0, 0, 0}, {S_ADD, R_ROWS, R_KN, R_J, 0}, {K_READ, 1, R_ROWS, 0, 1}, /* b[k * n + j].v */ {K_MULK, 2, 1, R_AV, 0}, {K_ADD, 0, 0, 2, 0}, /* s = s + ... */ }; double local[1], total = 0; uint32_t field[1] = {n}; M.field = field, M.local = local, M.n = n, M.first = 0, M.lane_step = 1; for (uint32_t i = 0; i < n; ++i) { local[0] = i; M.cells[0] = (Cells){(size_t)n * n, (void *)a, 1}, M.cells[1] = (Cells){(size_t)n * n, (void *)b, 1}; /* once a batch */ memset(M.own[0][0], 0, n * sizeof(uint32_t)), M.lane[0] = M.own[0][0]; /* let s = 0 */ for (uint32_t k = 0; k < n; ++k) { M.round = k; *work += n * (uint64_t)(sizeof(body) / sizeof(body[0]) + 1); if (!run(body, sizeof(body) / sizeof(body[0]), &M)) return NAN; } uint64_t sum = 0; for (uint32_t l = 0; l < n; ++l) sum += M.lane[0][l]; total += (double)sum, *last = M.lane[0][n - 1]; } return total; } /* The same kernels called straight, the scalars made in C. */ static double matmul_direct(uint32_t n, const uint8_t *a, const uint8_t *b, double *last) { static uint32_t s[2][LANES + 24], x[LANES + 24], y[LANES + 24]; /* s two (a result never its operand), as the program's */ double total = 0; for (uint32_t i = 0; i < n; ++i) { unsigned t = 0; memset(s[0], 0, n * sizeof(uint32_t)); for (uint32_t k = 0; k < n; ++k, t ^= 1) { k_widen(x, b + (size_t)k * n, n); if (!k_mulk(y, x, a[(size_t)i * n + k], n) || !k_add(s[t ^ 1], s[t], y, n)) return NAN; } uint64_t sum = 0; for (uint32_t l = 0; l < n; ++l) sum += s[t][l]; total += (double)sum, *last = s[t][n - 1]; } return total; } /* ---- knapsack */ /* Each item's rounds a batch, lanes in memory's order (the executor numbers them so after its first batch). */ static double knapsack_program(uint32_t items, uint32_t capacity, const uint32_t *weight, const uint32_t *worth, uint32_t *best, uint64_t *work) { static Machine M; enum { R_CAP, R_K, R_C, R_W, R_CW, R_V }; const Step body[] = { {S_FIELD, R_CAP, 0, 0, 0}, {S_LANE, R_K, 0, 0, 0}, {S_SUB, R_C, R_CAP, R_K, 0}, /* c = capacity - k (lanes: -k) */ {S_LOCAL, R_W, 0, 0, 0}, {S_SUB, R_CW, R_C, R_W, 0}, /* c - it.weight */ {K_READ, 1, R_CW, 0, 0}, {S_LOCAL, R_V, 0, 0, 1}, {K_ADDK, 2, 1, R_V, 0}, /* with = table[c - w].best + it.worth */ {K_READ, 3, R_C, 0, 0}, {K_GT, 4, 2, 3, 0}, /* with > table[c].best */ {K_IF, 0, 4, 1, 0}, {K_SET, 4, R_C, 2, 0}, /* table[c].best = with */ }; double local[2]; uint32_t field[1] = {capacity}; M.field = field, M.local = local; for (uint32_t i = 0; i < items; ++i) { local[0] = weight[i], local[1] = worth[i]; uint32_t rounds = capacity + 1 - weight[i]; for (uint32_t k0 = 0; k0 < rounds; k0 += LANES) { unsigned n = rounds - k0 < LANES ? rounds - k0 : LANES; /* lanes in memory's order: lane l is round k0 + n - 1 - l (c = capacity - k ascending in l) */ M.n = n, M.first = (double)(k0 + n - 1), M.lane_step = -1, M.set_count = 0; field[0] = capacity; M.cells[0] = (Cells){(size_t)capacity + 1, best, 4}; /* once a batch */ *work += n * (uint64_t)(sizeof(body) / sizeof(body[0]) + 1); if (!run(body, sizeof(body) / sizeof(body[0]), &M)) return NAN; for (unsigned t = 0; t < M.set_count; ++t) { /* the end: the set's rows against the reads' (distance w >= 0: none of these rounds reads what an earlier one sets), then set */ const Set *s = &M.sets[t]; double d = M.s[R_CW].value - s->rows.value; /* the read c - w against the set c, both step 1 */ if (!(d <= 0 && d == floor(d))) return NAN; uint32_t *to = best + (size_t)s->rows.value; if (s->mask) k_put_masked(to, s->value, s->mask, n); else memcpy(to, s->value, n * sizeof(uint32_t)); } } } return best[capacity]; } static double knapsack_direct(uint32_t items, uint32_t capacity, const uint32_t *weight, const uint32_t *worth, uint32_t *best) { static uint32_t with[LANES], gt[LANES]; for (uint32_t i = 0; i < items; ++i) { uint32_t w = weight[i], v = worth[i], rounds = capacity + 1 - w; for (uint32_t k0 = 0; k0 < rounds; k0 += LANES) { unsigned n = rounds - k0 < LANES ? rounds - k0 : LANES; size_t low = capacity - (k0 + n - 1); if (!k_addk(with, best + low - w, v, n)) return NAN; unsigned ones = k_gt(gt, with, best + low, n); if (!ones) continue; if (ones == n) memcpy(best + low, with, n * sizeof(uint32_t)); else k_put_masked(best + low, with, gt, n); } } return best[capacity]; } /* ---- */ static int cmp(const void *a, const void *b) { double x = *(const double *)a, y = *(const double *)b; return (x > y) - (x < y); } int main(int argc, char **argv) { const char *body = argc > 1 ? argv[1] : "all", *way = argc > 2 ? argv[2] : "all"; int rounds = argc > 3 ? atoi(argv[3]) : 11; bool all = !strcmp(way, "all"); double *ms = malloc(rounds * sizeof(double)), out[4]; if (!strcmp(body, "all") || !strcmp(body, "matmul")) { uint32_t n = 200, in[2] = {n, 0}; uint8_t *a = malloc((size_t)n * n), *b = malloc((size_t)n * n); minstd_x = 1; for (uint32_t i = 0; i < n * n; ++i) a[i] = (uint8_t)(minstd_next() % 100), b[i] = (uint8_t)(minstd_next() % 100); c_matmul(in, out); double want = out[0]; for (int w = 0; w < 3; ++w) { const char *name = w == 0 ? "program" : w == 1 ? "direct" : "c"; if (!all && strcmp(way, name)) continue; double got = 0, last = 0; uint64_t work = 0; for (int r = 0; r < rounds; ++r) { double t = now(); if (w == 0) got = matmul_program(n, a, b, &last, &work); else if (w == 1) got = matmul_direct(n, a, b, &last); else c_matmul(in, out), got = out[0]; ms[r] = now() - t; } qsort(ms, rounds, sizeof(double), cmp); printf("{\"body\":\"matmul\",\"way\":\"%s\",\"ms\":%.3f,\"same\":%s}\n", name, ms[rounds / 2], got == want ? "true" : "false"); } } if (!strcmp(body, "all") || !strcmp(body, "knapsack")) { uint32_t items = 100, capacity = 50000, in[2] = {items, capacity}; uint32_t *weight = malloc(4 * items), *worth = malloc(4 * items), *best = malloc(4 * (capacity + 1)); minstd_x = 1; for (uint32_t i = 0; i < items; ++i) weight[i] = minstd_next() % 1000 + 1, worth[i] = minstd_next() % 1000; c_knapsack(in, out); double want = out[0]; for (int w = 0; w < 3; ++w) { const char *name = w == 0 ? "program" : w == 1 ? "direct" : "c"; if (!all && strcmp(way, name)) continue; double got = 0; uint64_t work = 0; for (int r = 0; r < rounds; ++r) { memset(best, 0, 4 * (capacity + 1)); double t = now(); if (w == 0) got = knapsack_program(items, capacity, weight, worth, best, &work); else if (w == 1) got = knapsack_direct(items, capacity, weight, worth, best); else c_knapsack(in, out), got = out[0]; ms[r] = now() - t; } qsort(ms, rounds, sizeof(double), cmp); printf("{\"body\":\"knapsack\",\"way\":\"%s\",\"ms\":%.3f,\"same\":%s}\n", name, ms[rounds / 2], got == want ? "true" : "false"); } } return 0; }