diff --git a/.github/goose.sh b/.github/goose.sh new file mode 100755 index 000000000..a99438251 --- /dev/null +++ b/.github/goose.sh @@ -0,0 +1,12 @@ +#!/usr/bin/bash + +set -euo pipefail + +VERSION=v2026.2 +NAME=goose-$VERSION-linux-x64 +mkdir -p /tmp/goose +cd /tmp/goose +curl -sSLf -o $NAME.tar.gz https://github.com/aardappel/goose/releases/download/$VERSION/$NAME.tar.gz +tar -xzf $NAME.tar.gz +sudo ln -sf $(pwd)/$NAME/goose /usr/bin/goose +goose version diff --git a/.github/workflows/bench.yml b/.github/workflows/bench.yml index f20378144..826f6bd76 100644 --- a/.github/workflows/bench.yml +++ b/.github/workflows/bench.yml @@ -36,6 +36,7 @@ jobs: # elixir, # fortran, go, + goose, # hacklang, # haskell, haxe, @@ -157,6 +158,9 @@ jobs: - name: Install odin if: matrix.lang == 'odin' run: ./.github/odin.sh + - name: Install goose + if: matrix.lang == 'goose' + run: ./.github/goose.sh - uses: dlang-community/setup-dlang@v1 if: matrix.lang == 'd' with: diff --git a/bench/algorithm/binarytrees/1.goose b/bench/algorithm/binarytrees/1.goose new file mode 100644 index 000000000..21daf87ba --- /dev/null +++ b/bench/algorithm/binarytrees/1.goose @@ -0,0 +1,49 @@ +// Every tree is a pool of 8-byte nodes local to the function that builds it, +// linked by 4-byte self-relative references; the tree is freed when that +// function returns, by resetting one pointer. + +import std; +import os; + +let MIN_DEPTH = 4; + +struct Node { l: Node&?, r: Node&? } + +recursive fn make(pool: Node[>..]&, depth: i64) -> Node& { + if depth == 0 { return pool.push(Node {}); } + let l: Node& = make(pool, depth - 1); + let r: Node& = make(pool, depth - 1); + return pool.push(Node { l: l, r: r }); +} + +recursive fn check(n: Node&) -> i64 { + let l = n.l; + guard l else { return 1; } + let r = n.r; + guard r else { return 1; } + return 1 + check(l) + check(r); +} + +fn tree_check(depth: i64) -> i64 { + var pool: Node[>..] = []; + return check(make(pool, depth)); +} + +fn main() { + let a = args(); + var n = 10; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + let max_depth = max(MIN_DEPTH + 2, n); + print("stretch tree of depth ", max_depth + 1, "\t check: ", tree_check(max_depth + 1)); + var long: Node[>..] = []; + let long_root: Node& = make(long, max_depth); + var d = MIN_DEPTH; + while d <= max_depth { + let iterations = 1 << (max_depth - d + MIN_DEPTH); + var chk = 0; + for i in iterations { chk += tree_check(d); } + print(iterations, "\t trees of depth ", d, "\t check: ", chk); + d += 2; + } + print("long lived tree of depth ", max_depth, "\t check: ", check(long_root)); +} diff --git a/bench/algorithm/binarytrees/2-m.goose b/bench/algorithm/binarytrees/2-m.goose new file mode 100644 index 000000000..0b698a3ec --- /dev/null +++ b/bench/algorithm/binarytrees/2-m.goose @@ -0,0 +1,89 @@ +// 1.goose's trees, built on a pool of worker threads, one per hardware +// thread. The stretch tree is one job and each depth's iterations are split +// into 16 jobs of equal work; a job's trees are built one after another in +// the worker's own pool, and only the check sums come back. The main thread +// builds the long-lived tree while the workers run. + +import std; +import os; + +let MIN_DEPTH = 4; +let SPLIT = 16; + +struct Node { l: Node&?, r: Node&? } + +// count trees of the given depth; slot -1 is the stretch tree, slot i the +// depth MIN_DEPTH + 2i, and depth -1 tells a worker to stop. +struct Job { depth: i64, count: i64, slot: i64 } +struct Done { slot: i64, check: i64 } + +recursive fn make(pool: Node[>..]&, depth: i64) -> Node& { + if depth == 0 { return pool.push(Node {}); } + let l: Node& = make(pool, depth - 1); + let r: Node& = make(pool, depth - 1); + return pool.push(Node { l: l, r: r }); +} + +recursive fn check(n: Node&) -> i64 { + let l = n.l; + guard l else { return 1; } + let r = n.r; + guard r else { return 1; } + return 1 + check(l) + check(r); +} + +fn tree_check(depth: i64) -> i64 { + var pool: Node[>..] = []; + return check(make(pool, depth)); +} + +thread_fn worker() { + loop { + let job = qget(); + guard job.depth >= 0 else { break; } + var chk = 0; + for i in job.count { chk += tree_check(job.depth); } + qput(Done { slot: job.slot, check: chk }); + } +} + +fn main() { + let a = args(); + var n = 10; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + let max_depth = max(MIN_DEPTH + 2, n); + let depths = (max_depth - MIN_DEPTH) / 2 + 1; + + var jobs = 1; + qput(Job { depth: max_depth + 1, count: 1, slot: -1 }); + for i in depths { + let iterations = 1 << (max_depth - 2 * i); + let split = min(SPLIT, iterations); + for j in split { + qput(Job { depth: MIN_DEPTH + 2 * i, count: iterations / split, slot: i }); + } + jobs += split; + } + let workers = max(1, hardware_threads()); + for i in workers { qput(Job { depth: -1, count: 0, slot: 0 }); } + var ids: i64[>..] = []; + for i in workers { ids.push(thread_spawn(worker)); } + + var long: Node[>..] = []; + let long_root: Node& = make(long, max_depth); + + var stretch = 0; + var sums: i64[>..] = []; + sums.resize(depths, 0); + for k in jobs { + let d = qget(); + if d.slot < 0 { stretch = d.check; } else { sums[d.slot] += d.check; } + } + for id in ids { thread_wait(id); } + + print("stretch tree of depth ", max_depth + 1, "\t check: ", stretch); + for i in depths { + print(1 << (max_depth - 2 * i), "\t trees of depth ", MIN_DEPTH + 2 * i, "\t check: ", sums[i]); + } + print("long lived tree of depth ", max_depth, "\t check: ", check(long_root)); +} diff --git a/bench/algorithm/coro-prime-sieve/1.goose b/bench/algorithm/coro-prime-sieve/1.goose new file mode 100644 index 000000000..69c992207 --- /dev/null +++ b/bench/algorithm/coro-prime-sieve/1.goose @@ -0,0 +1,48 @@ +// Goose has no coroutines, and its threads cannot form this chain: there is +// one queue per element type, so n filters cannot each have a channel of +// their own, and n is an input of thousands where a worker is an OS thread +// with data stacks of its own. So the chain is the one the generator +// programs (2.cs, 2.d, 2.py) build, a generator followed by one lazy filter +// per prime found, each pulling the next value from the stage before it. +// A stage's state lives in an array instead of a heap-allocated frame: +// the generator's counter, and a filter's prime. Asking stage k for a value +// is a call and handing it back a return, as resume and suspend are in the +// Zig program; every value visits the filters in the order the channels +// would hand it on, and a new filter joins at the end of the chain. The +// values are u32, as in the Zig program. + +import std; +import os; + +// The next value of stage k, k == 0 being the generator. +recursive fn pull(primes: u32[:], k: i64, gen: u32&) -> u32 { + if k == 0 { + let v = gen; + gen += 1; + return v; + } + // The prime is read inside the loop: hoisted into a local, the compiler + // precomputes a reciprocal for the `%` with a 128-bit division per call, + // where a call here divides about once. + loop { + let v = pull(primes, k - 1, gen); + if v % primes[k - 1] != 0 { return v; } + } +} + +fn main() { + let a = args(); + var n = 100; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + var primes: u32[>..] = []; + var gen: u32 = 2; + // Goose's stdout is unbuffered, so the lines are collected and written + // at once, as the references' buffered writers do. + var out: u8[>..] = []; + for i in n { + let prime = pull(primes, primes.len, gen); + format(out, prime, "\n"); + primes.push(prime); + } + write_stdout(out); +} diff --git a/bench/algorithm/edigits/1.goose b/bench/algorithm/edigits/1.goose new file mode 100644 index 000000000..9a479f234 --- /dev/null +++ b/bench/algorithm/edigits/1.goose @@ -0,0 +1,864 @@ +// Digits of e by binary splitting, as the Rust programs compute them: +// sum_terms gives p / q = 1/1! + ... + 1/(k-1)!, and the digits are +// (p + q) * 10^(n-1) / q in decimal. The integers are 64-bit limbs, least +// significant first, in slices of one preallocated arena, with the +// algorithms of the Rust programs' libraries: multiplication by Toom-3 and +// Karatsuba down to the schoolbook, division by halves of the quotient down +// to the schoolbook (as GMP's dcpi1_div_qr), and the decimal conversion by +// dividing by 10^(19*2^i) recursively. + +import std; +import os; +import math; + +// Limbs from which Karatsuba and then Toom-3 beat the method below them. +let KARATSUBA = 24; +let KARATSUBA_SQR = 40; +let TOOM3 = 96; +let TOOM3_SQR = 96; +let LN_10 = 2.302585092994046; + +// The high half of a 64x64-bit product. clang compiles the pattern to one +// mulx, but only while it sees it on its own: inlined first, it rewrites the +// halves of an operand it knows something about (a constant, a square) and +// the pattern is gone. `recursive` is what keeps a Goose function out of line. +recursive fn mulhi(a: u64, b: u64) -> u64 { + let al = a & 0xffffffff; + let ah = a >> 32; + let bl = b & 0xffffffff; + let bh = b >> 32; + let ll = al * bl; + let lh = al * bh; + let hl = ah * bl; + let hh = ah * bh; + let mid = (ll >> 32) + (lh & 0xffffffff) + (hl & 0xffffffff); + return hh + (lh >> 32) + (hl >> 32) + (mid >> 32); +} + +// --- limb vectors ------------------------------------------------------------ + +fn top(a: u64[:]) -> i64 { + var n = a.len; + while n > 0 && a[n - 1] == 0 { n--; } + return n; +} + +fn norm(a: u64[:]) -> u64[:] { a[..top(a)] } + +fn zero(r: u64[:]) { for &x in r { x = 0; } } + +// r = a + b over r.len limbs; the carry out +fn add_n(r: u64[:], a: u64[:], b: u64[:]) -> u64 { + let n = r.len; + assert(a.len == n && b.len == n); + var c: u64 = 0; + for i in n { + let x = a[i] + b[i]; + let y = x + c; + c = if x < a[i] || y < x { 1 } else { 0 }; + r[i] = y; + } + return c; +} + +// r = a - b over r.len limbs; the borrow out +fn sub_n(r: u64[:], a: u64[:], b: u64[:]) -> u64 { + let n = r.len; + assert(a.len == n && b.len == n); + var c: u64 = 0; + for i in n { + let x = a[i] - b[i]; + let y = x - c; + c = if a[i] < b[i] || x < c { 1 } else { 0 }; + r[i] = y; + } + return c; +} + +fn add_1(r: u64[:], c: u64) -> u64 { + var k = c; + for &x in r { + if k == 0 { return 0; } + x += k; + k = if x < k { 1 } else { 0 }; + } + return k; +} + +fn sub_1(r: u64[:], c: u64) -> u64 { + var k = c; + for &x in r { + if k == 0 { return 0; } + let y = x - k; + k = if x < k { 1 } else { 0 }; + x = y; + } + return k; +} + +// r = a + b for a.len >= b.len and r.len == a.len; the carry out +fn add(r: u64[:], a: u64[:], b: u64[:]) -> u64 { + let n = b.len; + assert(r.len == a.len && a.len >= n); + let c = add_n(r[..n], a[..n], b); + if r.len == n { return c; } + copy_into(r[n..], a[n..]); + return add_1(r[n..], c); +} + +// r = a - b for a.len >= b.len and r.len == a.len; the borrow out +fn sub(r: u64[:], a: u64[:], b: u64[:]) -> u64 { + let n = b.len; + assert(r.len == a.len && a.len >= n); + let c = sub_n(r[..n], a[..n], b); + if r.len == n { return c; } + copy_into(r[n..], a[n..]); + return sub_1(r[n..], c); +} + +// r += b for r.len >= b.len; the carry out +fn add_in(r: u64[:], b: u64[:]) -> u64 { + let n = b.len; + let c = add_n(r[..n], r[..n], b); + return add_1(r[n..], c); +} + +// r -= b for r.len >= b.len; the borrow out +fn sub_in(r: u64[:], b: u64[:]) -> u64 { + let n = b.len; + let c = sub_n(r[..n], r[..n], b); + return sub_1(r[n..], c); +} + +fn cmp(a: u64[:], b: u64[:]) -> i64 { + var n = max(a.len, b.len); + while n > 0 { + n--; + let x: u64 = if n < a.len { a[n] } else { 0 }; + let y: u64 = if n < b.len { b[n] } else { 0 }; + if x != y { return if x < y { -1 } else { 1 }; } + } + return 0; +} + +// d = |x - y| for x.len == d.len >= y.len; whether x < y +fn absdiff(d: u64[:], x: u64[:], y: u64[:]) -> bool { + if cmp(x, y) >= 0 { + sub(d, x, y); + return false; + } + let n = y.len; + sub_n(d[..n], y, x[..n]); + zero(d[n..]); + return true; +} + +// r = a * m; the carry out +fn mul_1(r: u64[:], a: u64[:], m: u64) -> u64 { + assert(r.len == a.len); + var c: u64 = 0; + for i in a.len { + var lo = a[i] * m; + var hi = mulhi(a[i], m); + lo += c; + if lo < c { hi++; } + r[i] = lo; + c = hi; + } + return c; +} + +// r += a * m; the carry out +fn addmul_1(r: u64[:], a: u64[:], m: u64) -> u64 { + assert(r.len == a.len); + var c: u64 = 0; + for i in a.len { + var lo = a[i] * m; + var hi = mulhi(a[i], m); + lo += c; + if lo < c { hi++; } + let x = r[i] + lo; + if x < lo { hi++; } + r[i] = x; + c = hi; + } + return c; +} + +fn mul_basecase(r: u64[:], a: u64[:], b: u64[:]) { + let n = a.len; + assert(r.len == n + b.len && b.len >= 1); + r[n] = mul_1(r[..n], a, b[0]); + for j in 1..b.len { r[j + n] = addmul_1(r[j..j + n], a, b[j]); } +} + +// r = a^2: the products of distinct limbs once, doubled, and the squares of +// the limbs added +fn sqr_basecase(r: u64[:], a: u64[:]) { + let n = a.len; + assert(r.len == 2 * n && n >= 1); + r[0] = 0; + r[2 * n - 1] = 0; + if n > 1 { + r[n] = mul_1(r[1..n], a[1..], a[0]); + for i in 1..n - 1 { r[n + i] = addmul_1(r[2 * i + 1..n + i], a[i + 1..], a[i]); } + var c: u64 = 0; + for &x in r[1..2 * n] { + let y = (x << 1) | c; + c = x >> 63; + x = y; + } + } + var c: u64 = 0; + for i in n { + let lo = a[i] * a[i]; + let hi = mulhi(a[i], a[i]); + let x = r[2 * i] + lo; + let y = x + c; + let k: u64 = if x < lo || y < x { 1 } else { 0 }; + r[2 * i] = y; + let z = r[2 * i + 1] + hi; + let w = z + k; + c = if z < hi || w < z { 1 } else { 0 }; + r[2 * i + 1] = w; + } +} + +// x / 3 for x a multiple of 3, two's complement over x.len limbs +fn divexact3(x: u64[:]) { + var c: u64 = 0; + for &w in x { + let b: u64 = if w < c { 1 } else { 0 }; + let q = (w - c) * 0xAAAAAAAAAAAAAAAB; // the inverse of 3 modulo 2^64 + w = q; + c = b; + if q >= 0x5555555555555556 { c++; } + if q >= 0xAAAAAAAAAAAAAAAB { c++; } + } +} + +// x / 2 for x even, two's complement over x.len limbs +fn halve(x: u64[:]) { + let n = x.len; + for i in n - 1 { x[i] = (x[i] >> 1) | (x[i + 1] << 63); } + x[n - 1] = ((x[n - 1] as! i64) >> 1) as! u64; +} + +fn negate(x: u64[:]) { + for &w in x { w = ~w; } + add_1(x, 1); +} + +// x(1) = x0 + x1 + x2, |x(-1)| = |x0 - x1 + x2| and |x(-2)| = |x0 - 2 x1 + 4 x2| +// into k + 1 limbs each, for x0 and x1 of k limbs and x2 of at most k; the +// signs of x(-1) and x(-2), and t is scratch of k + 1 limbs. +fn toom3_eval(e1: u64[:], em1: u64[:], em2: u64[:], x0: u64[:], x1: u64[:], x2: u64[:], t: u64[:]) -> bool, bool { + let k = x0.len; + em1[k] = add(em1[..k], x0, x2); + e1[k] = em1[k] + add(e1[..k], em1[..k], x1); + let neg1 = cmp(em1, x1) < 0; + if neg1 { sub_n(em1[..k], x1, em1[..k]); } else { sub_in(em1, x1); } + zero(em2); + em2[x2.len] = shl_bits(em2[..x2.len], x2, 2); + em2[k] += add(em2[..k], em2[..k], x0); + t[k] = shl_bits(t[..k], x1, 1); + let neg2 = cmp(em2, t) < 0; + if neg2 { sub_n(em2, t, em2); } else { sub_n(em2, em2, t); } + return neg1, neg2; +} + +// The Toom-3 recombination of w0 = r[..2k], winf = r[4k..] and the values +// at 1, -1 and -2, the last two signed, each of 2k + 2 limbs, into r. +fn toom3_interpolate(r: u64[:], k: i64, w1: u64[:], wm1: u64[:], wm2: u64[:]) { + let w0 = r[..2 * k]; + let winf = r[4 * k..]; + // r3 = (w(-2) - w(1)) / 3, r1 = (w(1) - w(-1)) / 2, r2 = w(-1) - w(0) + sub_n(wm2, wm2, w1); + divexact3(wm2); + sub_n(w1, w1, wm1); + halve(w1); + sub_in(wm1, w0); + // r3 = (r2 - r3) / 2 + 2 w(inf), r2 = r2 + r1 - w(inf), r1 = r1 - r3 + sub_n(wm2, wm1, wm2); + halve(wm2); + add_in(wm2, winf); + add_in(wm2, winf); + add_n(wm1, wm1, w1); + sub_in(wm1, winf); + sub_n(w1, w1, wm2); + zero(r[2 * k..4 * k]); + for c, i in [w1, wm1, wm2] { + let rd = r[(i + 1) * k..]; + add_in(rd, c[..min(c.len, rd.len)]); // beyond r's end the limbs are zero + } +} + +// --- multiplication ----------------------------------------------------------- + +// r = a^2 by mul's methods with the operand shared, Karatsuba's step being +// a^2 = a1^2 B^2h + (a0^2 + a1^2 - d^2) B^h + a0^2 for d = |a0 - a1|. +recursive fn sqr(r: u64[:], a: u64[:], s: u64[:]) { + let n = a.len; + assert(r.len == 2 * n); + if n < KARATSUBA_SQR { sqr_basecase(r, a); return; } + if n >= TOOM3_SQR { + let k = (n + 2) / 3; + let l = 2 * k + 2; + let e1 = s[..k + 1]; + let em1 = s[k + 1..2 * k + 2]; + let em2 = s[2 * k + 2..3 * k + 3]; + let w1 = s[3 * k + 3..3 * k + 3 + l]; + let wm1 = s[3 * k + 3 + l..3 * k + 3 + 2 * l]; + let wm2 = s[3 * k + 3 + 2 * l..3 * k + 3 + 3 * l]; + let rest = s[3 * k + 3 + 3 * l..]; + toom3_eval(e1, em1, em2, a[..k], a[k..2 * k], a[2 * k..], rest[..k + 1]); + sqr(w1, e1, rest); + sqr(wm1, em1, rest); + sqr(wm2, em2, rest); + sqr(r[..2 * k], a[..k], rest); + sqr(r[4 * k..], a[2 * k..], rest); + toom3_interpolate(r, k, w1, wm1, wm2); + return; + } + let h = (n + 1) >> 1; + let a0 = a[..h]; + let a1 = a[h..]; + sqr(r[..2 * h], a0, s); + sqr(r[2 * h..], a1, s); + let d = s[..h]; + let v = s[h..3 * h]; + let t = s[3 * h..5 * h + 1]; + absdiff(d, a0, a1); + sqr(v, d, s[5 * h + 1..]); + t[2 * h] = add(t[..2 * h], r[..2 * h], r[2 * h..]); + sub_in(t, v); + let rh = r[h..]; + add_in(rh, t[..min(t.len, rh.len)]); +} + +// r = a * b, all a.len + b.len limbs of it, with s as scratch. +recursive fn mul(r: u64[:], a: u64[:], b: u64[:], s: u64[:]) { + if a.len < b.len { mul(r, b, a, s); return; } + let n = a.len; + let m = b.len; + assert(r.len == n + m && m >= 1); + if m < KARATSUBA { mul_basecase(r, a, b); return; } + if m >= TOOM3 && m > 2 * ((n + 2) / 3) { toom3(r, a, b, s); return; } + let h = (n + 1) >> 1; + if m > h { karatsuba(r, a, b, s); return; } + // b is at most half of a: a's pieces of b's length one at a time + mul(r[..2 * m], a[..m], b, s); + var off = m; + while off < n { + let len = min(m, n - off); + let t = s[..len + m]; + mul(t, a[off..off + len], b, s[len + m..]); + let c = add_n(r[off..off + m], r[off..off + m], t[..m]); + copy_into(r[off + m..off + m + len], t[m..]); + add_1(r[off + m..off + m + len], c); + off += len; + } +} + +// a = a1 B^h + a0 and b = b1 B^h + b0 with h = ceil(a.len / 2) < b.len: +// a b = a1 b1 B^2h + (a0 b0 + a1 b1 - (a0 - a1)(b0 - b1)) B^h + a0 b0 +fn karatsuba(r: u64[:], a: u64[:], b: u64[:], s: u64[:]) { + let h = (a.len + 1) >> 1; + let a0 = a[..h]; + let a1 = a[h..]; + let b0 = b[..h]; + let b1 = b[h..]; + mul(r[..2 * h], a0, b0, s); + mul(r[2 * h..], a1, b1, s); + let da = s[..h]; + let db = s[h..2 * h]; + let v = s[2 * h..4 * h]; + let t = s[4 * h..6 * h + 1]; + let na = absdiff(da, a0, a1); + let nb = absdiff(db, b0, b1); + mul(v, da, db, s[6 * h + 1..]); + t[2 * h] = add(t[..2 * h], r[..2 * h], r[2 * h..]); + if na == nb { sub_in(t, v); } else { add_in(t, v); } + let rh = r[h..]; + add_in(rh, t[..min(t.len, rh.len)]); // what lies beyond r is zero +} + +// a = a2 B^2k + a1 B^k + a0, likewise b, with k = ceil(a.len / 3) and +// 2k < b.len: the five products at 0, 1, -1, -2 and infinity, and the +// interpolation of Bodrato's sequence. +fn toom3(r: u64[:], a: u64[:], b: u64[:], s: u64[:]) { + let n = a.len; + let k = (n + 2) / 3; + let l = 2 * k + 2; + let e1 = s[..k + 1]; + let em1 = s[k + 1..2 * k + 2]; + let em2 = s[2 * k + 2..3 * k + 3]; + let f1 = s[3 * k + 3..4 * k + 4]; + let fm1 = s[4 * k + 4..5 * k + 5]; + let fm2 = s[5 * k + 5..6 * k + 6]; + let w1 = s[6 * k + 6..6 * k + 6 + l]; + let wm1 = s[6 * k + 6 + l..6 * k + 6 + 2 * l]; + let wm2 = s[6 * k + 6 + 2 * l..6 * k + 6 + 3 * l]; + let rest = s[6 * k + 6 + 3 * l..]; + let na1, na2 = toom3_eval(e1, em1, em2, a[..k], a[k..2 * k], a[2 * k..], rest[..k + 1]); + let nb1, nb2 = toom3_eval(f1, fm1, fm2, b[..k], b[k..2 * k], b[2 * k..], rest[..k + 1]); + mul(w1, e1, f1, rest); + mul(wm1, em1, fm1, rest); + if na1 != nb1 { negate(wm1); } + mul(wm2, em2, fm2, rest); + if na2 != nb2 { negate(wm2); } + mul(r[..2 * k], a[..k], b[..k], rest); + mul(r[4 * k..], a[2 * k..], b[2 * k..], rest); + toom3_interpolate(r, k, w1, wm1, wm2); +} + +// --- division ----------------------------------------------------------------- + +let DC_DIV = 48; // limbs from which dividing in halves beats the schoolbook + +fn shl_bits(r: u64[:], a: u64[:], k: i64) -> u64 { + assert(r.len == a.len); + if k == 0 { copy_into(r, a); return 0; } + var c: u64 = 0; + for i in a.len { + let x = a[i]; + r[i] = (x << k) | c; + c = x >> (64 - k); + } + return c; +} + +// r -= a * m; the borrow out +fn submul_1(r: u64[:], a: u64[:], m: u64) -> u64 { + assert(r.len == a.len); + var c: u64 = 0; + for i in a.len { + var lo = a[i] * m; + var hi = mulhi(a[i], m); + lo += c; + if lo < c { hi++; } + let x = r[i]; + if x < lo { hi++; } + r[i] = x - lo; + c = hi; + } + return c; +} + +// (u1 B + u0) / v for u1 < v and v's top bit set, by 32-bit digits +// (Hacker's Delight, divlu). +fn divlu(u1: u64, u0: u64, v: u64) -> u64 { + let b: u64 = 0x100000000; + let vn1 = v >> 32; + let vn0 = v & 0xffffffff; + let un1 = u0 >> 32; + let un0 = u0 & 0xffffffff; + var q1 = u1 / vn1; + var rhat = u1 - q1 * vn1; + while q1 >= b || q1 * vn0 > (rhat << 32) + un1 { + q1--; + rhat += vn1; + if rhat >= b { break; } + } + let un21 = (u1 << 32) + un1 - q1 * v; + var q0 = un21 / vn1; + rhat = un21 - q0 * vn1; + while q0 >= b || q0 * vn0 > (rhat << 32) + un0 { + q0--; + rhat += vn1; + if rhat >= b { break; } + } + return (q1 << 32) + q0; +} + +// floor((B^3 - 1) / (d1 B + d0)) - B, the reciprocal the 3-by-2 limb +// division multiplies by (Moller and Granlund; GMP's invert_pi1). +fn invert_pi1(d1: u64, d0: u64) -> u64 { + var v = divlu(~d1, 0xffffffffffffffff, d1); + var p = d1 * v; + p += d0; + if p < d0 { + v--; + if p >= d1 { v--; p -= d1; } + p -= d1; + } + let t0 = d0 * v; + let t1 = mulhi(d0, v); + p += t1; + if p < t1 { + v--; + if p >= d1 && (p > d1 || t0 >= d0) { v--; } + } + return v; +} + +// (n2 B^2 + n1 B + n0) / (d1 B + d0) for (n2, n1) < (d1, d0), with v from +// invert_pi1: the quotient and the two limbs of the remainder (GMP's +// udiv_qr_3by2). +fn div_3by2(n2: u64, n1: u64, n0: u64, d1: u64, d0: u64, v: u64) -> u64, u64, u64 { + var q = mulhi(n2, v); + var q0 = n2 * v; + let s0 = q0 + n1; + q += n2; + if s0 < q0 { q++; } + q0 = s0; + var r1 = n1 - d1 * q; + var r0 = n0 - d0; + r1 -= d1; + if n0 < d0 { r1--; } + let t0 = d0 * q; + let t1 = mulhi(d0, q); + if r0 < t0 { r1--; } + r0 -= t0; + r1 -= t1; + q++; + if r1 >= q0 { + q--; + let s = r0 + d0; + r1 += d1; + if s < r0 { r1++; } + r0 = s; + } + if r1 >= d1 && (r1 > d1 || r0 >= d0) { + q++; + if r0 < d0 { r1--; } + r0 -= d0; + r1 -= d1; + } + return q, r1, r0; +} + +fn cmp_n(a: u64[:], b: u64[:]) -> i64 { + assert(a.len == b.len); + var i = a.len; + while i > 0 { + i--; + if a[i] != b[i] { return if a[i] < b[i] { -1 } else { 1 }; } + } + return 0; +} + +// The schoolbook division of np by d, d.len >= 2 and normalized: the +// quotient into q (np.len - d.len limbs) and the remainder into np[..d.len]; +// returns the quotient's top limb, 0 or 1 (GMP's sbpi1_div_qr). +fn sb_div_qr(q: u64[:], np: u64[:], d: u64[:], v: u64) -> u64 { + let nn = np.len; + let dn = d.len; + let dm = dn - 2; + assert(q.len == nn - dn && dn >= 2); + var qh: u64 = 0; + let head = np[nn - dn..]; + if cmp_n(head, d) >= 0 { + sub_n(head, head, d); + qh = 1; + } + let d1 = d[dn - 1]; + let d0 = d[dn - 2]; + var idx = nn - 2; + var n1 = np[idx + 1]; + var i = nn - dn; + while i > 0 { + idx--; + var qi: u64 = 0; + if n1 == d1 && np[idx + 1] == d0 { + qi = 0xffffffffffffffff; + submul_1(np[idx - dm..idx + 2], d, qi); + n1 = np[idx + 1]; + } else { + let qq, r1, r0 = div_3by2(n1, np[idx + 1], np[idx], d1, d0, v); + qi = qq; + var n0 = r0; + n1 = r1; + let cy = submul_1(np[idx - dm..idx], d[..dm], qi); + let cy1: u64 = if n0 < cy { 1 } else { 0 }; + n0 -= cy; + let cy2 = n1 < cy1; + n1 -= cy1; + np[idx] = n0; + if cy2 { + n1 += d1 + add_n(np[idx - dm..idx + 1], np[idx - dm..idx + 1], d[..dm + 1]); + qi--; + } + } + i--; + q[i] = qi; + } + np[idx + 1] = n1; + return qh; +} + +// np (2n limbs) by d (n limbs, normalized) in two halves of the quotient, +// each the division of a top part by d's top half and a product to correct +// by: the quotient into q (n limbs), the remainder into np[..n]; returns the +// quotient's top limb (GMP's dcpi1_div_qr_n). +recursive fn dc_div_qr(q: u64[:], np: u64[:], d: u64[:], v: u64, tp: u64[:]) -> u64 { + let n = d.len; + if n < DC_DIV { return sb_div_qr(q, np, d, v); } + let lo = n >> 1; + let hi = n - lo; + var qh = dc_div_qr(q[lo..], np[2 * lo..], d[lo..], v, tp); + mul(tp[..n], q[lo..], d[..lo], tp[n..]); + var cy = sub_n(np[lo..lo + n], np[lo..lo + n], tp[..n]); + if qh != 0 { cy += sub_n(np[n..n + lo], np[n..n + lo], d[..lo]); } + while cy != 0 { + qh -= sub_1(q[lo..], 1); + cy -= add_n(np[lo..lo + n], np[lo..lo + n], d); + } + let ql = dc_div_qr(q[..lo], np[hi..hi + 2 * lo], d[hi..], v, tp); + mul(tp[..n], d[..hi], q[..lo], tp[n..]); + cy = sub_n(np[..n], np[..n], tp[..n]); + if ql != 0 { cy += sub_n(np[lo..lo + hi], np[lo..lo + hi], d[..hi]); } + while cy != 0 { + sub_1(q[..lo], 1); + cy -= add_n(np[..n], np[..n], d); + } + return qh; +} + +// A divisor shifted until its top bit is set, and the reciprocal of its +// top two limbs. +struct Divisor { shift: i64, d: u64[:], v: u64 } + +fn divisor(dst: u64[:], src: u64[:]) -> Divisor { + let m = src.len; + assert(dst.len == m); + let shift = clz(src[m - 1]); + shl_bits(dst, src, shift); + let v = if m >= 2 { invert_pi1(dst[m - 1], dst[m - 2]) } else { 0 }; + return Divisor { shift: shift, d: dst, v: v }; +} + +// q = num / div and rem = num % div, for num below div * B^div.len; q gets +// num.len - div.len + 1 limbs, rem div.len. +fn divrem(q: u64[:], rem: u64[:], num: u64[:], dv: Divisor, s: u64[:]) { + let m = dv.d.len; + assert(num.len <= 2 * m && rem.len == m); + let a = s[..2 * m]; + zero(a); + let c = shl_bits(a[..num.len], num, dv.shift); + if num.len < 2 * m { a[num.len] = c; } + let qq = s[2 * m..3 * m + 1]; + qq[m] = dc_div_qr(qq[..m], a, dv.d, dv.v, s[3 * m + 1..]); + let sh = dv.shift; + for i in m { + var w = a[i] >> sh; + if sh > 0 && i + 1 < m { w |= a[i + 1] << (64 - sh); } + rem[i] = w; + } + zero(q); + let qn = min(q.len, qq.len); + copy_into(q[..qn], qq[..qn]); +} + +// --- decimal conversion ---------------------------------------------------------- + +let CHUNK = 19; // 10^19 < 2^64 +let TEN19: u64 = 10000000000000000000; + +// The low digits of a (of at most four limbs) right-aligned in out, with +// leading zeros, by division by 10^9 on 32-bit halves. +fn base_digits(out: u8[:], a: u64[:]) { + var w: u32[8] = [0; 8]; + var n = 0; + for i in a.len { + w[2 * i] = a[i] as! u32; + w[2 * i + 1] = (a[i] >> 32) as! u32; + n = 2 * i + 2; + } + var pos = out.len; + while pos > 0 { + var rem: u64 = 0; + var j = n; + while j > 0 { + j--; + let cur = (rem << 32) | w[j]; + w[j] = (cur / 1000000000) as! u32; + rem = cur % 1000000000; + } + while n > 0 && w[n - 1] == 0 { n--; } + for d in 9 { + if pos == 0 { break; } + pos--; + out[pos] = ('0' + rem % 10) as! u8; + rem /= 10; + } + } +} + +// The digits of a, of at most CHUNK * 2^(i+1) of them, right-aligned in out, +// with leading zeros; pows[i] is 10^(CHUNK * 2^i). +recursive fn to_digits(out: u8[:], a: u64[:], pows: Divisor[:], i: i64, s: u64[:]) { + let an = norm(a); + if an.len <= 2 || i < 0 { + assert(an.len <= 4); + base_digits(out, an); + return; + } + let dv = pows[i]; + let m = dv.d.len; + if an.len < m { + to_digits(out, an, pows, i - 1, s); + return; + } + let ql = an.len - m + 1; + let q = s[..ql]; + let r = s[ql..ql + m]; + divrem(q, r, an, dv, s[ql + m..]); + let low = CHUNK << i; + to_digits(out[..out.len - low], q, pows, i - 1, s[ql + m..]); + to_digits(out[out.len - low..], r, pows, i - 1, s[ql + m..]); +} + +// --- the series ------------------------------------------------------------------ + +fn bitlen(v: i64) -> i64 { 64 - clz(v as u64) } + +// Limbs enough for p and q of sum_terms(a, b) and its intermediate results. +fn space(a: i64, b: i64) -> i64 { 2 * (((b - a) * bitlen(b) + 63) / 64 + 2) + 8 } + +// a! (1/(a+1)! + ... + 1/b!) as p / q, with q = (a+1) ... b: p and then q +// into the front of out, returning their lengths. +recursive fn sum_terms(a: i64, b: i64, out: u64[:], s: u64[:]) -> i64, i64 { + if b == a + 1 { + out[0] = 1; + out[1] = b as u64; + return 1, 1; + } + let mid = (a + b) / 2; + let ls = space(a, mid); + let rs = space(mid, b); + let left = s[..ls]; + let right = s[ls..ls + rs]; + let rest = s[ls + rs..]; + let pln, qln = sum_terms(a, mid, left, rest); + let prn, qrn = sum_terms(mid, b, right, rest); + let pl = left[..pln]; + let ql = left[pln..pln + qln]; + let pr = right[..prn]; + let qr = right[prn..prn + qrn]; + // p = pl qr + pr, q = ql qr + var pn = pln + qrn; + mul(out[..pn], pl, qr, rest); + let c = add_in(out[..pn], pr); + if c != 0 { out[pn] = c; pn++; } + pn = top(out[..pn]); + let qn = qln + qrn; + mul(out[pn..pn + qn], ql, qr, rest); + return pn, top(out[pn..pn + qn]); +} + +// The smallest k with log10(k!) >= n + 50, by Stirling's approximation. +fn terms_needed(n: i64) -> i64 { + fn ok(k: i64) -> bool { + if k <= 0 { return false; } + let kf = k as f64; + let ln_k_factorial = kf * log(kf) - kf + 0.5 * log(TAU * kf); + return ln_k_factorial / LN_10 >= n + 50; + } + var a = 0; + var b = 1; + while !ok(b) { a = b; b *= 2; } + while b - a > 1 { + let m = a + (b - a) / 2; + if ok(m) { b = m; } else { a = m; } + } + return b; +} + +fn main() { + let a = args(); + var n = 27; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + let k = terms_needed(n); + + let size = space(0, k - 1); // limbs of p and q + var mem: u64[>..] = []; + mem.resize(12 * size + 4096, 0); + let pq = mem[..size]; + let s = mem[size..]; + + let pn, qn = sum_terms(0, k - 1, pq, s); + // p += q + let qv = pq[pn..pn + qn]; + let pv = s[..max(pn, qn) + 1]; + zero(pv); + copy_into(pv[..pn], pq[..pn]); + add_in(pv, qv); + let p = norm(pv); + let q = norm(qv); + let work = s[pv.len..]; + + // 10^(n-1), by squaring from the top bit of the exponent + let e = n - 1; + let pw_cap = e * 3322 / 1000 / 64 + 4; // log2(10) < 3.322 + let pw = work[..pw_cap]; + let tmp = work[pw_cap..3 * pw_cap]; + let work2 = work[3 * pw_cap..]; + zero(pw); + pw[0] = 1; + var pl = 1; + var bit = bitlen(e) - 1; + while bit >= 0 { + sqr(tmp[..2 * pl], pw[..pl], work2); + pl = top(tmp[..2 * pl]); + copy_into(pw[..pl], tmp[..pl]); + if ((e >> bit) & 1) != 0 { + let c = mul_1(pw[..pl], pw[..pl], 10); + if c != 0 { pw[pl] = c; pl++; } + } + bit--; + } + // answer = p 10^(n-1) / q + let num = work2[..p.len + pl]; + mul(num, p, pw[..pl], work2[p.len + pl..]); + let numn = norm(num); + let w3 = work2[p.len + pl..]; + let qdiv = divisor(w3[..q.len], q); + let w4 = w3[q.len..]; + let ansl = numn.len - q.len + 1; + let ans = w4[..ansl]; + let remd = w4[ansl..ansl + q.len]; + let w5 = w4[ansl + q.len..]; + divrem(ans, remd, numn, qdiv, w5); + let answer = norm(ans); + + // the powers 10^(19 2^i) the conversion divides by + var pows: Divisor[>..] = []; + var plen = 1; + var cur = w5[..1]; + cur[0] = TEN19; + var free = w5[1..]; + var chunk = CHUNK; + while chunk < n { + pows.push(divisor(free[..plen], cur)); + free = free[plen..]; + chunk *= 2; + if chunk >= n { break; } + let sq = free[..2 * plen]; + sqr(sq, cur, free[2 * plen..]); + plen = top(sq); + cur = sq[..plen]; + free = free[2 * plen..]; + } + var text: u8[>..] = []; + text.resize(chunk, '0'); + to_digits(text, answer, pows, pows.len - 1, free); + var start = 0; + while start < text.len - 1 && text[start] == '0' { start++; } + let digits = text[start..]; + + var out: u8[>..] = []; + var i = 0; + while i < n { + if i + 10 <= n { + out.append(digits[i..i + 10]); + format(out, "\t:", i + 10, "\n"); + } else { + out.append(digits[i..n]); + for j in 10 - n % 10 { out.push(' '); } + format(out, "\t:", n, "\n"); + } + i += 10; + } + write_stdout(out); +} diff --git a/bench/algorithm/fannkuch-redux/1-m.goose b/bench/algorithm/fannkuch-redux/1-m.goose new file mode 100644 index 000000000..e48d6d99e --- /dev/null +++ b/bench/algorithm/fannkuch-redux/1-m.goose @@ -0,0 +1,159 @@ +// The single-threaded program's permutations and flips (see 1.goose), with +// the permutation indices cut into blocks that worker threads take from a +// queue; each block's checksum and largest flip count come back through +// another queue. + +import std; +import os; + +let ONES: u64 = 0xFFFFFFFFFFFFFFFF; + +// HI[c]: the top 4c bits, which a flip with s = 4c leaves in place. +let HI: u64[16] = [0x0000000000000000, 0xF000000000000000, 0xFF00000000000000, 0xFFF0000000000000, + 0xFFFF000000000000, 0xFFFFF00000000000, 0xFFFFFF0000000000, 0xFFFFFFF000000000, + 0xFFFFFFFF00000000, 0xFFFFFFFFF0000000, 0xFFFFFFFFFF000000, 0xFFFFFFFFFFF00000, + 0xFFFFFFFFFFFF0000, 0xFFFFFFFFFFFFF000, 0xFFFFFFFFFFFFFF00, 0xFFFFFFFFFFFFFFF0]; + +fn nibble_reverse(x: u64) -> u64 { + var y = ((x >> 4) & 0x0F0F0F0F0F0F0F0F) | ((x & 0x0F0F0F0F0F0F0F0F) << 4); + y = ((y >> 8) & 0x00FF00FF00FF00FF) | ((y & 0x00FF00FF00FF00FF) << 8); + y = ((y >> 16) & 0x0000FFFF0000FFFF) | ((y & 0x0000FFFF0000FFFF) << 16); + return (y >> 32) | (y << 32); +} + +// The flip counts of the permutations a and b. +fn flips2(a: u64, b: u64) -> i64, i64 { + var xa = a; + var ra = nibble_reverse(a); + var sa = (a & 15) << 2; + var ka = 0; + var xb = b; + var rb = nibble_reverse(b); + var sb = (b & 15) << 2; + var kb = 0; + while (sa & sb) != 60 { + ka += if sa != 60 { 1 } else { 0 }; + kb += if sb != 60 { 1 } else { 0 }; + let ha = ra >> sa; + let ya = (xa & HI[sa >> 2]) | ha; + ra = (ra & ~(ONES << sa)) | (xa << sa); + xa = ya; + sa = (ha & 15) << 2; + let hb = rb >> sb; + let yb = (xb & HI[sb >> 2]) | hb; + rb = (rb & ~(ONES << sb)) | (xb << sb); + xb = yb; + sb = (hb & 15) << 2; + } + return ka, kb; +} + +// Rotates the first k elements left by one. +fn rotate(x: u64, k: i64) -> u64 { + let low = ONES >> (64 - (k << 2)); + return (x & ~low) | ((x >> 4) & (low >> 4)) | ((x & 15) << ((k - 1) << 2)); +} + +// The next permutation in the counting order: the first two elements swap, +// and each counter that wraps rotates one more element in. +fn next_perm(perm: u64, count: i64[16]&) -> u64 { + var x = (perm & ~0xFF) | ((perm & 15) << 4) | ((perm >> 4) & 15); + var i = 1; + while i < 15 { + count[i] += 1; + if count[i] <= i { break; } + count[i] = 0; + i++; + x = rotate(x, i + 1); + } + return x; +} + +fn factorial(n: i64) -> i64 { + var f = 1; + for i in 2..n + 1 { f *= i; } + return f; +} + +// The permutation with index `start`, and its counters. +fn start_perm(n: i64, start: i64, count: i64[16]&) -> u64 { + var p: u8[16] = [0; 16]; + for i in 16 { p[i] = i as u8; } + var idx = start; + var i = n - 1; + while i > 0 { + let f = factorial(i); + let d = idx / f; + idx %= f; + count[i] = d; + let t = p; + for j in i + 1 { p[j] = t[(j + d) % (i + 1)]; } + i--; + } + var perm: u64 = 0; + for j in 16 { perm |= (p[j] as u64) << (j << 2); } + return ~perm; +} + +// The checksum and the largest flip count of the permutations with indices +// start .. start + len, start and len even. +fn run_block(n: i64, start: i64, len: i64) -> i64, i64 { + var count: i64[16] = [0; 16]; + var perm = start_perm(n, start, count); + var checksum = 0; + var max_flips = 0; + for pair in len >> 1 { + let a = perm; + perm = next_perm(perm, count); + let b = perm; + perm = next_perm(perm, count); + let ka, kb = flips2(a, b); + checksum += ka - kb; + max_flips = max(max_flips, max(ka, kb)); + } + return checksum, max_flips; +} + +struct Job { start: i64, len: i64 } // len 0: no more blocks +struct Part { checksum: i64, max_flips: i64 } + +thread_fn worker(n: i64) { + loop { + let job = qget(); + guard job.len > 0 else { break; } + let checksum, max_flips = run_block(n, job.start, job.len); + qput(Part { checksum: checksum, max_flips: max_flips }); + } +} + +fn main() { + let a = args(); + var n = 7; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + assert(n >= 2 && n <= 16); + let total = factorial(n); + let workers = max(1, hardware_threads()); + // Several blocks a worker, of an even length, so that blocks with more + // flips than others even out. + var size = (total + workers * 4 - 1) / (workers * 4); + size += size & 1; + var blocks = 0; + var start = 0; + while start < total { + qput(Job { start: start, len: min(size, total - start) }); + start += size; + blocks++; + } + for i in workers { qput(Job { start: 0, len: 0 }); } + var ids: i64[>..] = []; + for i in workers { ids.push(thread_spawn(worker, n)); } + var checksum = 0; + var max_flips = 0; + for i in blocks { + let p = qget(); + checksum += p.checksum; + max_flips = max(max_flips, p.max_flips); + } + for id in ids { thread_wait(id); } + print(checksum, "\nPfannkuchen(", n, ") = ", max_flips); +} diff --git a/bench/algorithm/fannkuch-redux/1.goose b/bench/algorithm/fannkuch-redux/1.goose new file mode 100644 index 000000000..df4bef080 --- /dev/null +++ b/bench/algorithm/fannkuch-redux/1.goose @@ -0,0 +1,129 @@ +// A permutation of up to 16 elements is one u64, element i in nibble i (as +// the SIMD programs keep it in one 128-bit register, a byte an element), each +// element stored complemented, as 15 - v. A flip of the first f + 1 elements +// moves bits between the word and a nibble-reversed copy of it, kept up to +// date alongside, with shifts by s = 4 * (15 - f): the low nibble of the +// complemented word times four. That makes s == 60 mean that 0 is first, and +// a "flip" with s == 60 moves nothing. Two permutations are flipped at once, +// until both are done, for the instruction-level parallelism and so that the +// unpredictable end of the flipping loop is paid once per pair. + +import std; +import os; + +let ONES: u64 = 0xFFFFFFFFFFFFFFFF; + +// HI[c]: the top 4c bits, which a flip with s = 4c leaves in place. +let HI: u64[16] = [0x0000000000000000, 0xF000000000000000, 0xFF00000000000000, 0xFFF0000000000000, + 0xFFFF000000000000, 0xFFFFF00000000000, 0xFFFFFF0000000000, 0xFFFFFFF000000000, + 0xFFFFFFFF00000000, 0xFFFFFFFFF0000000, 0xFFFFFFFFFF000000, 0xFFFFFFFFFFF00000, + 0xFFFFFFFFFFFF0000, 0xFFFFFFFFFFFFF000, 0xFFFFFFFFFFFFFF00, 0xFFFFFFFFFFFFFFF0]; + +fn nibble_reverse(x: u64) -> u64 { + var y = ((x >> 4) & 0x0F0F0F0F0F0F0F0F) | ((x & 0x0F0F0F0F0F0F0F0F) << 4); + y = ((y >> 8) & 0x00FF00FF00FF00FF) | ((y & 0x00FF00FF00FF00FF) << 8); + y = ((y >> 16) & 0x0000FFFF0000FFFF) | ((y & 0x0000FFFF0000FFFF) << 16); + return (y >> 32) | (y << 32); +} + +// The flip counts of the permutations a and b. +fn flips2(a: u64, b: u64) -> i64, i64 { + var xa = a; + var ra = nibble_reverse(a); + var sa = (a & 15) << 2; + var ka = 0; + var xb = b; + var rb = nibble_reverse(b); + var sb = (b & 15) << 2; + var kb = 0; + while (sa & sb) != 60 { + ka += if sa != 60 { 1 } else { 0 }; + kb += if sb != 60 { 1 } else { 0 }; + let ha = ra >> sa; + let ya = (xa & HI[sa >> 2]) | ha; + ra = (ra & ~(ONES << sa)) | (xa << sa); + xa = ya; + sa = (ha & 15) << 2; + let hb = rb >> sb; + let yb = (xb & HI[sb >> 2]) | hb; + rb = (rb & ~(ONES << sb)) | (xb << sb); + xb = yb; + sb = (hb & 15) << 2; + } + return ka, kb; +} + +// Rotates the first k elements left by one. +fn rotate(x: u64, k: i64) -> u64 { + let low = ONES >> (64 - (k << 2)); + return (x & ~low) | ((x >> 4) & (low >> 4)) | ((x & 15) << ((k - 1) << 2)); +} + +// The next permutation in the counting order: the first two elements swap, +// and each counter that wraps rotates one more element in. +fn next_perm(perm: u64, count: i64[16]&) -> u64 { + var x = (perm & ~0xFF) | ((perm & 15) << 4) | ((perm >> 4) & 15); + var i = 1; + while i < 15 { + count[i] += 1; + if count[i] <= i { break; } + count[i] = 0; + i++; + x = rotate(x, i + 1); + } + return x; +} + +fn factorial(n: i64) -> i64 { + var f = 1; + for i in 2..n + 1 { f *= i; } + return f; +} + +// The permutation with index `start`, and its counters. +fn start_perm(n: i64, start: i64, count: i64[16]&) -> u64 { + var p: u8[16] = [0; 16]; + for i in 16 { p[i] = i as u8; } + var idx = start; + var i = n - 1; + while i > 0 { + let f = factorial(i); + let d = idx / f; + idx %= f; + count[i] = d; + let t = p; + for j in i + 1 { p[j] = t[(j + d) % (i + 1)]; } + i--; + } + var perm: u64 = 0; + for j in 16 { perm |= (p[j] as u64) << (j << 2); } + return ~perm; +} + +// The checksum and the largest flip count of the permutations with indices +// start .. start + len, start and len even. +fn run_block(n: i64, start: i64, len: i64) -> i64, i64 { + var count: i64[16] = [0; 16]; + var perm = start_perm(n, start, count); + var checksum = 0; + var max_flips = 0; + for pair in len >> 1 { + let a = perm; + perm = next_perm(perm, count); + let b = perm; + perm = next_perm(perm, count); + let ka, kb = flips2(a, b); + checksum += ka - kb; + max_flips = max(max_flips, max(ka, kb)); + } + return checksum, max_flips; +} + +fn main() { + let a = args(); + var n = 7; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + assert(n >= 2 && n <= 16); + let checksum, max_flips = run_block(n, 0, factorial(n)); + print(checksum, "\nPfannkuchen(", n, ") = ", max_flips); +} diff --git a/bench/algorithm/fasta/1-m.goose b/bench/algorithm/fasta/1-m.goose new file mode 100644 index 000000000..3aa853700 --- /dev/null +++ b/bench/algorithm/fasta/1-m.goose @@ -0,0 +1,123 @@ +// The single-threaded program (see 1.goose) with the work split three ways: +// one worker makes the repeated sequence, another the two random ones, both +// in 256 KB blocks that come back through queues, and the main thread writes +// the blocks out as they arrive, so that writing overlaps with computing. +// The random numbers stay one sequence, made by one thread in order. + +import std; +import os; +import math; + +let LINE = 60; +let BLOCK = 1 << 18; +let IM = 139968; + +let ALU = "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGGGAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGACCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAATACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCAGCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGGAGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCCAGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA"; + +// The benchmark's naive linear congruential generator. +inline fn next_random(s: u64) -> u64 { (s * 3877 + 29573) % IM } + +// A block of output; the last of each part is marked. +struct Fixed { last: bool, data: u8[] } +struct Random { last: bool, data: u8[] } + +fn send_fixed(out: u8[>..]&) { qput(Fixed { last: false, data: copy(out) }); } +fn send_random(out: u8[>..]&) { qput(Random { last: false, data: copy(out) }); } + +// Send is one of the two above. +fn flush_if_full(out: u8[>..]&) { + if out.len >= BLOCK { + Send(out); + out.clear(); + } +} + +fn repeat_fasta(out: u8[>..]&, seq: u8[:], count: i64) { + // The sequence followed by a line's worth of its start, so that every + // line is one slice of it. + var padded: u8[>..] = []; + padded.append(seq); + padded.append(seq[..LINE]); + var off = 0; + var left = count; + while left > 0 { + let len = min(LINE, left); + out.append(padded[off..off + len]); + out.push('\n'); + off += len; + if off >= seq.len { off -= seq.len; } + left -= len; + flush_if_full(out, Send); + } +} + +// Appends count nucleotides of the distribution, returning the generator's +// state to go on from. +fn random_fasta(out: u8[>..]&, seed: u64, probs: f64[:], chars: u8[:], count: i64) -> u64 { + // Padded to 16 with thresholds no random number reaches. + var thresholds: i32[16] = [0x7FFFFFFF; 16]; + var table: u8[16] = [0; 16]; + var cum = 0.0; + for p, i in probs { + cum += p; + thresholds[i] = ceil(cum * IM) as i32; + table[i] = chars[i]; + } + var s = seed; + var left = count; + while left > 0 { + let len = min(LINE, left); + let start = out.len; + out.resize(start + len + 1, '\n'); + for &c in out[start..start + len] { + s = next_random(s); + let r = s as! i32; + var k = 0; + for t in thresholds { k += if t <= r { 1 } else { 0 }; } + c = table[k & 15]; + } + left -= len; + flush_if_full(out, Send); + } + return s; +} + +thread_fn fixed_part(n: i64) { + var out: u8[>..] = []; + out.append(">ONE Homo sapiens alu\n"); + repeat_fasta(out, ALU, 2 * n, send_fixed); + qput(Fixed { last: true, data: copy(out) }); +} + +thread_fn random_part(n: i64) { + var out: u8[>..] = []; + var seed: u64 = 42; + out.append(">TWO IUB ambiguity codes\n"); + seed = random_fasta(out, seed, [0.27, 0.12, 0.12, 0.27, 0.02, 0.02, 0.02, 0.02, + 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02], + "acgtBDHKMNRSVWY", 3 * n, send_random); + out.append(">THREE Homo sapiens frequency\n"); + seed = random_fasta(out, seed, [0.3029549426680, 0.1979883004921, 0.1975473066391, 0.3015094502008], + "acgt", 5 * n, send_random); + qput(Random { last: true, data: copy(out) }); +} + +fn main() { + let a = args(); + var n = 1000; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + let fixed = thread_spawn(fixed_part, n); + let random = thread_spawn(random_part, n); + loop { + let b = qget(); + write_stdout(b.data); + if b.last { break; } + } + loop { + let b = qget(); + write_stdout(b.data); + if b.last { break; } + } + thread_wait(fixed); + thread_wait(random); +} diff --git a/bench/algorithm/fasta/1.goose b/bench/algorithm/fasta/1.goose new file mode 100644 index 000000000..3c31eec6f --- /dev/null +++ b/bench/algorithm/fasta/1.goose @@ -0,0 +1,93 @@ +// The output is built in one buffer, written out and reused every 1 MB: +// stdout is unbuffered, so each write is a system call. The cumulative +// probabilities become integer thresholds, and each random number's +// nucleotide is the number of thresholds at or below it: a linear search +// without branches, which the C compiler does in vector registers. + +import std; +import os; +import math; + +let LINE = 60; +let FLUSH = 1 << 20; +let IM = 139968; + +let ALU = "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGGGAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGACCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAATACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCAGCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGGAGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCCAGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA"; + +// The benchmark's naive linear congruential generator. +inline fn next_random(s: u64) -> u64 { (s * 3877 + 29573) % IM } + +fn flush_if_full(out: u8[>..]&) { + if out.len >= FLUSH { + write_stdout(out); + out.clear(); + } +} + +fn repeat_fasta(out: u8[>..]&, seq: u8[:], count: i64) { + // The sequence followed by a line's worth of its start, so that every + // line is one slice of it. + var padded: u8[>..] = []; + padded.append(seq); + padded.append(seq[..LINE]); + var off = 0; + var left = count; + while left > 0 { + let len = min(LINE, left); + out.append(padded[off..off + len]); + out.push('\n'); + off += len; + if off >= seq.len { off -= seq.len; } + left -= len; + flush_if_full(out); + } +} + +// Appends count nucleotides of the distribution, returning the generator's +// state to go on from. +fn random_fasta(out: u8[>..]&, seed: u64, probs: f64[:], chars: u8[:], count: i64) -> u64 { + // Padded to 16 with thresholds no random number reaches. + var thresholds: i32[16] = [0x7FFFFFFF; 16]; + var table: u8[16] = [0; 16]; + var cum = 0.0; + for p, i in probs { + cum += p; + thresholds[i] = ceil(cum * IM) as i32; + table[i] = chars[i]; + } + var s = seed; + var left = count; + while left > 0 { + let len = min(LINE, left); + let start = out.len; + out.resize(start + len + 1, '\n'); + for &c in out[start..start + len] { + s = next_random(s); + let r = s as! i32; + var k = 0; + for t in thresholds { k += if t <= r { 1 } else { 0 }; } + c = table[k & 15]; + } + left -= len; + flush_if_full(out); + } + return s; +} + +fn main() { + let a = args(); + var n = 1000; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + var out: u8[>..] = []; + out.append(">ONE Homo sapiens alu\n"); + repeat_fasta(out, ALU, 2 * n); + var seed: u64 = 42; + out.append(">TWO IUB ambiguity codes\n"); + seed = random_fasta(out, seed, [0.27, 0.12, 0.12, 0.27, 0.02, 0.02, 0.02, 0.02, + 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02], + "acgtBDHKMNRSVWY", 3 * n); + out.append(">THREE Homo sapiens frequency\n"); + seed = random_fasta(out, seed, [0.3029549426680, 0.1979883004921, 0.1975473066391, 0.3015094502008], + "acgt", 5 * n); + write_stdout(out); +} diff --git a/bench/algorithm/helloworld/1.goose b/bench/algorithm/helloworld/1.goose new file mode 100644 index 000000000..a3297aecf --- /dev/null +++ b/bench/algorithm/helloworld/1.goose @@ -0,0 +1,6 @@ +import os; + +fn main() { + let a = args(); + print("Hello world ", if a.len > 1 { a[1] } else { "" }, "!"); +} diff --git a/bench/algorithm/json-serde/1.goose b/bench/algorithm/json-serde/1.goose new file mode 100644 index 000000000..a1ed2d6f2 --- /dev/null +++ b/bench/algorithm/json-serde/1.goose @@ -0,0 +1,255 @@ +// The documents are read with json's reader straight into typed records, as +// the serde derives of the Rust programs read them, and written back with +// json's writer. Records are fixed-size and kept in arenas: a feature's +// polygon is a run of rings, a ring a run of points, and a document a run of +// features, each run a start and a count in the arena below it. Strings are +// views of the text, or of the decoded text for those that held escapes. + +import std; +import os; +import json; + +struct Run { first: i64, count: i64 } +struct Properties { name: const u8[:] } +struct Geometry { kind: const u8[:], coordinates: Run } // rings +struct Feature { kind: const u8[:], properties: Properties, geometry: Geometry } +struct GeoData { kind: const u8[:], features: Run } + +// --- numbers -------------------------------------------------------------------- + +// As the Rust programs' MyF64 writes a number: integral values as an i64 +// (saturated, as Rust's `as` is), others as json writes them. +fn format_num(out: u8[>..]&, x: f64) { + if x != x || x - x != 0.0 { json::format_number(out, x); return; } + if x >= 9223372036854775807.0 { format(out, 9223372036854775807); return; } + if x <= -9223372036854775808.0 { format(out, -9223372036854775807 - 1); return; } + let i = x as! i64; + if i == x { format(out, i); return; } + json::format_number(out, x); +} + +// --- MD5 (RFC 1321) --------------------------------------------------------------- + +// G's two terms have no bits in common, so they can be added, which lets +// the one without b be summed before b is known. +inline fn rotl(x: u32, s: u32) -> u32 { (x << s) | (x >> (32 - s)) } +inline fn ff(a: u32, b: u32, c: u32, d: u32, m: u32, s: u32, k: u32) -> u32 { b + rotl(a + (d ^ (b & (c ^ d))) + m + k, s) } +inline fn gg(a: u32, b: u32, c: u32, d: u32, m: u32, s: u32, k: u32) -> u32 { b + rotl(a + m + k + (c & ~d) + (b & d), s) } +inline fn hh(a: u32, b: u32, c: u32, d: u32, m: u32, s: u32, k: u32) -> u32 { b + rotl(a + (b ^ c ^ d) + m + k, s) } +inline fn ii(a: u32, b: u32, c: u32, d: u32, m: u32, s: u32, k: u32) -> u32 { b + rotl(a + (c ^ (b | ~d)) + m + k, s) } + +// One 64-byte block into the state h. +fn md5_block(h: u32[4]&, blk: u8[:]) { + assert(blk.len == 64); + var m: u32[16] = [0; 16]; + for j in 16 { + m[j] = blk[4 * j] | (blk[4 * j + 1] as u32) << 8 | (blk[4 * j + 2] as u32) << 16 | (blk[4 * j + 3] as u32) << 24; + } + var a = h[0]; + var b = h[1]; + var c = h[2]; + var d = h[3]; + a = ff(a, b, c, d, m[0], 7, 0xd76aa478); d = ff(d, a, b, c, m[1], 12, 0xe8c7b756); + c = ff(c, d, a, b, m[2], 17, 0x242070db); b = ff(b, c, d, a, m[3], 22, 0xc1bdceee); + a = ff(a, b, c, d, m[4], 7, 0xf57c0faf); d = ff(d, a, b, c, m[5], 12, 0x4787c62a); + c = ff(c, d, a, b, m[6], 17, 0xa8304613); b = ff(b, c, d, a, m[7], 22, 0xfd469501); + a = ff(a, b, c, d, m[8], 7, 0x698098d8); d = ff(d, a, b, c, m[9], 12, 0x8b44f7af); + c = ff(c, d, a, b, m[10], 17, 0xffff5bb1); b = ff(b, c, d, a, m[11], 22, 0x895cd7be); + a = ff(a, b, c, d, m[12], 7, 0x6b901122); d = ff(d, a, b, c, m[13], 12, 0xfd987193); + c = ff(c, d, a, b, m[14], 17, 0xa679438e); b = ff(b, c, d, a, m[15], 22, 0x49b40821); + + a = gg(a, b, c, d, m[1], 5, 0xf61e2562); d = gg(d, a, b, c, m[6], 9, 0xc040b340); + c = gg(c, d, a, b, m[11], 14, 0x265e5a51); b = gg(b, c, d, a, m[0], 20, 0xe9b6c7aa); + a = gg(a, b, c, d, m[5], 5, 0xd62f105d); d = gg(d, a, b, c, m[10], 9, 0x02441453); + c = gg(c, d, a, b, m[15], 14, 0xd8a1e681); b = gg(b, c, d, a, m[4], 20, 0xe7d3fbc8); + a = gg(a, b, c, d, m[9], 5, 0x21e1cde6); d = gg(d, a, b, c, m[14], 9, 0xc33707d6); + c = gg(c, d, a, b, m[3], 14, 0xf4d50d87); b = gg(b, c, d, a, m[8], 20, 0x455a14ed); + a = gg(a, b, c, d, m[13], 5, 0xa9e3e905); d = gg(d, a, b, c, m[2], 9, 0xfcefa3f8); + c = gg(c, d, a, b, m[7], 14, 0x676f02d9); b = gg(b, c, d, a, m[12], 20, 0x8d2a4c8a); + + a = hh(a, b, c, d, m[5], 4, 0xfffa3942); d = hh(d, a, b, c, m[8], 11, 0x8771f681); + c = hh(c, d, a, b, m[11], 16, 0x6d9d6122); b = hh(b, c, d, a, m[14], 23, 0xfde5380c); + a = hh(a, b, c, d, m[1], 4, 0xa4beea44); d = hh(d, a, b, c, m[4], 11, 0x4bdecfa9); + c = hh(c, d, a, b, m[7], 16, 0xf6bb4b60); b = hh(b, c, d, a, m[10], 23, 0xbebfbc70); + a = hh(a, b, c, d, m[13], 4, 0x289b7ec6); d = hh(d, a, b, c, m[0], 11, 0xeaa127fa); + c = hh(c, d, a, b, m[3], 16, 0xd4ef3085); b = hh(b, c, d, a, m[6], 23, 0x04881d05); + a = hh(a, b, c, d, m[9], 4, 0xd9d4d039); d = hh(d, a, b, c, m[12], 11, 0xe6db99e5); + c = hh(c, d, a, b, m[15], 16, 0x1fa27cf8); b = hh(b, c, d, a, m[2], 23, 0xc4ac5665); + + a = ii(a, b, c, d, m[0], 6, 0xf4292244); d = ii(d, a, b, c, m[7], 10, 0x432aff97); + c = ii(c, d, a, b, m[14], 15, 0xab9423a7); b = ii(b, c, d, a, m[5], 21, 0xfc93a039); + a = ii(a, b, c, d, m[12], 6, 0x655b59c3); d = ii(d, a, b, c, m[3], 10, 0x8f0ccc92); + c = ii(c, d, a, b, m[10], 15, 0xffeff47d); b = ii(b, c, d, a, m[1], 21, 0x85845dd1); + a = ii(a, b, c, d, m[8], 6, 0x6fa87e4f); d = ii(d, a, b, c, m[15], 10, 0xfe2ce6e0); + c = ii(c, d, a, b, m[6], 15, 0xa3014314); b = ii(b, c, d, a, m[13], 21, 0x4e0811a1); + a = ii(a, b, c, d, m[4], 6, 0xf7537e82); d = ii(d, a, b, c, m[11], 10, 0xbd3af235); + c = ii(c, d, a, b, m[2], 15, 0x2ad7d2bb); b = ii(b, c, d, a, m[9], 21, 0xeb86d391); + + h[0] += a; + h[1] += b; + h[2] += c; + h[3] += d; +} + +// The digest of data, in lower-case hex. +fn md5_hex(data: u8[:]) -> u8[>..] { + var h: u32[4] = [0x67452301, 0xefcdab89, 0x98badcfe, 0x10325476]; + let full = data.len / 64 * 64; + var at = 0; + while at < full { + md5_block(h, data[at..at + 64]); + at += 64; + } + // The rest, a 0x80, zeros up to 8 bytes short of a block, the length in bits. + var tail: u8[128] = [0; 128]; + let rest = data.len - full; + copy_into(tail[..rest], data[full..]); + tail[rest] = 0x80; + let end = if rest < 56 { 64 } else { 128 }; + let bits = (data.len as u64) * 8; + for j in 8 { tail[end - 8 + j] = (bits >> (8 * j)) as! u8; } + md5_block(h, tail[..64]); + if end == 128 { md5_block(h, tail[64..]); } + var out: u8[>..] = []; + for w in h { + for j in 4 { format_uint(out, (w >> (8 * j)) & 0xff, 16, 2, '0'); } + } + return out; +} + +// --- main ------------------------------------------------------------------------- + +fn main() { + let argv = args(); + var path: u8[>..] = []; + path.append(if argv.len > 1 { argv[1] } else { "sample" }); + path.append(".json"); + var n = 10; + if argv.len > 2 { let v, ok = parse_int(argv[2]); if ok { n = v; } } + var text: u8[>..] = []; + if !read_file(path, text) { abort("cannot read the input"); } + + var decoded: u8[>..] = []; // strings that held escapes + // The strings read are views of the reader's text, so the reader is + // one that outlives the records; each document resets it. + var r = json::reader(text); + var points: f64[>..] = []; // x and y, by pairs + var rings: Run[>..] = []; // runs of points + var features: Feature[>..] = []; + var docs: GeoData[>..] = []; + + fn read_ring(r: json::Reader&) -> Run { + let first = points.len / 2; + json::each_element(r) { + var k = 0; + json::each_element(r) { + if k < 2 { points.push(json::read_number(r)); } else { json::skip(r); } + k++; + }; + if k != 2 { abort("a point is two numbers"); } + }; + return Run { first: first, count: points.len / 2 - first }; + } + + fn read_geometry(r: json::Reader&) -> Geometry { + var g = Geometry { kind: "", coordinates: Run { 0, 0 } }; + json::each_member(r) { key => + if key == "type" { + g.kind = json::read_string(r, decoded); + } else if key == "coordinates" { + let first = rings.len; + json::each_element(r) { rings.push(read_ring(r)); }; + g.coordinates = Run { first: first, count: rings.len - first }; + } else { + json::skip(r); + } + }; + return g; + } + + fn read_feature(r: json::Reader&) -> Feature { + var f = Feature { kind: "", properties: Properties { name: "" }, geometry: Geometry { kind: "", coordinates: Run { 0, 0 } } }; + json::each_member(r) { key => + if key == "type" { + f.kind = json::read_string(r, decoded); + } else if key == "properties" { + json::each_member(r) { pkey => + if pkey == "name" { f.properties.name = json::read_string(r, decoded); } else { json::skip(r); } + }; + } else if key == "geometry" { + f.geometry = read_geometry(r); + } else { + json::skip(r); + } + }; + return f; + } + + fn read_doc() -> GeoData { + r = json::reader(text); + var d = GeoData { kind: "", features: Run { 0, 0 } }; + json::each_member(r) { key => + if key == "type" { + d.kind = json::read_string(r, decoded); + } else if key == "features" { + let first = features.len; + json::each_element(r) { features.push(read_feature(r)); }; + d.features = Run { first: first, count: features.len - first }; + } else { + json::skip(r); + } + }; + if !json::finish(r) { abort(json::message(r)); } + return d; + } + + fn write_doc(out: u8[>..]&, d: GeoData) { + out.append("{\"type\":"); + json::format_string(out, d.kind); + out.append(",\"features\":["); + for fi in d.features.count { + let f = features[d.features.first + fi]; + if fi > 0 { out.push(','); } + out.append("{\"type\":"); + json::format_string(out, f.kind); + out.append(",\"properties\":{\"name\":"); + json::format_string(out, f.properties.name); + out.append("},\"geometry\":{\"type\":"); + json::format_string(out, f.geometry.kind); + out.append(",\"coordinates\":["); + let c = f.geometry.coordinates; + for ri in c.count { + let ring = rings[c.first + ri]; + if ri > 0 { out.push(','); } + out.push('['); + for pi in ring.count { + let at = 2 * (ring.first + pi); + if pi > 0 { out.push(','); } + out.push('['); + format_num(out, points[at]); + out.push(','); + format_num(out, points[at + 1]); + out.push(']'); + } + out.push(']'); + } + out.append("]}}"); + } + out.append("]}"); + } + + var out: u8[>..] = []; + let first = read_doc(); + write_doc(out, first); + print(md5_hex(out)); + + for i in n { docs.push(read_doc()); } + var all: u8[>..] = []; + all.push('['); + for d, i in docs { + if i > 0 { all.push(','); } + write_doc(all, d); + } + all.push(']'); + print(md5_hex(all)); +} diff --git a/bench/algorithm/knucleotide/1-m.goose b/bench/algorithm/knucleotide/1-m.goose new file mode 100644 index 000000000..f8f56a22a --- /dev/null +++ b/bench/algorithm/knucleotide/1-m.goose @@ -0,0 +1,215 @@ +// The single-threaded program's reading and counting (see 1.goose), with +// each k counted by a worker thread of its own: the workers take jobs from a +// queue, the largest first, and send back the text of their part of the +// output. Each worker gets a copy of the coded sequence as it starts, packed +// four codes a byte to keep the copies small. + +import std; +import os; + +// A slot holds 0 when it is unused, else the key, a set bit and the count in +// the low COUNT_BITS bits: one compare of the bits above the count tests the +// key, and the count increments in place. +let COUNT_BITS = 27; +let COUNT_MASK: u64 = (1 << COUNT_BITS) - 1; +let CHUNK = 1 << 14; + +fn home(key: u64) -> i64 { (key ^ (key >> 7)) as! i64 } + +fn code_of(c: u8) -> u8 { (c >> 1) & 3 } + +// The codes of the sequence after the line starting with ">THREE", written +// over the start of text. +fn read_three(text: u8[:]) -> u8[:] { + let header = byte_set(">"); + var at = 0; + loop { + let k = find_any(text[at..], header); + assert(k >= 0); + at += k; + if starts_with(text[at..], ">THREE") { break; } + at += 1; + } + at += find(text[at..], "\n") + 1; + var end = find_any(text[at..], header); + end = if end < 0 { text.len } else { at + end }; + // Every byte is stored, and the next one overwrites a newline's. + var n = 0; + for c in text[at..end] { + text[n] = code_of(c); + n += if c != '\n' { 1 } else { 0 }; + } + return text[..n]; +} + +// Doubles the table, putting every entry in again. +fn grow(slots: u64[>..]&) { + var old: u64[>..] = []; + old.append(slots); + slots.clear(); + slots.resize(old.len * 2, 0); + let mask = slots.len - 1; + for v in old { + guard v != 0; + var i = home(v >> (COUNT_BITS + 1)) & mask; + while slots[i] != 0 { i = (i + 1) & mask; } + slots[i] = v; + } +} + +// Four codes a byte, the first in the top bits. +fn pack(seq: u8[:]) -> u8[>..] { + var packed: u8[>..] = []; + packed.resize((seq.len + 3) / 4, 0); + for c, i in seq { packed[i >> 2] |= c << (6 - 2 * (i & 3)); } + return packed; +} + +// The table of the k-nucleotides of the n codes packed, which are unpacked +// a chunk at a time into a buffer that stays in the cache. +fn count(packed: u8[:], n: i64, k: i64) -> u64[>..] { + var slots: u64[>..] = []; + slots.resize(16, 0); + var used = 0; + var chunk: u8[>..] = []; + chunk.resize(CHUNK, 0); + let key_mask: u64 = (1 << (2 * k)) - 1; + // The codes so far, the last k of them the key (see 1.goose). + var codes: u64 = 0; + var pos = 0; + while pos < n { + let m = min(CHUNK, n - pos); + for b, i in packed[pos >> 2..(pos + m + 3) >> 2] { + let d = chunk[4 * i..4 * i + 4]; + d[0] = b >> 6; + d[1] = (b >> 4) & 3; + d[2] = (b >> 2) & 3; + d[3] = b & 3; + } + var from = 0; + if pos == 0 { + for c in chunk[..k - 1] { codes = (codes << 2) | c; } + from = k - 1; + } + for c in chunk[from..m] { + codes = (codes << 2) | c; + let key = codes & key_mask; + let tag = (key << 1) | 1; + let mask = slots.len - 1; + var i = home(key); + loop { + let s = i & mask; + let v = slots[s]; + if v >> COUNT_BITS == tag { slots[s] = v + 1; break; } + if v == 0 { + slots[s] = (tag << COUNT_BITS) | 1; + used++; + if used * 2 > slots.len { grow(slots); } + break; + } + i++; + } + } + pos += m; + } + return slots; +} + +fn decode(out: u8[>..]&, key: u64, k: i64) { + for i in k { + out.push("ACTG"[((key >> (2 * (k - 1 - i))) & 3) as! i64]); + } +} + +struct Entry { key: u64, count: i64 } + +fn frequencies(out: u8[>..]&, packed: u8[:], n: i64, k: i64) { + let slots = count(packed, n, k); + var entries: Entry[>..] = []; + var total = 0; + for v in slots { + guard v != 0; + let n = (v & COUNT_MASK) as! i64; + entries.push(Entry { v >> (COUNT_BITS + 1), n }); + total += n; + } + sort(entries) { a, b => a.count > b.count || (a.count == b.count && a.key > b.key) }; + for e in entries { + decode(out, e.key, k); + out.push(' '); + format_flt(out, (e.count as f32 * 100.0) / total, 3); + out.push('\n'); + } + out.push('\n'); +} + +fn occurrences(out: u8[>..]&, packed: u8[:], n: i64, name: u8[:]) { + let slots = count(packed, n, name.len); + var key: u64 = 0; + for c in name { key = (key << 2) | code_of(c); } + let tag = (key << 1) | 1; + var n: u64 = 0; + let mask = slots.len - 1; + var i = home(key) & mask; + while slots[i] != 0 { + if slots[i] >> COUNT_BITS == tag { n = slots[i] & COUNT_MASK; break; } + i = (i + 1) & mask; + } + format(out, n, "\t", name, "\n"); +} + +// A part of the output: the frequencies of the k-nucleotides where name is +// empty, else the count of name; k == 0 tells a worker to stop. +struct Job { part: i64, k: i64, name: u8[..18] } +struct Report { part: i64, text: u8[] } + +thread_fn worker(packed: u8[], n: i64) { + loop { + let job = qget(); + guard job.k > 0 else { break; } + var out: u8[>..] = []; + if job.name.len == 0 { + frequencies(out, packed, n, job.k); + } else { + occurrences(out, packed, n, job.name); + } + qput(Report { part: job.part, text: copy(out) }); + } +} + +fn main() { + let a = args(); + assert(a.len > 1); + var text: u8[>..] = []; + if !read_file(a[1], text) { abort("cannot read the input"); } + let seq = read_three(text); + let packed = pack(seq); + let names = ["GGT", "GGTA", "GGTATT", "GGTATTTTAATT", "GGTATTTTAATTTATAGT"]; + let parts = 2 + names.len; + for i in names.len { + let name = names[names.len - 1 - i]; + qput(Job { part: parts - 1 - i, k: name.len, name: name }); + } + qput(Job { part: 1, k: 2, name: "" }); + qput(Job { part: 0, k: 1, name: "" }); + let workers = clamp(hardware_threads(), 1, parts); + for i in workers { qput(Job { part: 0, k: 0, name: "" }); } + var ids: i64[>..] = []; + for i in workers { ids.push(thread_spawn(worker, packed[..], seq.len)); } + // The parts arrive in any order; at[p] is where part p's text went. + var reports: u8[>..] = []; + var at: i64[>..] = []; + at.resize(parts, 0); + var len: i64[>..] = []; + len.resize(parts, 0); + for i in parts { + let r = qget(); + at[r.part] = reports.len; + len[r.part] = r.text.len; + reports.append(r.text); + } + for id in ids { thread_wait(id); } + var out: u8[>..] = []; + for p in parts { out.append(reports[at[p]..at[p] + len[p]]); } + write_stdout(out); +} diff --git a/bench/algorithm/knucleotide/1.goose b/bench/algorithm/knucleotide/1.goose new file mode 100644 index 000000000..1e9dbf583 --- /dev/null +++ b/bench/algorithm/knucleotide/1.goose @@ -0,0 +1,153 @@ +// The third sequence is read into memory once and coded two bits a +// nucleotide (A 0, C 1, T 2, G 3: bits 1 and 2 of the letter, either case), +// in place over the file's own bytes. Each k-nucleotide is then the last k +// codes packed into an integer, counted in a hash table of this program's +// own, as C's program has: open addressing and linear probing over a +// power-of-two array of u64 slots, each holding a key and its count, which +// starts at 16 slots and doubles at half load. Half rather than the usual +// three quarters, as most lookups then end at the first slot they try, +// which spares a mispredicted branch on most of the others. (2.goose counts +// with the library's dictionary instead.) + +import std; +import os; + +// A slot holds 0 when it is unused, else the key, a set bit and the count in +// the low COUNT_BITS bits: one compare of the bits above the count tests the +// key, and the count increments in place. +let COUNT_BITS = 27; +let COUNT_MASK: u64 = (1 << COUNT_BITS) - 1; + +fn home(key: u64) -> i64 { (key ^ (key >> 7)) as! i64 } + +fn code_of(c: u8) -> u8 { (c >> 1) & 3 } + +// The codes of the sequence after the line starting with ">THREE", written +// over the start of text. +fn read_three(text: u8[:]) -> u8[:] { + let header = byte_set(">"); + var at = 0; + loop { + let k = find_any(text[at..], header); + assert(k >= 0); + at += k; + if starts_with(text[at..], ">THREE") { break; } + at += 1; + } + at += find(text[at..], "\n") + 1; + var end = find_any(text[at..], header); + end = if end < 0 { text.len } else { at + end }; + // Every byte is stored, and the next one overwrites a newline's. + var n = 0; + for c in text[at..end] { + text[n] = code_of(c); + n += if c != '\n' { 1 } else { 0 }; + } + return text[..n]; +} + +// Doubles the table, putting every entry in again. +fn grow(slots: u64[>..]&) { + var old: u64[>..] = []; + old.append(slots); + slots.clear(); + slots.resize(old.len * 2, 0); + let mask = slots.len - 1; + for v in old { + guard v != 0; + var i = home(v >> (COUNT_BITS + 1)) & mask; + while slots[i] != 0 { i = (i + 1) & mask; } + slots[i] = v; + } +} + +// The table of the k-nucleotides of seq. +fn count(seq: u8[:], k: i64) -> u64[>..] { + var slots: u64[>..] = []; + slots.resize(16, 0); + var used = 0; + let key_mask: u64 = (1 << (2 * k)) - 1; + // The codes so far, the last k of them the key: shifting the earlier + // ones out of the word rather than masking them off keeps the masking + // out of the chain from one position to the next. + var codes: u64 = 0; + for c in seq[..k - 1] { codes = (codes << 2) | c; } + for c in seq[k - 1..] { + codes = (codes << 2) | c; + let key = codes & key_mask; + let tag = (key << 1) | 1; + let mask = slots.len - 1; + var i = home(key); + loop { + let s = i & mask; + let v = slots[s]; + if v >> COUNT_BITS == tag { slots[s] = v + 1; break; } + if v == 0 { + slots[s] = (tag << COUNT_BITS) | 1; + used++; + if used * 2 > slots.len { grow(slots); } + break; + } + i++; + } + } + return slots; +} + +fn decode(out: u8[>..]&, key: u64, k: i64) { + for i in k { + out.push("ACTG"[((key >> (2 * (k - 1 - i))) & 3) as! i64]); + } +} + +struct Entry { key: u64, count: i64 } + +fn frequencies(out: u8[>..]&, seq: u8[:], k: i64) { + let slots = count(seq, k); + var entries: Entry[>..] = []; + var total = 0; + for v in slots { + guard v != 0; + let n = (v & COUNT_MASK) as! i64; + entries.push(Entry { v >> (COUNT_BITS + 1), n }); + total += n; + } + sort(entries) { a, b => a.count > b.count || (a.count == b.count && a.key > b.key) }; + for e in entries { + decode(out, e.key, k); + out.push(' '); + format_flt(out, (e.count as f32 * 100.0) / total, 3); + out.push('\n'); + } + out.push('\n'); +} + +fn occurrences(out: u8[>..]&, seq: u8[:], name: u8[:]) { + let slots = count(seq, name.len); + var key: u64 = 0; + for c in name { key = (key << 2) | code_of(c); } + let tag = (key << 1) | 1; + var n: u64 = 0; + let mask = slots.len - 1; + var i = home(key) & mask; + while slots[i] != 0 { + if slots[i] >> COUNT_BITS == tag { n = slots[i] & COUNT_MASK; break; } + i = (i + 1) & mask; + } + format(out, n, "\t", name, "\n"); +} + +fn main() { + let a = args(); + assert(a.len > 1); + var text: u8[>..] = []; + if !read_file(a[1], text) { abort("cannot read the input"); } + let seq = read_three(text); + var out: u8[>..] = []; + frequencies(out, seq, 1); + frequencies(out, seq, 2); + for name in ["GGT", "GGTA", "GGTATT", "GGTATTTTAATT", "GGTATTTTAATTTATAGT"] { + occurrences(out, seq, name); + } + write_stdout(out); +} diff --git a/bench/algorithm/knucleotide/2.goose b/bench/algorithm/knucleotide/2.goose new file mode 100644 index 000000000..d58b574da --- /dev/null +++ b/bench/algorithm/knucleotide/2.goose @@ -0,0 +1,100 @@ +// The third sequence is read into memory once and coded two bits a +// nucleotide (A 0, C 1, T 2, G 3: bits 1 and 2 of the letter, either case), +// in place over the file's own bytes. Each k-nucleotide is then the last k +// codes packed into a u64, the key of the library's dictionary counting +// them. (1.goose counts with a hash table of its own.) + +import std; +import os; +import dictionary; + +// The key, with a hash of its own: std's hash of a u64 scatters the small +// dense keys of short k-nucleotides as it would any, into the collisions of +// a random hash, where this one keeps them apart and is cheaper. +struct Code { v: u64 } +fn hash(c: Code) -> u64 { c.v ^ (c.v >> 7) } + +fn code_of(c: u8) -> u8 { (c >> 1) & 3 } + +// The codes of the sequence after the line starting with ">THREE", written +// over the start of text. +fn read_three(text: u8[:]) -> u8[:] { + let header = byte_set(">"); + var at = 0; + loop { + let k = find_any(text[at..], header); + assert(k >= 0); + at += k; + if starts_with(text[at..], ">THREE") { break; } + at += 1; + } + at += find(text[at..], "\n") + 1; + var end = find_any(text[at..], header); + end = if end < 0 { text.len } else { at + end }; + // Every byte is stored, and the next one overwrites a newline's. + var n = 0; + for c in text[at..end] { + text[n] = code_of(c); + n += if c != '\n' { 1 } else { 0 }; + } + return text[..n]; +} + +fn count(seq: u8[:], k: i64) -> dictionary { + var counts = dictionary {}; + let mask: u64 = (1 << (2 * k)) - 1; + var codes: u64 = 0; + for c in seq[..k - 1] { codes = (codes << 2) | c; } + for c in seq[k - 1..] { + codes = (codes << 2) | c; + let key = codes & mask; + counts.update(Code { key }, 0) { it += 1; }; + } + return counts; +} + +fn decode(out: u8[>..]&, key: u64, k: i64) { + for i in k { + out.push("ACTG"[((key >> (2 * (k - 1 - i))) & 3) as! i64]); + } +} + +struct Entry { key: u64, count: i32 } + +fn frequencies(out: u8[>..]&, seq: u8[:], k: i64) { + var counts = count(seq, k); + var entries: Entry[>..] = []; + var total = 0; + counts.each() { key, n => entries.push(Entry { key.v, n }); total += n; }; + sort(entries) { a, b => a.count > b.count || (a.count == b.count && a.key > b.key) }; + for e in entries { + decode(out, e.key, k); + out.push(' '); + format_flt(out, (e.count as f32 * 100.0) / total, 3); + out.push('\n'); + } + out.push('\n'); +} + +fn occurrences(out: u8[>..]&, seq: u8[:], name: u8[:]) { + var counts = count(seq, name.len); + var key: u64 = 0; + for c in name { key = (key << 2) | code_of(c); } + let n = counts.get(Code { key }); + format(out, if n { n } else { 0 }, "\t", name, "\n"); +} + +fn main() { + let a = args(); + assert(a.len > 1); + var text: u8[>..] = []; + if !read_file(a[1], text) { abort("cannot read the input"); } + let seq = read_three(text); + var out: u8[>..] = []; + frequencies(out, seq, 1); + frequencies(out, seq, 2); + for name in ["GGT", "GGTA", "GGTATT", "GGTATTTTAATT", "GGTATTTTAATTTATAGT"] { + occurrences(out, seq, name); + } + write_stdout(out); +} diff --git a/bench/algorithm/lru/1.goose b/bench/algorithm/lru/1.goose new file mode 100644 index 000000000..2f75bf091 --- /dev/null +++ b/bench/algorithm/lru/1.goose @@ -0,0 +1,91 @@ +// The cache's entries are 16-byte nodes in one pool, on a circular list +// through a sentinel, least recently used first; the links, and the +// dictionary's values, are 4-byte offsets from the pool's base. A full +// cache evicts by reusing its least recently used node for the new key. + +import std; +import os; +import dictionary; + +struct Node { key: u32, val: u32, prev: Node&, next: Node& } + +reusable var pool: Node[>..] = []; +let head = pool.alloc_index(Node { key: 0, val: 0, prev: self, next: self }); + +struct Lcg { seed: u32 } + +fn next(r: Lcg&) -> u32 { + r.seed = (1103515245 * r.seed + 12345) % 2147483648; + return r.seed; +} + +inline fn unlink(n: Node&) { + n.prev.next .= n.next; + n.next.prev .= n.prev; +} + +// n becomes the most recently used. +inline fn link_last(n: Node&) { + let h .= pool[head]; + n.prev .= h.prev; + n.next .= h; + h.prev.next .= n; + h.prev .= n; +} + +fn main() { + let a = args(); + var size = 100; + var n = 1000; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { size = v; } } + if a.len > 2 { let v, ok = parse_int(a[2]); if ok { n = v; } } + let modular = (size * 10) as u32; + var map = dictionary?> {}; + map.reserve(size); + var count = 0; + + fn cache_put(key: u32, val: u32) { + let f = map.get(key); + if f { + f.val = val; + unlink(f); + link_last(f); + return; + } + if count == size { + let lru .= pool[head].next; + map.remove(lru.key); + lru.key = key; + lru.val = val; + unlink(lru); + link_last(lru); + map.insert(key, lru); + return; + } + let e .= pool.alloc_ref(Node { key: key, val: val, prev: pool[head], next: pool[head] }); + link_last(e); + map.insert(key, e); + count++; + } + + fn cache_get(key: u32) -> bool { + let f = map.get(key); + guard f else { return false; } + unlink(f); + link_last(f); + return true; + } + + var rng0 = Lcg { 0 }; + var rng1 = Lcg { 1 }; + var hit = 0; + var missed = 0; + for i in n { + let n0 = rng0.next() % modular; + cache_put(n0, n0); + let n1 = rng1.next() % modular; + if cache_get(n1) { hit++; } else { missed++; } + } + print(hit); + print(missed); +} diff --git a/bench/algorithm/lru/2.goose b/bench/algorithm/lru/2.goose new file mode 100644 index 000000000..cd3ac0008 --- /dev/null +++ b/bench/algorithm/lru/2.goose @@ -0,0 +1,136 @@ +// 1.goose with the dictionary replaced by an open-addressing table made for +// the cache: 8-byte slots holding a key and a 4-byte link into the node +// pool (null for an empty slot), a multiplicative hash, linear probing and +// backward-shift deletion. Nine lookups in ten miss, and a miss probes to +// the end of its run, so the table is kept at most a quarter full: 32 KB +// for 1000 entries. At half full the program takes twice as long. + +import std; +import os; + +struct Node { key: u32, val: u32, prev: Node&, next: Node& } +struct Slot { key: u32, node: Node&? } + +reusable var pool: Node[>..] = []; +let head = pool.alloc_index(Node { key: 0, val: 0, prev: self, next: self }); + +struct Lcg { seed: u32 } + +fn next(r: Lcg&) -> u32 { + r.seed = (1103515245 * r.seed + 12345) % 2147483648; + return r.seed; +} + +inline fn unlink(n: Node&) { + n.prev.next .= n.next; + n.next.prev .= n.prev; +} + +// n becomes the most recently used. +inline fn link_last(n: Node&) { + let h .= pool[head]; + n.prev .= h.prev; + n.next .= h; + h.prev.next .= n; + h.prev .= n; +} + +fn main() { + let a = args(); + var size = 100; + var n = 1000; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { size = v; } } + if a.len > 2 { let v, ok = parse_int(a[2]); if ok { n = v; } } + let modular = (size * 10) as u32; + + var bits: u32 = 3; + while (1 << bits) < 4 * size { bits++; } + let shift = 32 - bits; + var slots: Slot[>..] = []; + push_n(slots, Slot { key: 0, node: null }, 1 << bits); + + // Slot indices run on past the end of the table and wrap here. + fn slot(i: i64) -> Slot& { slots[i & (slots.len - 1)] } + fn home(key: u32) -> i64 { (key * 2654435769) >> shift } + + // The slot holding key, or the empty slot where it goes. + fn find(key: u32) -> i64 { + var i = home(key); + loop { + let s = slot(i); + guard s.node && s.key != key else { return i; } + i++; + } + } + + // Empties slot `at`, moving the entries after it in its run that may + // move back into the hole. + fn remove_at(at: i64) { + var i = at; + var j = at; + loop { + j++; + let s = slot(j); + guard s.node else { break; } + if ((j - home(s.key)) & (slots.len - 1)) >= j - i { + let hole .= slot(i); + hole = s; + i = j; + } + } + slot(i).node .= null; + } + + var count = 0; + + fn cache_put(key: u32, val: u32) { + let at = find(key); + let f = slot(at).node; + if f { + f.val = val; + unlink(f); + link_last(f); + return; + } + if count == size { + // The new key goes into the empty slot found for it before the + // old one's slot is emptied, which may move it. + let lru .= pool[head].next; + let old = find(lru.key); + lru.key = key; + lru.val = val; + unlink(lru); + link_last(lru); + let s .= slot(at); + s = Slot { key: key, node: lru }; + remove_at(old); + return; + } + let e .= pool.alloc_ref(Node { key: key, val: val, prev: pool[head], next: pool[head] }); + link_last(e); + let s .= slot(at); + s = Slot { key: key, node: e }; + count++; + } + + fn cache_get(key: u32) -> bool { + let f = slot(find(key)).node; + guard f else { return false; } + unlink(f); + link_last(f); + return true; + } + + var rng0 = Lcg { 0 }; + var rng1 = Lcg { 1 }; + var hit = 0; + var missed = 0; + for i in n { + let n0 = rng0.next() % modular; + cache_put(n0, n0); + let n1 = rng1.next() % modular; + if cache_get(n1) { hit++; } else { missed++; } + } + print(hit); + print(missed); +} diff --git a/bench/algorithm/mandelbrot/1-m.goose b/bench/algorithm/mandelbrot/1-m.goose new file mode 100644 index 000000000..2f516048c --- /dev/null +++ b/bench/algorithm/mandelbrot/1-m.goose @@ -0,0 +1,223 @@ +// The rows of 1.goose, rendered by worker threads a few rows to a job. The +// main program puts the rows into the image as they come back and hashes +// the image as far as it is complete while the rest is rendered. + +import std; +import os; + +let MAX_ITER = 50; +let BLOCK = 32; + +// The bits of the pixels at cr in the row at ci, eight to a byte with the +// leftmost pixel in the high bit. +inline fn mbrot(cr: f64[BLOCK], ci: f64) -> u8[BLOCK / 8] { + var bytes: u8[BLOCK / 8] = [0; BLOCK / 8]; + var zr: f64[BLOCK] = [0.0; BLOCK]; + var zi: f64[BLOCK] = [0.0; BLOCK]; + var tr: f64[BLOCK] = [0.0; BLOCK]; + var ti: f64[BLOCK] = [0.0; BLOCK]; + for i in MAX_ITER / 5 { + for j in 5 { + for k in BLOCK { + zi[k] = (zr[k] + zr[k]) * zi[k] + ci; + zr[k] = tr[k] - ti[k] + cr[k]; + tr[k] = zr[k] * zr[k]; + ti[k] = zi[k] * zi[k]; + } + } + var inside = false; + for k in BLOCK { + if tr[k] + ti[k] <= 4.0 { inside = true; } + } + if !inside { return bytes; } + } + for b in BLOCK / 8 { + var bits: u8 = 0; + for k in 8 { + if tr[b * 8 + k] + ti[b * 8 + k] <= 4.0 { bits |= 0x80 >> k; } + } + bytes[b] = bits; + } + return bytes; +} + +// One row of the image, BLOCK / 8 bytes to a block of xs. +simd fn render_row(out: u8[:], xs: f64[BLOCK][:], ci: f64) { + assert(out.len == xs.len * (BLOCK / 8)); + for cr, i in xs { + copy_into(out[i * (BLOCK / 8)..(i + 1) * (BLOCK / 8)], mbrot(cr, ci)); + } +} + +// The real parts of a row's pixels, in blocks; the last block runs past size. +fn x_coords(size: i64) -> f64[BLOCK][>..] { + let inv = 2.0 / size; + var xs: f64[BLOCK][>..] = []; + for x in (size + BLOCK - 1) / BLOCK { + var cr: f64[BLOCK] = [0.0; BLOCK]; + for k in BLOCK { cr[k] = (x * BLOCK + k) * inv - 1.5; } + xs.push(cr); + } + return xs; +} + +// MD5 (RFC 1321). +const MD5_S: u32[64] = [ + 7, 12, 17, 22, 7, 12, 17, 22, 7, 12, 17, 22, 7, 12, 17, 22, + 5, 9, 14, 20, 5, 9, 14, 20, 5, 9, 14, 20, 5, 9, 14, 20, + 4, 11, 16, 23, 4, 11, 16, 23, 4, 11, 16, 23, 4, 11, 16, 23, + 6, 10, 15, 21, 6, 10, 15, 21, 6, 10, 15, 21, 6, 10, 15, 21, +]; +const MD5_K: u32[64] = [ + 0xd76aa478, 0xe8c7b756, 0x242070db, 0xc1bdceee, 0xf57c0faf, 0x4787c62a, 0xa8304613, 0xfd469501, + 0x698098d8, 0x8b44f7af, 0xffff5bb1, 0x895cd7be, 0x6b901122, 0xfd987193, 0xa679438e, 0x49b40821, + 0xf61e2562, 0xc040b340, 0x265e5a51, 0xe9b6c7aa, 0xd62f105d, 0x02441453, 0xd8a1e681, 0xe7d3fbc8, + 0x21e1cde6, 0xc33707d6, 0xf4d50d87, 0x455a14ed, 0xa9e3e905, 0xfcefa3f8, 0x676f02d9, 0x8d2a4c8a, + 0xfffa3942, 0x8771f681, 0x6d9d6122, 0xfde5380c, 0xa4beea44, 0x4bdecfa9, 0xf6bb4b60, 0xbebfbc70, + 0x289b7ec6, 0xeaa127fa, 0xd4ef3085, 0x04881d05, 0xd9d4d039, 0xe6db99e5, 0x1fa27cf8, 0xc4ac5665, + 0xf4292244, 0x432aff97, 0xab9423a7, 0xfc93a039, 0x655b59c3, 0x8f0ccc92, 0xffeff47d, 0x85845dd1, + 0x6fa87e4f, 0xfe2ce6e0, 0xa3014314, 0x4e0811a1, 0xf7537e82, 0xbd3af235, 0x2ad7d2bb, 0xeb86d391, +]; + +fn md5_block(h: u32[4]&, blk: u8[:]) { + assert(blk.len == 64); + var w: u32[16] = [0; 16]; + for i in 16 { + w[i] = blk[4 * i] | (blk[4 * i + 1] as u32 << 8) | (blk[4 * i + 2] as u32 << 16) | + (blk[4 * i + 3] as u32 << 24); + } + var a = h[0]; + var b = h[1]; + var c = h[2]; + var d = h[3]; + for i in 16 { + let t = a + ((b & c) | (~b & d)) + MD5_K[i] + w[i]; + a = d; + d = c; + c = b; + b = b + ((t << MD5_S[i]) | (t >> (32 - MD5_S[i]))); + } + for i in 16..32 { + let t = a + ((d & b) | (~d & c)) + MD5_K[i] + w[(5 * i + 1) % 16]; + a = d; + d = c; + c = b; + b = b + ((t << MD5_S[i]) | (t >> (32 - MD5_S[i]))); + } + for i in 32..48 { + let t = a + (b ^ c ^ d) + MD5_K[i] + w[(3 * i + 5) % 16]; + a = d; + d = c; + c = b; + b = b + ((t << MD5_S[i]) | (t >> (32 - MD5_S[i]))); + } + for i in 48..64 { + let t = a + (c ^ (b | ~d)) + MD5_K[i] + w[(7 * i) % 16]; + a = d; + d = c; + c = b; + b = b + ((t << MD5_S[i]) | (t >> (32 - MD5_S[i]))); + } + h[0] += a; + h[1] += b; + h[2] += c; + h[3] += d; +} + +struct Md5 { h: u32[4], len: i64 } + +fn md5() -> Md5 { Md5 { h: [0x67452301, 0xefcdab89, 0x98badcfe, 0x10325476], len: 0 } } + +// Hashes data, a whole number of blocks. +fn md5_blocks(m: Md5&, data: u8[:]) { + assert(data.len % 64 == 0); + var i = 0; + while i < data.len { + md5_block(m.h, data[i..i + 64]); + i += 64; + } + m.len += data.len; +} + +// Hashes rest, less than a block, and the padding, and appends the digest +// in hex to out. +fn md5_finish(m: Md5&, rest: u8[:], out: u8[>..]&) { + assert(rest.len < 64); + var tail: u8[128] = [0; 128]; + for k in rest.len { tail[k] = rest[k]; } + tail[rest.len] = 0x80; + let end = if rest.len < 56 { 64 } else { 128 }; + let bitlen = ((m.len + rest.len) as u64) << 3; + for k in 8 { tail[end - 8 + k] = (bitlen >> (8 * k)) as! u8; } + md5_block(m.h, tail[0..64]); + if end == 128 { md5_block(m.h, tail[64..128]); } + for word in m.h { + for k in 4 { format_uint(out, (word >> (8 * k)) & 0xff, 16, 2, '0'); } + } +} + +let ROWS_PER_JOB = 8; + +struct Job { y0: i64, y1: i64 } // y0 < 0: no more work +struct Rows { y0: i64, bits: u8[] } + +fn render_rows(xs: f64[BLOCK][:], size: i64, y0: i64, y1: i64) -> u8[>..] { + let inv = 2.0 / size; + let row_bytes = size / 8; + var out: u8[>..] = []; + out.resize((y1 - y0) * row_bytes, 0); + var row: u8[>..] = []; + row.resize(xs.len * (BLOCK / 8), 0); + for y in y0..y1 { + render_row(row, xs, y * inv - 1.0); + copy_into(out[(y - y0) * row_bytes..(y - y0 + 1) * row_bytes], row[..row_bytes]); + } + return out; +} + +thread_fn worker(size: i64) { + let xs = x_coords(size); + loop { + let job = qget(); + guard job.y0 >= 0 else { break; } + qput(Rows { y0: job.y0, bits: render_rows(xs, size, job.y0, job.y1) }); + } +} + +fn main() { + let a = args(); + var n = 200; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + let size = (n + 7) / 8 * 8; + let row_bytes = size / 8; + let jobs = (size + ROWS_PER_JOB - 1) / ROWS_PER_JOB; + let workers = max(1, hardware_threads()); + for j in jobs { qput(Job { y0: j * ROWS_PER_JOB, y1: min(size, (j + 1) * ROWS_PER_JOB) }); } + for w in workers { qput(Job { y0: -1, y1: -1 }); } + var ids: i64[>..] = []; + for w in workers { ids.push(thread_spawn(worker, size)); } + var image: u8[>..] = []; + image.resize(size * row_bytes, 0); + var ready: bool[>..] = []; + ready.resize(jobs, false); + var next = 0; + var hashed = 0; + var m = md5(); + for k in jobs { + let r = qget(); + let at = r.y0 * row_bytes; + copy_into(image[at..at + r.bits.len], r.bits); + ready[r.y0 / ROWS_PER_JOB] = true; + while next < jobs && ready[next] { next++; } + let upto = min(size, next * ROWS_PER_JOB) * row_bytes / 64 * 64; + if upto > hashed { + md5_blocks(m, image[hashed..upto]); + hashed = upto; + } + } + for id in ids { thread_wait(id); } + var out: u8[>..] = []; + format(out, "P4\n", size, " ", size, "\n"); + md5_finish(m, image[hashed..], out); + print(out); +} diff --git a/bench/algorithm/mandelbrot/1.goose b/bench/algorithm/mandelbrot/1.goose new file mode 100644 index 000000000..ec74dd1a4 --- /dev/null +++ b/bench/algorithm/mandelbrot/1.goose @@ -0,0 +1,185 @@ +// Pixels are iterated a block of BLOCK at a time, a lane each, so that the +// iteration is a few independent vector operations and the block stops as +// soon as every pixel in it has escaped, checked every five iterations as in +// the other programs. A pixel that has escaped stays escaped (its |z| only +// grows, to infinity and then NaN), so stopping a block early or late +// changes no bit of the image. The image is hashed with MD5, written here. + +import std; +import os; + +let MAX_ITER = 50; +let BLOCK = 32; + +// The bits of the pixels at cr in the row at ci, eight to a byte with the +// leftmost pixel in the high bit. +inline fn mbrot(cr: f64[BLOCK], ci: f64) -> u8[BLOCK / 8] { + var bytes: u8[BLOCK / 8] = [0; BLOCK / 8]; + var zr: f64[BLOCK] = [0.0; BLOCK]; + var zi: f64[BLOCK] = [0.0; BLOCK]; + var tr: f64[BLOCK] = [0.0; BLOCK]; + var ti: f64[BLOCK] = [0.0; BLOCK]; + for i in MAX_ITER / 5 { + for j in 5 { + for k in BLOCK { + zi[k] = (zr[k] + zr[k]) * zi[k] + ci; + zr[k] = tr[k] - ti[k] + cr[k]; + tr[k] = zr[k] * zr[k]; + ti[k] = zi[k] * zi[k]; + } + } + var inside = false; + for k in BLOCK { + if tr[k] + ti[k] <= 4.0 { inside = true; } + } + if !inside { return bytes; } + } + for b in BLOCK / 8 { + var bits: u8 = 0; + for k in 8 { + if tr[b * 8 + k] + ti[b * 8 + k] <= 4.0 { bits |= 0x80 >> k; } + } + bytes[b] = bits; + } + return bytes; +} + +// One row of the image, BLOCK / 8 bytes to a block of xs. +simd fn render_row(out: u8[:], xs: f64[BLOCK][:], ci: f64) { + assert(out.len == xs.len * (BLOCK / 8)); + for cr, i in xs { + copy_into(out[i * (BLOCK / 8)..(i + 1) * (BLOCK / 8)], mbrot(cr, ci)); + } +} + +// The real parts of a row's pixels, in blocks; the last block runs past size. +fn x_coords(size: i64) -> f64[BLOCK][>..] { + let inv = 2.0 / size; + var xs: f64[BLOCK][>..] = []; + for x in (size + BLOCK - 1) / BLOCK { + var cr: f64[BLOCK] = [0.0; BLOCK]; + for k in BLOCK { cr[k] = (x * BLOCK + k) * inv - 1.5; } + xs.push(cr); + } + return xs; +} + +// MD5 (RFC 1321). +const MD5_S: u32[64] = [ + 7, 12, 17, 22, 7, 12, 17, 22, 7, 12, 17, 22, 7, 12, 17, 22, + 5, 9, 14, 20, 5, 9, 14, 20, 5, 9, 14, 20, 5, 9, 14, 20, + 4, 11, 16, 23, 4, 11, 16, 23, 4, 11, 16, 23, 4, 11, 16, 23, + 6, 10, 15, 21, 6, 10, 15, 21, 6, 10, 15, 21, 6, 10, 15, 21, +]; +const MD5_K: u32[64] = [ + 0xd76aa478, 0xe8c7b756, 0x242070db, 0xc1bdceee, 0xf57c0faf, 0x4787c62a, 0xa8304613, 0xfd469501, + 0x698098d8, 0x8b44f7af, 0xffff5bb1, 0x895cd7be, 0x6b901122, 0xfd987193, 0xa679438e, 0x49b40821, + 0xf61e2562, 0xc040b340, 0x265e5a51, 0xe9b6c7aa, 0xd62f105d, 0x02441453, 0xd8a1e681, 0xe7d3fbc8, + 0x21e1cde6, 0xc33707d6, 0xf4d50d87, 0x455a14ed, 0xa9e3e905, 0xfcefa3f8, 0x676f02d9, 0x8d2a4c8a, + 0xfffa3942, 0x8771f681, 0x6d9d6122, 0xfde5380c, 0xa4beea44, 0x4bdecfa9, 0xf6bb4b60, 0xbebfbc70, + 0x289b7ec6, 0xeaa127fa, 0xd4ef3085, 0x04881d05, 0xd9d4d039, 0xe6db99e5, 0x1fa27cf8, 0xc4ac5665, + 0xf4292244, 0x432aff97, 0xab9423a7, 0xfc93a039, 0x655b59c3, 0x8f0ccc92, 0xffeff47d, 0x85845dd1, + 0x6fa87e4f, 0xfe2ce6e0, 0xa3014314, 0x4e0811a1, 0xf7537e82, 0xbd3af235, 0x2ad7d2bb, 0xeb86d391, +]; + +fn md5_block(h: u32[4]&, blk: u8[:]) { + assert(blk.len == 64); + var w: u32[16] = [0; 16]; + for i in 16 { + w[i] = blk[4 * i] | (blk[4 * i + 1] as u32 << 8) | (blk[4 * i + 2] as u32 << 16) | + (blk[4 * i + 3] as u32 << 24); + } + var a = h[0]; + var b = h[1]; + var c = h[2]; + var d = h[3]; + for i in 16 { + let t = a + ((b & c) | (~b & d)) + MD5_K[i] + w[i]; + a = d; + d = c; + c = b; + b = b + ((t << MD5_S[i]) | (t >> (32 - MD5_S[i]))); + } + for i in 16..32 { + let t = a + ((d & b) | (~d & c)) + MD5_K[i] + w[(5 * i + 1) % 16]; + a = d; + d = c; + c = b; + b = b + ((t << MD5_S[i]) | (t >> (32 - MD5_S[i]))); + } + for i in 32..48 { + let t = a + (b ^ c ^ d) + MD5_K[i] + w[(3 * i + 5) % 16]; + a = d; + d = c; + c = b; + b = b + ((t << MD5_S[i]) | (t >> (32 - MD5_S[i]))); + } + for i in 48..64 { + let t = a + (c ^ (b | ~d)) + MD5_K[i] + w[(7 * i) % 16]; + a = d; + d = c; + c = b; + b = b + ((t << MD5_S[i]) | (t >> (32 - MD5_S[i]))); + } + h[0] += a; + h[1] += b; + h[2] += c; + h[3] += d; +} + +struct Md5 { h: u32[4], len: i64 } + +fn md5() -> Md5 { Md5 { h: [0x67452301, 0xefcdab89, 0x98badcfe, 0x10325476], len: 0 } } + +// Hashes data, a whole number of blocks. +fn md5_blocks(m: Md5&, data: u8[:]) { + assert(data.len % 64 == 0); + var i = 0; + while i < data.len { + md5_block(m.h, data[i..i + 64]); + i += 64; + } + m.len += data.len; +} + +// Hashes rest, less than a block, and the padding, and appends the digest +// in hex to out. +fn md5_finish(m: Md5&, rest: u8[:], out: u8[>..]&) { + assert(rest.len < 64); + var tail: u8[128] = [0; 128]; + for k in rest.len { tail[k] = rest[k]; } + tail[rest.len] = 0x80; + let end = if rest.len < 56 { 64 } else { 128 }; + let bitlen = ((m.len + rest.len) as u64) << 3; + for k in 8 { tail[end - 8 + k] = (bitlen >> (8 * k)) as! u8; } + md5_block(m.h, tail[0..64]); + if end == 128 { md5_block(m.h, tail[64..128]); } + for word in m.h { + for k in 4 { format_uint(out, (word >> (8 * k)) & 0xff, 16, 2, '0'); } + } +} + +fn main() { + let a = args(); + var n = 200; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + let size = (n + 7) / 8 * 8; + let inv = 2.0 / size; + let xs = x_coords(size); + let row_bytes = size / 8; + var image: u8[>..] = []; + image.resize(size * row_bytes, 0); + var row: u8[>..] = []; + row.resize(xs.len * (BLOCK / 8), 0); + for y in size { + render_row(row, xs, y * inv - 1.0); + copy_into(image[y * row_bytes..(y + 1) * row_bytes], row[..row_bytes]); + } + var out: u8[>..] = []; + format(out, "P4\n", size, " ", size, "\n"); + var m = md5(); + let full = image.len / 64 * 64; + md5_blocks(m, image[..full]); + md5_finish(m, image[full..], out); + print(out); +} diff --git a/bench/algorithm/merkletrees/1.goose b/bench/algorithm/merkletrees/1.goose new file mode 100644 index 000000000..82a679cfd --- /dev/null +++ b/bench/algorithm/merkletrees/1.goose @@ -0,0 +1,98 @@ +// The problem's rules forbid arenas: every node is allocated on its own and +// freed on its own. The nodes live in one `reusable` pool, a grow-only array +// with a freelist, so `alloc_ref` takes the slot a `free` last handed back, +// the way a size-class allocator does, and the links are 4-byte offsets +// from the pool's base. A tree is freed node by node, children first, as +// Rust drops its boxes. + +import std; +import os; + +let MIN_DEPTH = 4; + +// hash and value are the reference's two Options, each a flag and a value. +struct Node { + hash: i64, + value: i64, + l: Node&?, + r: Node&?, + has_hash: bool, + has_value: bool, +} + +reusable var pool: Node[>..] = []; + +recursive fn make(depth: i64) -> Node& { + if depth == 0 { return pool.alloc_ref(Node { value: 1, has_value: true, .. }); } + let l: Node& = make(depth - 1); + let r: Node& = make(depth - 1); + return pool.alloc_ref(Node { l: l, r: r, .. }); +} + +recursive fn free_tree(n: Node&) { + let l = n.l; + if l { free_tree(l); } + let r = n.r; + if r { free_tree(r); } + pool.free(pool.index_of(n)); +} + +fn get_hash(n: Node&) -> i64 { if n.has_hash { n.hash } else { -1 } } + +recursive fn cal_hash(n: Node&) { + guard !n.has_hash; + if n.has_value { + n.hash = n.value; + n.has_hash = true; + return; + } + let l = n.l; + guard l; + let r = n.r; + guard r; + cal_hash(l); + cal_hash(r); + n.hash = get_hash(l) + get_hash(r); + n.has_hash = true; +} + +recursive fn check(n: Node&) -> bool { + guard n.has_hash else { return false; } + if n.has_value { return true; } + let l = n.l; + guard l else { return false; } + let r = n.r; + guard r else { return false; } + return check(l) && check(r); +} + +fn main() { + let a = args(); + var n = 10; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + let max_depth = max(MIN_DEPTH + 2, n); + { + let stretch: Node& = make(max_depth + 1); + cal_hash(stretch); + print("stretch tree of depth ", max_depth + 1, "\t root hash: ", get_hash(stretch), + " check: ", check(stretch)); + free_tree(stretch); + } + let long_lived: Node& = make(max_depth); + var d = MIN_DEPTH; + while d <= max_depth { + let iterations = 1 << (max_depth - d + MIN_DEPTH); + var sum = 0; + for i in iterations { + let t: Node& = make(d); + cal_hash(t); + sum += get_hash(t); + free_tree(t); + } + print(iterations, "\t trees of depth ", d, "\t root hash sum: ", sum); + d += 2; + } + cal_hash(long_lived); + print("long lived tree of depth ", max_depth, "\t root hash: ", get_hash(long_lived), + " check: ", check(long_lived)); +} diff --git a/bench/algorithm/nbody/1.goose b/bench/algorithm/nbody/1.goose new file mode 100644 index 000000000..084962ba2 --- /dev/null +++ b/bench/algorithm/nbody/1.goose @@ -0,0 +1,110 @@ +// The bodies are a fixed array the run keeps in a local, and the pair loops +// have constant bounds, so the C compiler unrolls them completely and keeps +// every body in registers, as the C++ program by Jeffrey Beu gets its +// templates to do. (2.goose lays the bodies out for vector registers.) + +import std; +import os; +import math; + +let SOLAR_MASS = 4.0 * PI * PI; +let DAYS_PER_YEAR = 365.24; +let DT = 0.01; + +struct Body { x: f64, y: f64, z: f64, vx: f64, vy: f64, vz: f64, mass: f64 } + +fn body(x: f64, y: f64, z: f64, vx: f64, vy: f64, vz: f64, mass: f64) -> Body { + Body { x, y, z, vx * DAYS_PER_YEAR, vy * DAYS_PER_YEAR, vz * DAYS_PER_YEAR, mass * SOLAR_MASS } +} + +fn solar_system() -> Body[5] { + [ + body(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0), + body(4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01, + 1.66007664274403694e-03, 7.69901118419740425e-03, -6.90460016972063023e-05, + 9.54791938424326609e-04), + body(8.34336671824457987e+00, 4.12479856412430479e+00, -4.03523417114321381e-01, + -2.76742510726862411e-03, 4.99852801234917238e-03, 2.30417297573763929e-05, + 2.85885980666130812e-04), + body(1.28943695621391310e+01, -1.51111514016986312e+01, -2.23307578892655734e-01, + 2.96460137564761618e-03, 2.37847173959480950e-03, -2.96589568540237556e-05, + 4.36624404335156298e-05), + body(1.53796971148509165e+01, -2.59193146099879641e+01, 1.79258772950371181e-01, + 2.68067772490389322e-03, 1.62824170038242295e-03, -9.51592254519715870e-05, + 5.15138902046611451e-05), + ] +} + +fn offset_momentum(b: Body[5]&) { + var px = 0.0; + var py = 0.0; + var pz = 0.0; + for q in b { + px += q.vx * q.mass; + py += q.vy * q.mass; + pz += q.vz * q.mass; + } + b[0].vx = -px / SOLAR_MASS; + b[0].vy = -py / SOLAR_MASS; + b[0].vz = -pz / SOLAR_MASS; +} + +fn energy(b: Body[5]&) -> f64 { + var e = 0.0; + for i in 5 { + e += 0.5 * b[i].mass * (b[i].vx * b[i].vx + b[i].vy * b[i].vy + b[i].vz * b[i].vz); + for j in i + 1..5 { + let dx = b[i].x - b[j].x; + let dy = b[i].y - b[j].y; + let dz = b[i].z - b[j].z; + e -= b[i].mass * b[j].mass / sqrt(dx * dx + dy * dy + dz * dz); + } + } + return e; +} + +fn advance(b: Body[5]&, steps: i64) { + var p = b; + for s in steps { + for i in 5 { + for j in i + 1..5 { + let dx = p[i].x - p[j].x; + let dy = p[i].y - p[j].y; + let dz = p[i].z - p[j].z; + let d2 = dx * dx + dy * dy + dz * dz; + let mag = DT / (d2 * sqrt(d2)); + let mi = p[i].mass * mag; + let mj = p[j].mass * mag; + p[i].vx -= dx * mj; + p[i].vy -= dy * mj; + p[i].vz -= dz * mj; + p[j].vx += dx * mi; + p[j].vy += dy * mi; + p[j].vz += dz * mi; + } + } + for &q in p { + q.x += DT * q.vx; + q.y += DT * q.vy; + q.z += DT * q.vz; + } + } + b = p; +} + +fn print_energy(b: Body[5]&) { + var line: u8[>..] = []; + format_flt(line, energy(b), 9); + print(line); +} + +fn main() { + let a = args(); + var n = 1000; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + var b = solar_system(); + offset_momentum(b); + print_energy(b); + advance(b, n); + print_energy(b); +} diff --git a/bench/algorithm/nbody/2.goose b/bench/algorithm/nbody/2.goose new file mode 100644 index 000000000..940a55022 --- /dev/null +++ b/bench/algorithm/nbody/2.goose @@ -0,0 +1,194 @@ +// The five bodies sit in the lanes of 8-wide arrays, one array per +// coordinate. Since five is odd, every pair is (i, i + 1) or (i, i + 2) mod 5 +// for exactly one i, so a step is two rounds of lane-wise arithmetic on +// whole arrays, each between the bodies and a rotated copy of them: one +// vector square root and one vector division a round, and nothing but +// rotations moving data between lanes. Lanes 5 to 7 are padding, computed +// along and never read into lanes 0 to 4. +// +// As in the C program by Christoph Bauer, the velocities are kept multiplied +// by the time step and the masses by its square while advancing, which takes +// a multiplication off every step. + +import std; +import os; +import math; + +let SOLAR_MASS = 4.0 * PI * PI; +let DAYS_PER_YEAR = 365.24; +let DT = 0.01; + +struct Bodies { + x: f64[8], y: f64[8], z: f64[8], + vx: f64[8], vy: f64[8], vz: f64[8], + mass: f64[8], +} + +fn solar_system() -> Bodies { + var b = Bodies { .. }; + let p: f64[7][5] = [ + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0], + [4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01, + 1.66007664274403694e-03, 7.69901118419740425e-03, -6.90460016972063023e-05, + 9.54791938424326609e-04], + [8.34336671824457987e+00, 4.12479856412430479e+00, -4.03523417114321381e-01, + -2.76742510726862411e-03, 4.99852801234917238e-03, 2.30417297573763929e-05, + 2.85885980666130812e-04], + [1.28943695621391310e+01, -1.51111514016986312e+01, -2.23307578892655734e-01, + 2.96460137564761618e-03, 2.37847173959480950e-03, -2.96589568540237556e-05, + 4.36624404335156298e-05], + [1.53796971148509165e+01, -2.59193146099879641e+01, 1.79258772950371181e-01, + 2.68067772490389322e-03, 1.62824170038242295e-03, -9.51592254519715870e-05, + 5.15138902046611451e-05], + ]; + for i in 5 { + b.x[i] = p[i][0]; + b.y[i] = p[i][1]; + b.z[i] = p[i][2]; + b.vx[i] = p[i][3] * DAYS_PER_YEAR; + b.vy[i] = p[i][4] * DAYS_PER_YEAR; + b.vz[i] = p[i][5] * DAYS_PER_YEAR; + b.mass[i] = p[i][6] * SOLAR_MASS; + } + return b; +} + +fn offset_momentum(b: Bodies&) { + var px = 0.0; + var py = 0.0; + var pz = 0.0; + for i in 5 { + px += b.vx[i] * b.mass[i]; + py += b.vy[i] * b.mass[i]; + pz += b.vz[i] * b.mass[i]; + } + b.vx[0] = -px / SOLAR_MASS; + b.vy[0] = -py / SOLAR_MASS; + b.vz[0] = -pz / SOLAR_MASS; +} + +fn energy(b: Bodies&) -> f64 { + var e = 0.0; + for i in 5 { + e += 0.5 * b.mass[i] * (b.vx[i] * b.vx[i] + b.vy[i] * b.vy[i] + b.vz[i] * b.vz[i]); + for j in i + 1..5 { + let dx = b.x[i] - b.x[j]; + let dy = b.y[i] - b.y[j]; + let dz = b.z[i] - b.z[j]; + e -= b.mass[i] * b.mass[j] / sqrt(dx * dx + dy * dy + dz * dz); + } + } + return e; +} + +// Lane l of up1(a) holds a[l + 1 mod 5], of down1(a) a[l - 1 mod 5], and so on. +inline fn up1(a: f64[8]) -> f64[8] { [a[1], a[2], a[3], a[4], a[0], a[5], a[6], a[7]] } +inline fn up2(a: f64[8]) -> f64[8] { [a[2], a[3], a[4], a[0], a[1], a[5], a[6], a[7]] } +inline fn down1(a: f64[8]) -> f64[8] { [a[4], a[0], a[1], a[2], a[3], a[5], a[6], a[7]] } +inline fn down2(a: f64[8]) -> f64[8] { [a[3], a[4], a[0], a[1], a[2], a[5], a[6], a[7]] } + +simd fn advance(b: Bodies&, steps: i64) { + var x = b.x; + var y = b.y; + var z = b.z; + var vx = b.vx; + var vy = b.vy; + var vz = b.vz; + var m = b.mass; + for l in 8 { + vx[l] *= DT; + vy[l] *= DT; + vz[l] *= DT; + m[l] *= DT * DT; + } + let m1 = up1(m); + let m2 = up2(m); + for s in steps { + let x1 = up1(x); + let y1 = up1(y); + let z1 = up1(z); + let x2 = up2(x); + let y2 = up2(y); + let z2 = up2(z); + // Pair (l, l + 1) gives body l + 1 g1 and takes f1 from body l; + // pair (l, l + 2) gives body l + 2 g2 and takes f2 from body l. + var f1x: f64[8] = [0.0; 8]; + var f1y: f64[8] = [0.0; 8]; + var f1z: f64[8] = [0.0; 8]; + var g1x: f64[8] = [0.0; 8]; + var g1y: f64[8] = [0.0; 8]; + var g1z: f64[8] = [0.0; 8]; + var f2x: f64[8] = [0.0; 8]; + var f2y: f64[8] = [0.0; 8]; + var f2z: f64[8] = [0.0; 8]; + var g2x: f64[8] = [0.0; 8]; + var g2y: f64[8] = [0.0; 8]; + var g2z: f64[8] = [0.0; 8]; + for l in 8 { + let dx = x[l] - x1[l]; + let dy = y[l] - y1[l]; + let dz = z[l] - z1[l]; + let d2 = dx * dx + dy * dy + dz * dz; + let mag = 1.0 / (d2 * sqrt(d2)); + let ex = x[l] - x2[l]; + let ey = y[l] - y2[l]; + let ez = z[l] - z2[l]; + let e2 = ex * ex + ey * ey + ez * ez; + let mag2 = 1.0 / (e2 * sqrt(e2)); + f1x[l] = dx * m1[l] * mag; + f1y[l] = dy * m1[l] * mag; + f1z[l] = dz * m1[l] * mag; + g1x[l] = dx * m[l] * mag; + g1y[l] = dy * m[l] * mag; + g1z[l] = dz * m[l] * mag; + f2x[l] = ex * m2[l] * mag2; + f2y[l] = ey * m2[l] * mag2; + f2z[l] = ez * m2[l] * mag2; + g2x[l] = ex * m[l] * mag2; + g2y[l] = ey * m[l] * mag2; + g2z[l] = ez * m[l] * mag2; + } + let h1x = down1(g1x); + let h1y = down1(g1y); + let h1z = down1(g1z); + let h2x = down2(g2x); + let h2y = down2(g2y); + let h2z = down2(g2z); + for l in 8 { + vx[l] += (h1x[l] - f1x[l]) + (h2x[l] - f2x[l]); + vy[l] += (h1y[l] - f1y[l]) + (h2y[l] - f2y[l]); + vz[l] += (h1z[l] - f1z[l]) + (h2z[l] - f2z[l]); + x[l] += vx[l]; + y[l] += vy[l]; + z[l] += vz[l]; + } + } + for l in 8 { + vx[l] /= DT; + vy[l] /= DT; + vz[l] /= DT; + } + b.x = x; + b.y = y; + b.z = z; + b.vx = vx; + b.vy = vy; + b.vz = vz; +} + +fn print_energy(b: Bodies&) { + var line: u8[>..] = []; + format_flt(line, energy(b), 9); + print(line); +} + +fn main() { + let a = args(); + var n = 1000; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + var b = solar_system(); + offset_momentum(b); + print_energy(b); + advance(b, n); + print_energy(b); +} diff --git a/bench/algorithm/nsieve/1.goose b/bench/algorithm/nsieve/1.goose new file mode 100644 index 000000000..45d455b7a --- /dev/null +++ b/bench/algorithm/nsieve/1.goose @@ -0,0 +1,30 @@ +import std; +import os; + +// The flags are one byte each, as in the C, C++ and first Rust programs. + +fn nsieve(n: i64) { + var flags: bool[>..] = []; + flags.resize(n, true); + var count = 0; + for i in 2..n { + if flags[i] { + count++; + var j = i << 1; + while j < n { flags[j] = false; j += i; } + } + } + var line: u8[>..] = []; + line.append("Primes up to "); + format_int(line, n, 10, 8, ' '); + line.push(' '); + format_int(line, count, 10, 8, ' '); + print(line); +} + +fn main() { + let a = args(); + var n = 4; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + for i in 3 { nsieve(10000 << (n - i)); } +} diff --git a/bench/algorithm/nsieve/2.goose b/bench/algorithm/nsieve/2.goose new file mode 100644 index 000000000..3e20473a4 --- /dev/null +++ b/bench/algorithm/nsieve/2.goose @@ -0,0 +1,34 @@ +import std; +import os; + +// The flags are bits of u32 words, a set bit marking a composite, as in the +// second Rust program. + +fn nsieve(n: i64) { + var flags: u32[>..] = []; + flags.resize((n + 31) >> 5, 0); + var count = 0; + for i in 2..n { + if (flags[i >> 5] >> i) & 1 == 0 { + count++; + var j = i << 1; + while j < n { + flags[j >> 5] |= (1 as u32) << j; + j += i; + } + } + } + var line: u8[>..] = []; + line.append("Primes up to "); + format_int(line, n, 10, 8, ' '); + line.push(' '); + format_int(line, count, 10, 8, ' '); + print(line); +} + +fn main() { + let a = args(); + var n = 4; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + for i in 3 { nsieve(10000 << (n - i)); } +} diff --git a/bench/algorithm/pidigits/1.goose b/bench/algorithm/pidigits/1.goose new file mode 100644 index 000000000..34bb29543 --- /dev/null +++ b/bench/algorithm/pidigits/1.goose @@ -0,0 +1,214 @@ +// The spigot of the Rust programs: n1 / d and n2 / d bound the digits still +// to come, a term of the series narrows them, and a digit is emitted once +// both quotients agree. The integers are 64-bit limbs, least significant +// first, with the shortcuts of GMP's functions: the products by small +// factors and the sums of a step in one pass over the limbs, and a quotient +// known to be small taken from the leading limbs, with a pass over all of +// them only where those cannot decide it. + +import std; +import os; + +// The high half of a 64x64-bit product. clang compiles the pattern to one +// mulx, but only while it sees it on its own: inlined first, it rewrites the +// halves of an operand it knows something about (2k + 1, 10) and the pattern +// is gone. `recursive` is what keeps a Goose function out of line. +recursive fn mulhi(a: u64, b: u64) -> u64 { + let al = a & 0xffffffff; + let ah = a >> 32; + let bl = b & 0xffffffff; + let bh = b >> 32; + let ll = al * bl; + let lh = al * bh; + let hl = ah * bl; + let hh = ah * bh; + let mid = (ll >> 32) + (lh & 0xffffffff) + (hl & 0xffffffff); + return hh + (lh >> 32) + (hl >> 32) + (mid >> 32); +} + +fn bits(a: u64[:]) -> i64 { + var i = a.len - 1; + while i >= 0 { + if a[i] != 0 { return i * 64 + 64 - clz(a[i]); } + i--; + } + return 0; +} + +// The 64 bits of a from bit s up. +fn bits_at(a: u64[:], s: i64) -> u64 { + let w = s >> 6; + let b = s & 63; + var x = a[w] >> b; + if b > 0 && w + 1 < a.len { x |= a[w + 1] << (64 - b); } + return x; +} + +// Whether q * d <= n. +fn fits(n: u64[:], d: u64[:], q: u64) -> bool { + var c: u64 = 0; + var borrow = false; + for i in max(n.len, d.len + 1) { + var lo: u64 = 0; + var hi: u64 = 0; + if i < d.len { lo = d[i] * q; hi = mulhi(d[i], q); } + lo += c; + if lo < c { hi++; } + c = hi; + let x: u64 = if i < n.len { n[i] } else { 0 }; + let y = x - lo; + let b1 = x < lo; + borrow = b1 || (borrow && y == 0); + } + return !borrow; +} + +// n / d, for a quotient below 2^32. The leading 64 bits of n, and the bits +// of d at the same place, bound it to [lo, hi]; only where those differ +// does it take a pass over the limbs. +fn quotient(n: u64[:], d: u64[:]) -> u64 { + let nb = bits(n); + if nb < bits(d) { return 0; } + if nb <= 64 { return bits_at(n, 0) / bits_at(d, 0); } + let s = nb - 64; + let nt = bits_at(n, s); + let dt = bits_at(d, s); + assert(dt >> 32 != 0); + let lo = if dt == 0xffffffffffffffff { 0 } else { nt / (dt + 1) }; + var hi = nt / dt; + if nt % dt == dt - 1 { hi++; } + var q = hi; + while q > lo && !fits(n, d, q) { q--; } + return q; +} + +// n1, n2, d = n1 (2k - 1) + 2 n2, n1 (k - 1) + n2 (k + 2), d (2k + 1), all of +// the same length, in one pass over their limbs. +fn next_term(n1: u64[>..]&, n2: u64[>..]&, d: u64[>..]&, k: u64) { + let len = n1.len; + assert(n2.len == len && d.len == len); + let a1 = 2 * k - 1; + let a2 = k - 1; + let b2 = k + 2; + let m = 2 * k + 1; + var c1: u64 = 0; + var c2: u64 = 0; + var c3: u64 = 0; + for i in len { + let x = n1[i]; + let y = n2[i]; + var lo1 = x * a1; + var hi1 = mulhi(x, a1); + let y2 = y << 1; + lo1 += y2; + if lo1 < y2 { hi1++; } + hi1 += y >> 63; + lo1 += c1; + if lo1 < c1 { hi1++; } + var lo2 = x * a2; + var hi2 = mulhi(x, a2); + let ylo = y * b2; + let yhi = mulhi(y, b2); + lo2 += ylo; + if lo2 < ylo { hi2++; } + hi2 += yhi; + lo2 += c2; + if lo2 < c2 { hi2++; } + let z = d[i]; + var lo3 = z * m; + var hi3 = mulhi(z, m); + lo3 += c3; + if lo3 < c3 { hi3++; } + n1[i] = lo1; + n2[i] = lo2; + d[i] = lo3; + c1 = hi1; + c2 = hi2; + c3 = hi3; + } + if c1 != 0 || c2 != 0 || c3 != 0 { + n1.push(c1); + n2.push(c2); + d.push(c3); + } +} + +// n1, n2 = (n1 - q d) 10, (n2 - q d) 10, for q d at most either, in one pass. +fn eliminate_digit(n1: u64[>..]&, n2: u64[>..]&, d: u64[>..]&, q: u64) { + let len = n1.len; + assert(n2.len == len && d.len == len); + var c: u64 = 0; + var borrow1 = false; + var borrow2 = false; + var c1: u64 = 0; + var c2: u64 = 0; + for i in len { + var lo = d[i] * q; + var hi = mulhi(d[i], q); + lo += c; + if lo < c { hi++; } + c = hi; + let x1 = n1[i]; + let y1 = x1 - lo; + let z1 = if borrow1 { y1 - 1 } else { y1 }; + borrow1 = x1 < lo || (borrow1 && y1 == 0); + var p1 = z1 * 10; + let h1 = mulhi(z1, 10); + p1 += c1; + c1 = if p1 < c1 { h1 + 1 } else { h1 }; + n1[i] = p1; + let x2 = n2[i]; + let y2 = x2 - lo; + let z2 = if borrow2 { y2 - 1 } else { y2 }; + borrow2 = x2 < lo || (borrow2 && y2 == 0); + var p2 = z2 * 10; + let h2 = mulhi(z2, 10); + p2 += c2; + c2 = if p2 < c2 { h2 + 1 } else { h2 }; + n2[i] = p2; + } + if c1 != 0 || c2 != 0 { + n1.push(c1); + n2.push(c2); + d.push(0); + } + while n1.len > 1 && n1[n1.len - 1] == 0 && n2[n2.len - 1] == 0 && d[d.len - 1] == 0 { + n1.pop(); + n2.pop(); + d.pop(); + } +} + +fn main() { + let a = args(); + var digits_to_print = 27; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { digits_to_print = v; } } + var out: u8[>..] = []; + var digits_printed = 0; + var k: u64 = 1; + var n1: u64[>..] = [4]; + var n2: u64[>..] = [3]; + var d: u64[>..] = [1]; + loop { + let u = quotient(n1, d); + let v = quotient(n2, d); + if u == v { + out.push(('0' + u) as! u8); + digits_printed++; + let m = digits_printed % 10; + if m == 0 { format(out, "\t:", digits_printed, "\n"); } + if digits_printed >= digits_to_print { + if m > 0 { + for i in 10 - m { out.push(' '); } + format(out, "\t:", digits_printed, "\n"); + } + break; + } + eliminate_digit(n1, n2, d, u); + } else { + next_term(n1, n2, d, k); + k++; + } + } + write_stdout(out); +} diff --git a/bench/algorithm/regex-redux/1.goose b/bench/algorithm/regex-redux/1.goose new file mode 100644 index 000000000..bb9fb9281 --- /dev/null +++ b/bench/algorithm/regex-redux/1.goose @@ -0,0 +1,61 @@ +import std; +import os; +import regex; + +let VARIANTS = [ + "agggtaaa|tttaccct", + "[cgt]gggtaaa|tttaccc[acg]", + "a[act]ggtaaa|tttacc[agt]t", + "ag[act]gtaaa|tttac[agt]ct", + "agg[act]taaa|ttta[agt]cct", + "aggg[acg]aaa|ttt[cgt]ccct", + "agggt[cgt]aa|tt[acg]accct", + "agggta[cgt]a|t[acg]taccct", + "agggtaa[cgt]|[acg]ttaccct", +]; + +let SUBST_PATTERNS = ["tHa[Nt]", "aND|caN|Ha[DS]|WaS", "a[NSt]|BY", "<[^>]*>", "\\|[^|][^|]*\\|"]; +let SUBST_WITH = ["<4>", "<3>", "<2>", "|", "-"]; + +// text with every match of re replaced by `with`, appended to out. +fn replace_all(out: u8[>..]&, re: regex::Regex&, text: u8[:], with: u8[:]) { + var last = 0; + regex::each_match(re, text) { s, e => + out.append(text[last..s]); + out.append(with); + last = e; + }; + out.append(text[last..]); +} + +fn main() { + let argv = args(); + var input: u8[>..] = []; + if !read_file(if argv.len > 1 { argv[1] } else { "25000_in" }, input) { abort("cannot read the input"); } + + var seq: u8[>..] = []; + let strip = regex::compile(">[^\n]*\n|\n"); + replace_all(seq, strip, input, ""); + + for v in VARIANTS { + let re = regex::compile(v); + print(v, " ", regex::count(re, seq)); + } + + // The substitutions go back and forth between two buffers, ending in a. + var a: u8[>..] = []; + var b: u8[>..] = []; + for p, i in SUBST_PATTERNS { + let re = regex::compile(p); + if i == 0 { + replace_all(a, re, seq, SUBST_WITH[i]); + } else if i % 2 == 1 { + b.clear(); + replace_all(b, re, a, SUBST_WITH[i]); + } else { + a.clear(); + replace_all(a, re, b, SUBST_WITH[i]); + } + } + print("\n", input.len, "\n", seq.len, "\n", a.len); +} diff --git a/bench/algorithm/regex-redux/2-m.goose b/bench/algorithm/regex-redux/2-m.goose new file mode 100644 index 000000000..bd77f6737 --- /dev/null +++ b/bench/algorithm/regex-redux/2-m.goose @@ -0,0 +1,90 @@ +// 1.goose with the variant counts on worker threads while the main thread +// does the substitutions, which take about as long as three of the counts. +// Threads share no memory, so each worker is handed its own copy of the +// sequence when it is spawned, and takes variants off a queue. + +import std; +import os; +import regex; + +let VARIANTS = [ + "agggtaaa|tttaccct", + "[cgt]gggtaaa|tttaccc[acg]", + "a[act]ggtaaa|tttacc[agt]t", + "ag[act]gtaaa|tttac[agt]ct", + "agg[act]taaa|ttta[agt]cct", + "aggg[acg]aaa|ttt[cgt]ccct", + "agggt[cgt]aa|tt[acg]accct", + "agggta[cgt]a|t[acg]taccct", + "agggtaa[cgt]|[acg]ttaccct", +]; + +let SUBST_PATTERNS = ["tHa[Nt]", "aND|caN|Ha[DS]|WaS", "a[NSt]|BY", "<[^>]*>", "\\|[^|][^|]*\\|"]; +let SUBST_WITH = ["<4>", "<3>", "<2>", "|", "-"]; + +let MAX_WORKERS = 3; + +struct Job { k: i32 } // a variant to count; -1: no more +struct Found { k: i32, n: i64 } + +// text with every match of re replaced by `with`, appended to out. +fn replace_all(out: u8[>..]&, re: regex::Regex&, text: u8[:], with: u8[:]) { + var last = 0; + regex::each_match(re, text) { s, e => + out.append(text[last..s]); + out.append(with); + last = e; + }; + out.append(text[last..]); +} + +thread_fn counter(seq: u8[]) { + loop { + let job = qget(); + guard job.k >= 0 else { break; } + let re = regex::compile(VARIANTS[job.k]); + qput(Found { k: job.k, n: regex::count(re, seq) }); + } +} + +fn main() { + let argv = args(); + var input: u8[>..] = []; + if !read_file(if argv.len > 1 { argv[1] } else { "25000_in" }, input) { abort("cannot read the input"); } + + var seq: u8[>..] = []; + let strip = regex::compile(">[^\n]*\n|\n"); + replace_all(seq, strip, input, ""); + + let workers = clamp(hardware_threads() - 1, 1, MAX_WORKERS); + for k in VARIANTS.len { qput(Job { k: k as i32 }); } + for w in workers { qput(Job { k: -1 }); } + var ids: i64[>..] = []; + for w in workers { ids.push(thread_spawn(counter, copy(seq))); } + + // The substitutions go back and forth between two buffers, ending in a. + var a: u8[>..] = []; + var b: u8[>..] = []; + for p, i in SUBST_PATTERNS { + let re = regex::compile(p); + if i == 0 { + replace_all(a, re, seq, SUBST_WITH[i]); + } else if i % 2 == 1 { + b.clear(); + replace_all(b, re, a, SUBST_WITH[i]); + } else { + a.clear(); + replace_all(a, re, b, SUBST_WITH[i]); + } + } + + var counts: i64[>..] = []; + counts.resize(VARIANTS.len, 0); + for k in VARIANTS.len { + let f = qget(); + counts[f.k] = f.n; + } + for id in ids { thread_wait(id); } + for v, k in VARIANTS { print(v, " ", counts[k]); } + print("\n", input.len, "\n", seq.len, "\n", a.len); +} diff --git a/bench/algorithm/secp256k1/1.goose b/bench/algorithm/secp256k1/1.goose new file mode 100644 index 000000000..5de6e8947 --- /dev/null +++ b/bench/algorithm/secp256k1/1.goose @@ -0,0 +1,540 @@ +// secp256k1 point multiplication with the algorithm of the hand-written +// programs (1.rs, after 1.ts): Jacobian coordinates, the scalar split by the +// GLV endomorphism into two halves of about 128 bits that share one chain of +// doublings, and an extended-Euclid inversion back to affine coordinates. +// Field elements are four 64-bit limbs, kept below p = 2^256 - 2^32 - 977; +// the scalar split and the inversion use 384-bit two's complement integers. +// The constants are the Rust program's, written out as limbs. + +import std; +import os; + +type Fe = u64[4]; +type Big = u64[6]; + +struct Jp { x: Fe, y: Fe, z: Fe } + +let PC: u64 = 0x1000003D1; // 2^256 - p + +let P: Fe = [0xFFFFFFFEFFFFFC2F, 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF]; +let N: Big = [0xBFD25E8CD0364141, 0xBAAEDCE6AF48A03B, 0xFFFFFFFFFFFFFFFE, 0xFFFFFFFFFFFFFFFF, 0, 0]; +let BETA: Fe = [0xc1396c28719501ee, 0x9cf0497512f58995, 0x6e64479eac3434e9, 0x7ae96a2b657c0710]; +let A1: Big = [0xe86c90e49284eb15, 0x3086d221a7d46bcd, 0, 0, 0, 0]; +let MINUS_B1: Big = [0x6f547fa90abfe4c3, 0xe4437ed6010e8828, 0, 0, 0, 0]; +let A2: Big = [0x57c1108d9d44cfd8, 0x14ca50f7a8e2f3f6, 1, 0, 0, 0]; +let PRIVATE_KEY: Big = [0xD9132D70B5BC7693, 0xB64592E3FE0AB0DA, 0x4FCA3EF970FF4D38, 0x2DEE927079283C3C, 0, 0]; +let GX: Fe = [0x59F2815B16F81798, 0x029BFCDB2DCE28D9, 0x55A06295CE870B07, 0x79BE667EF9DCBBAC]; +let GY: Fe = [0x9C47D08FFB10D4B8, 0xFD17B448A6855419, 0x5DA4FBFC0E1108A8, 0x483ADA7726A3C465]; + +// The high half of a 64x64-bit product. clang compiles the pattern to one +// mulx, but only while it sees it on its own: inlined first, it rewrites the +// halves of an operand it knows something about (a constant, a square) and +// the pattern is gone. `recursive` is what keeps a Goose function out of line. +recursive fn mulhi(a: u64, b: u64) -> u64 { + let al = a & 0xffffffff; + let ah = a >> 32; + let bl = b & 0xffffffff; + let bh = b >> 32; + let ll = al * bl; + let lh = al * bh; + let hl = ah * bl; + let hh = ah * bh; + let mid = (ll >> 32) + (lh & 0xffffffff) + (hl & 0xffffffff); + return hh + (lh >> 32) + (hl >> 32) + (mid >> 32); +} + +// --- the field --------------------------------------------------------------- + +// t (eight limbs) modulo p, using 2^256 = PC (mod p). +inline fn reduce(t: u64[8]) -> Fe { + var r: Fe = [0; 4]; + var c: u64 = 0; + for i in 4 { + let h = t[i + 4]; + var lo = h * PC; + var hi = mulhi(h, PC); + lo += c; + if lo < c { hi++; } + let x = t[i] + lo; + if x < lo { hi++; } + r[i] = x; + c = hi; + } + // c < 2^34, and c * PC < 2^68 + var lo = c * PC; + var hi = mulhi(c, PC); + r[0] += lo; + if r[0] < lo { hi++; } + r[1] += hi; + var carry = r[1] < hi; + for i in 2..4 { + if carry { r[i]++; carry = r[i] == 0; } + } + if carry { r[0] += PC; } // what wraps past 2^256 is tiny, so this cannot carry + return canon(r); +} + +// r modulo p, for r < 2^256. +inline fn canon(r: Fe) -> Fe { + var s: Fe = r; + s[0] += PC; + var carry = s[0] < PC; + for i in 1..4 { + if carry { s[i]++; carry = s[i] == 0; } + } + if carry { return s; } + return r; +} + +inline fn fe_mul(a: Fe, b: Fe) -> Fe { + var t: u64[8] = [0; 8]; + for i in 4 { + var c: u64 = 0; + let ai = a[i]; + for j in 4 { + var lo = ai * b[j]; + var hi = mulhi(ai, b[j]); + lo += c; + if lo < c { hi++; } + let x = t[i + j] + lo; + if x < lo { hi++; } + t[i + j] = x; + c = hi; + } + t[i + 4] = c; + } + return reduce(t); +} + +// The products of distinct limbs once, doubled, and the squares of the limbs. +inline fn fe_sqr(a: Fe) -> Fe { + var t: u64[8] = [0; 8]; + for i in 3 { + var c: u64 = 0; + let ai = a[i]; + for j in i + 1..4 { + var lo = ai * a[j]; + var hi = mulhi(ai, a[j]); + lo += c; + if lo < c { hi++; } + let x = t[i + j] + lo; + if x < lo { hi++; } + t[i + j] = x; + c = hi; + } + t[i + 4] = c; + } + var top: u64 = 0; + for k in 1..8 { + let x = t[k]; + t[k] = (x << 1) | top; + top = x >> 63; + } + var c: u64 = 0; + for i in 4 { + var lo = a[i] * a[i]; + var hi = mulhi(a[i], a[i]); + lo += c; + if lo < c { hi++; } + let x = t[2 * i] + lo; + if x < lo { hi++; } + t[2 * i] = x; + let y = t[2 * i + 1] + hi; + c = if y < hi { 1 } else { 0 }; + t[2 * i + 1] = y; + } + return reduce(t); +} + +inline fn fe_add(a: Fe, b: Fe) -> Fe { + var s: Fe = [0; 4]; + var c = false; + for i in 4 { + let x = a[i] + b[i]; + let c1 = x < a[i]; + let y = if c { x + 1 } else { x }; + s[i] = y; + c = c1 || (c && y == 0); + } + if c { // a + b - 2^256 + PC is below p + s[0] += PC; + var k = s[0] < PC; + for i in 1..4 { if k { s[i]++; k = s[i] == 0; } } + return s; + } + return canon(s); +} + +inline fn fe_sub(a: Fe, b: Fe) -> Fe { + var s: Fe = [0; 4]; + var borrow = false; + for i in 4 { + let x = a[i] - b[i]; + let b1 = a[i] < b[i]; + let y = if borrow { x - 1 } else { x }; + s[i] = y; + borrow = b1 || (borrow && x == 0); + } + if borrow { // a - b + 2^256 - PC = a - b + p + let old = s[0]; + s[0] -= PC; + var k = old < PC; + for i in 1..4 { if k { k = s[i] == 0; s[i]--; } } + } + return s; +} + +fn fe_is_zero(a: Fe) -> bool { (a[0] | a[1] | a[2] | a[3]) == 0 } + +fn fe_neg(a: Fe) -> Fe { fe_sub([0; 4], a) } + +// --- 384-bit integers --------------------------------------------------------- + +fn big(a: Fe) -> Big { [a[0], a[1], a[2], a[3], 0, 0] } + +fn fe(a: Big) -> Fe { [a[0], a[1], a[2], a[3]] } + +fn big_add(a: Big, b: Big) -> Big { + var s: Big = [0; 6]; + var c: u64 = 0; + for i in 6 { + let x = a[i] + b[i]; + let y = x + c; + c = if x < a[i] || y < x { 1 } else { 0 }; + s[i] = y; + } + return s; +} + +fn big_neg(a: Big) -> Big { + var s: Big = [0; 6]; + var c: u64 = 1; + for i in 6 { + let y = ~a[i] + c; + c = if y < c { 1 } else { 0 }; + s[i] = y; + } + return s; +} + +fn big_sub(a: Big, b: Big) -> Big { big_add(a, big_neg(b)) } + +// The product modulo 2^384, which is the product itself for every operand +// pair the program forms. +fn big_mul(a: Big, b: Big) -> Big { + var t: Big = [0; 6]; + for i in 6 { + var c: u64 = 0; + let ai = a[i]; + if ai == 0 { continue; } + for j in 6 - i { + var lo = ai * b[j]; + var hi = mulhi(ai, b[j]); + lo += c; + if lo < c { hi++; } + let x = t[i + j] + lo; + if x < lo { hi++; } + t[i + j] = x; + c = hi; + } + } + return t; +} + +fn big_is_neg(a: Big) -> bool { (a[5] >> 63) != 0 } + +fn big_is_zero(a: Big) -> bool { (a[0] | a[1] | a[2] | a[3] | a[4] | a[5]) == 0 } + +// Unsigned comparison: -1, 0 or 1. +fn big_cmp(a: Big, b: Big) -> i64 { + var i = 5; + while i >= 0 { + if a[i] != b[i] { return if a[i] < b[i] { -1 } else { 1 }; } + i--; + } + return 0; +} + +fn big_bits(a: Big) -> i64 { + var i = 5; + while i >= 0 { + if a[i] != 0 { return i * 64 + 64 - clz(a[i]); } + i--; + } + return 0; +} + +fn big_shl(a: Big, k: i64) -> Big { + var s: Big = [0; 6]; + let w = k >> 6; + let b = k & 63; + for i in w..6 { + var x = a[i - w] << b; + if b > 0 && i - w > 0 { x |= a[i - w - 1] >> (64 - b); } + s[i] = x; + } + return s; +} + +fn big_shr1(a: Big) -> Big { + var s: Big = [0; 6]; + for i in 6 { + var x = a[i] >> 1; + if i < 5 { x |= a[i + 1] << 63; } + s[i] = x; + } + return s; +} + +// Quotient and remainder of non-negative a and positive d, by shifting and +// subtracting. +fn big_divmod(a: Big, d: Big) -> Big, Big { + var q: Big = [0; 6]; + var r = a; + let shift = big_bits(a) - big_bits(d); + if shift < 0 { return q, r; } + var m = big_shl(d, shift); + var i = shift; + while i >= 0 { + if big_cmp(r, m) >= 0 { + r = big_sub(r, m); + q[i >> 6] |= (1 as u64) << (i & 63); + } + m = big_shr1(m); + i--; + } + return q, r; +} + +// The Euclidean remainder of a signed a modulo a positive m. +fn big_rem(a: Big, m: Big) -> Big { + if !big_is_neg(a) { + let q, r = big_divmod(a, m); + return r; + } + let q, r = big_divmod(big_neg(a), m); + if big_is_zero(r) { return r; } + return big_sub(m, r); +} + +fn fe_bits(a: Fe) -> i64 { + var i = 3; + while i >= 0 { + if a[i] != 0 { return i * 64 + 64 - clz(a[i]); } + i--; + } + return 0; +} + +// The 64 bits of a from bit s up. +fn fe_bits_at(a: Fe, s: i64) -> u64 { + let w = s >> 6; + let b = s & 63; + var x = a[w] >> b; + if b > 0 && w < 3 { x |= a[w + 1] << (64 - b); } + return x; +} + +fn big_mul_u64(a: Big, q: u64) -> Big { + var t: Big = [0; 6]; + var c: u64 = 0; + for i in 6 { + var lo = a[i] * q; + var hi = mulhi(a[i], q); + lo += c; + if lo < c { hi++; } + t[i] = lo; + c = hi; + } + return t; +} + +// b / a and b % a where the quotient is below 2^32, from the leading 64 bits +// of both and a correction of at most a few steps; false where the quotient +// may be larger. +fn divmod_short(b: Fe, a: Fe) -> u64, Fe, bool { + let nb = fe_bits(b); + if nb - fe_bits(a) > 31 { return 0, b, false; } + let s = max(nb - 64, 0); + var q = fe_bits_at(b, s) / fe_bits_at(a, s); + // r = b - q * a, five limbs wide + var r: u64[5] = [b[0], b[1], b[2], b[3], 0]; + var c: u64 = 0; + var borrow = false; + for i in 5 { + var lo: u64 = 0; + var hi: u64 = 0; + if i < 4 { lo = a[i] * q; hi = mulhi(a[i], q); } + lo += c; + if lo < c { hi++; } + let x = r[i] - lo; + let b1 = r[i] < lo; + let y = if borrow { x - 1 } else { x }; + borrow = b1 || (borrow && x == 0); + r[i] = y; + c = hi; + } + while (r[4] >> 63) != 0 { // q was too large + q--; + var k: u64 = 0; + for i in 5 { + let ai = if i < 4 { a[i] } else { 0 }; + let x = r[i] + ai; + let y = x + k; + k = if x < ai || y < x { 1 } else { 0 }; + r[i] = y; + } + } + var rem: Fe = [r[0], r[1], r[2], r[3]]; + while fe_bits(rem) >= fe_bits(a) && big_cmp(big(rem), big(a)) >= 0 { + q++; + rem = fe(big_sub(big(rem), big(a))); + } + return q, rem, true; +} + +fn invert(number: Fe) -> Fe { + var a = number; + var b = P; + var x: Big = [0; 6]; + var y: Big = [1, 0, 0, 0, 0, 0]; + var u: Big = [1, 0, 0, 0, 0, 0]; + var v: Big = [0; 6]; + while !fe_is_zero(a) { + var m: Big = [0; 6]; + var n: Big = [0; 6]; + let qs, rs, short = divmod_short(b, a); + var r = rs; + if short { + m = big_sub(x, big_mul_u64(u, qs)); + n = big_sub(y, big_mul_u64(v, qs)); + } else { + let q, rb = big_divmod(big(b), big(a)); + r = fe(rb); + m = big_sub(x, big_mul(u, q)); + n = big_sub(y, big_mul(v, q)); + } + b = a; + a = r; + x = u; + y = v; + u = m; + v = n; + } + return fe(big_rem(x, big(P))); +} + +fn div_nearest(a: Big, b: Big) -> Big { + let q, r = big_divmod(big_add(a, big_shr1(b)), b); + return q; +} + +// k = k1 + k2 * lambda (mod n), with k1 and k2 of about 128 bits, each with +// the sign it is to be applied with. +fn split_scalar_endo(k: Big) -> bool, Big, bool, Big { + let b2 = A1; + let c1 = div_nearest(big_mul(b2, k), N); + let c2 = div_nearest(big_mul(MINUS_B1, k), N); + var k1 = big_rem(big_sub(big_sub(k, big_mul(c1, A1)), big_mul(c2, A2)), N); + var k2 = big_rem(big_sub(big_mul(c1, MINUS_B1), big_mul(c2, b2)), N); + let pow_2_128: Big = [0, 0, 1, 0, 0, 0]; + let k1neg = big_cmp(k1, pow_2_128) > 0; + let k2neg = big_cmp(k2, pow_2_128) > 0; + if k1neg { k1 = big_sub(N, k1); } + if k2neg { k2 = big_sub(N, k2); } + return k1neg, k1, k2neg, k2; +} + +// --- the curve ----------------------------------------------------------------- + +fn jp_zero() -> Jp { Jp { x: [0; 4], y: [1, 0, 0, 0], z: [0; 4] } } + +fn jp_double(p: Jp) -> Jp { + let a = fe_sqr(p.x); + let b = fe_sqr(p.y); + let c = fe_sqr(b); + let t = fe_sub(fe_sub(fe_sqr(fe_add(p.x, b)), a), c); + let d = fe_add(t, t); + let e = fe_add(fe_add(a, a), a); + let f = fe_sqr(e); + let x3 = fe_sub(f, fe_add(d, d)); + let c2 = fe_add(c, c); + let c4 = fe_add(c2, c2); + let y3 = fe_sub(fe_mul(e, fe_sub(d, x3)), fe_add(c4, c4)); + let yz = fe_mul(p.y, p.z); + return Jp { x: x3, y: y3, z: fe_add(yz, yz) }; +} + +fn jp_add(p: Jp, q: Jp) -> Jp { + if fe_is_zero(q.x) || fe_is_zero(q.y) { return p; } + if fe_is_zero(p.x) || fe_is_zero(p.y) { return q; } + let z1z1 = fe_sqr(p.z); + let z2z2 = fe_sqr(q.z); + let u1 = fe_mul(p.x, z2z2); + let u2 = fe_mul(q.x, z1z1); + let s1 = fe_mul(fe_mul(p.y, q.z), z2z2); + let s2 = fe_mul(fe_mul(q.y, p.z), z1z1); + let h = fe_sub(u2, u1); + let r = fe_sub(s2, s1); + if fe_is_zero(h) { + if fe_is_zero(r) { return jp_double(p); } + return jp_zero(); + } + let hh = fe_sqr(h); + let hhh = fe_mul(h, hh); + let v = fe_mul(u1, hh); + let x3 = fe_sub(fe_sub(fe_sqr(r), hhh), fe_add(v, v)); + let y3 = fe_sub(fe_mul(r, fe_sub(v, x3)), fe_mul(s1, hhh)); + let z3 = fe_mul(fe_mul(p.z, q.z), h); + return Jp { x: x3, y: y3, z: z3 }; +} + +fn bit(k: Big, i: i64) -> bool { ((k[i >> 6] >> (i & 63)) & 1) != 0 } + +fn jp_mul(p: Jp, scalar: Big) -> Jp { + let k1neg, k1, k2neg, k2 = split_scalar_endo(scalar); + var k1p = jp_zero(); + var k2p = jp_zero(); + var d = p; + let bits = max(big_bits(k1), big_bits(k2)); + for i in bits { + if bit(k1, i) { k1p = jp_add(k1p, d); } + if bit(k2, i) { k2p = jp_add(k2p, d); } + d = jp_double(d); + } + if k1neg { k1p.y = fe_neg(k1p.y); } + if k2neg { k2p.y = fe_neg(k2p.y); } + k2p.x = fe_mul(k2p.x, BETA); + return jp_add(k1p, k2p); +} + +fn point_mul(x: Fe, y: Fe, scalar: Big) -> Fe, Fe { + let r = jp_mul(Jp { x: x, y: y, z: [1, 0, 0, 0] }, scalar); + let inv_z = invert(r.z); + let inv_z_pow = fe_sqr(inv_z); + return fe_mul(r.x, inv_z_pow), fe_mul(fe_mul(r.y, inv_z), inv_z_pow); +} + +fn format_hex(out: u8[>..]&, a: Fe) { + var i = 3; + while i > 0 && a[i] == 0 { i--; } + format_uint(out, a[i], 16, 0, '0'); + while i > 0 { + i--; + format_uint(out, a[i], 16, 16, '0'); + } +} + +fn main() { + let a = args(); + var n = 1; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + var x = GX; + var y = GY; + for i in n { x, y = point_mul(x, y, PRIVATE_KEY); } + var line: u8[>..] = []; + format_hex(line, x); + line.push(','); + format_hex(line, y); + print(line); +} diff --git a/bench/algorithm/spectral-norm/1-m.goose b/bench/algorithm/spectral-norm/1-m.goose new file mode 100644 index 000000000..3f093ed83 --- /dev/null +++ b/bench/algorithm/spectral-norm/1-m.goose @@ -0,0 +1,110 @@ +// The products of 1.goose, with each product's rows dealt out to worker +// threads in as many slices as there are workers. A worker shares no memory +// with the main program, so a slice's job carries the vector it multiplies, +// and its result comes back as the slice's part of the product. + +import std; +import os; +import math; + +// dst[k] = row lo + k of A src, or with transposed, of A' src. +simd fn times(dst: f64[:], src: f64[:], lo: i64, transposed: bool) { + let n = src.len; + var i = 0; + while i < dst.len { + var a: f64[8] = [0.0; 8]; + var step: f64[8] = [0.0; 8]; + var sum: f64[8] = [0.0; 8]; + for k in 8 { + let r = (lo + i + k) as f64; + // A(r, 0) = r (r + 1) / 2 + r + 1, and A(0, r) = r (r + 1) / 2 + 1. + a[k] = r * (r + 1.0) * 0.5 + 1.0; + step[k] = r; + if transposed { step[k] += 1.0; } else { a[k] += r; } + sum[k] = src[0] / a[k]; + } + for j in 1..n { + let s = src[j]; + for k in 8 { + step[k] += 1.0; + a[k] += step[k]; + sum[k] += s / a[k]; + } + } + for k in 8 { + if i + k < dst.len { dst[i + k] = sum[k]; } + } + i += 8; + } +} + +struct Job { lo: i64, hi: i64, transposed: bool, src: f64[] } // hi < 0: no more work +struct Part { lo: i64, vals: f64[] } + +fn rows(job: Job&) -> f64[>..] { + var out: f64[>..] = []; + out.resize(job.hi - job.lo, 0.0); + times(out, job.src, job.lo, job.transposed); + return out; +} + +thread_fn worker() { + loop { + let job = qget(); + guard job.hi >= 0 else { break; } + qput(Part { lo: job.lo, vals: rows(job) }); + } +} + +fn product(dst: f64[:], src: f64[:], transposed: bool, slice: i64) { + let n = src.len; + var parts = 0; + var lo = 0; + while lo < n { + let hi = min(lo + slice, n); + qput(Job { lo: lo, hi: hi, transposed: transposed, src: src }); + parts++; + lo = hi; + } + for k in parts { + let p = qget(); + copy_into(dst[p.lo..p.lo + p.vals.len], p.vals); + } +} + +fn times_ata(dst: f64[:], src: f64[:], tmp: f64[:], slice: i64) { + product(tmp, src, false, slice); + product(dst, tmp, true, slice); +} + +fn main() { + let a = args(); + var n = 100; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + let workers = max(1, hardware_threads()); + // A multiple of the eight rows times works on at once. + let slice = ((n + workers - 1) / workers + 7) / 8 * 8; + var ids: i64[>..] = []; + for w in workers { ids.push(thread_spawn(worker)); } + var u: f64[>..] = []; + var v: f64[>..] = []; + var tmp: f64[>..] = []; + u.resize(n, 1.0); + v.resize(n, 0.0); + tmp.resize(n, 0.0); + for i in 10 { + times_ata(v, u, tmp, slice); + times_ata(u, v, tmp, slice); + } + for w in workers { qput(Job { lo: 0, hi: -1, transposed: false, src: [] }); } + for id in ids { thread_wait(id); } + var vbv = 0.0; + var vv = 0.0; + for i in n { + vbv += u[i] * v[i]; + vv += v[i] * v[i]; + } + var line: u8[>..] = []; + format_flt(line, sqrt(vbv / vv), 9); + print(line); +} diff --git a/bench/algorithm/spectral-norm/1.goose b/bench/algorithm/spectral-norm/1.goose new file mode 100644 index 000000000..b414187a9 --- /dev/null +++ b/bench/algorithm/spectral-norm/1.goose @@ -0,0 +1,72 @@ +// A product works on eight rows at once, one row a lane, so that each column +// is one vector division. Within a row the entries of A are made by addition +// as the column advances, as in the C++ program by Bela Pecsek: going from +// column j - 1 to j adds r + j to A(r, j) and r + j + 1 to A(j, r), all of +// them integers that doubles hold exactly. Each row's sum runs over the +// columns in order, so every product is the same to the bit as a plain +// row-by-row loop computes it. + +import std; +import os; +import math; + +// dst[k] = row lo + k of A src, or with transposed, of A' src. +simd fn times(dst: f64[:], src: f64[:], lo: i64, transposed: bool) { + let n = src.len; + var i = 0; + while i < dst.len { + var a: f64[8] = [0.0; 8]; + var step: f64[8] = [0.0; 8]; + var sum: f64[8] = [0.0; 8]; + for k in 8 { + let r = (lo + i + k) as f64; + // A(r, 0) = r (r + 1) / 2 + r + 1, and A(0, r) = r (r + 1) / 2 + 1. + a[k] = r * (r + 1.0) * 0.5 + 1.0; + step[k] = r; + if transposed { step[k] += 1.0; } else { a[k] += r; } + sum[k] = src[0] / a[k]; + } + for j in 1..n { + let s = src[j]; + for k in 8 { + step[k] += 1.0; + a[k] += step[k]; + sum[k] += s / a[k]; + } + } + for k in 8 { + if i + k < dst.len { dst[i + k] = sum[k]; } + } + i += 8; + } +} + +fn times_ata(dst: f64[:], src: f64[:], tmp: f64[:]) { + times(tmp, src, 0, false); + times(dst, tmp, 0, true); +} + +fn main() { + let a = args(); + var n = 100; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + var u: f64[>..] = []; + var v: f64[>..] = []; + var tmp: f64[>..] = []; + u.resize(n, 1.0); + v.resize(n, 0.0); + tmp.resize(n, 0.0); + for i in 10 { + times_ata(v, u, tmp); + times_ata(u, v, tmp); + } + var vbv = 0.0; + var vv = 0.0; + for i in n { + vbv += u[i] * v[i]; + vv += v[i] * v[i]; + } + var line: u8[>..] = []; + format_flt(line, sqrt(vbv / vv), 9); + print(line); +} diff --git a/bench/algorithm/spectral-norm/2-m.goose b/bench/algorithm/spectral-norm/2-m.goose new file mode 100644 index 000000000..b81ca413f --- /dev/null +++ b/bench/algorithm/spectral-norm/2-m.goose @@ -0,0 +1,145 @@ +// The products of 2.goose, with each product's rows dealt out to worker +// threads in as many slices as there are workers. A worker shares no memory +// with the main program, so a slice's job carries the vector it multiplies, +// and its result comes back as the slice's part of the product. + +import std; +import os; +import math; + +// dst[k] = row lo + k of A src, or with transposed, of A' src. +simd fn times(dst: f64[:], src: f64[:], lo: i64, transposed: bool) { + let cols = src[1..]; + var i = 0; + while i < dst.len { + // a holds the entry of the last column done, and d0 to d3 what the + // next four columns' entries add to it. + var a: f64[8] = [0.0; 8]; + var d0: f64[8] = [0.0; 8]; + var d1: f64[8] = [0.0; 8]; + var d2: f64[8] = [0.0; 8]; + var d3: f64[8] = [0.0; 8]; + var sum: f64[8] = [0.0; 8]; + for k in 8 { + let r = (lo + i + k) as f64; + // A(r, 0) = r (r + 1) / 2 + r + 1, and A(0, r) = r (r + 1) / 2 + 1. + a[k] = r * (r + 1.0) * 0.5 + 1.0; + var step = r + 1.0; + if transposed { step += 1.0; } else { a[k] += r; } + sum[k] = src[0] / a[k]; + d0[k] = step; + d1[k] = 2.0 * step + 1.0; + d2[k] = 3.0 * step + 3.0; + d3[k] = 4.0 * step + 6.0; + } + var j = 0; + while j + 4 <= cols.len { + let s0 = cols[j]; + let s1 = cols[j + 1]; + let s2 = cols[j + 2]; + let s3 = cols[j + 3]; + for k in 8 { + let a0 = a[k] + d0[k]; + let a1 = a[k] + d1[k]; + let a2 = a[k] + d2[k]; + let a3 = a[k] + d3[k]; + a[k] = a3; + d0[k] += 4.0; + d1[k] += 8.0; + d2[k] += 12.0; + d3[k] += 16.0; + // s0 / a0 + s1 / a1 + s2 / a2 + s3 / a3, over a0 a1 a2 a3. + let p01 = a0 * a1; + let p012 = p01 * a2; + let n01 = s0 * a1 + s1 * a0; + let n012 = n01 * a2 + s2 * p01; + sum[k] += (n012 * a3 + s3 * p012) / (p012 * a3); + } + j += 4; + } + while j < cols.len { + let s = cols[j]; + for k in 8 { + a[k] += d0[k]; + d0[k] += 1.0; + sum[k] += s / a[k]; + } + j++; + } + for k in 8 { + if i + k < dst.len { dst[i + k] = sum[k]; } + } + i += 8; + } +} + +struct Job { lo: i64, hi: i64, transposed: bool, src: f64[] } // hi < 0: no more work +struct Part { lo: i64, vals: f64[] } + +fn rows(job: Job&) -> f64[>..] { + var out: f64[>..] = []; + out.resize(job.hi - job.lo, 0.0); + times(out, job.src, job.lo, job.transposed); + return out; +} + +thread_fn worker() { + loop { + let job = qget(); + guard job.hi >= 0 else { break; } + qput(Part { lo: job.lo, vals: rows(job) }); + } +} + +fn product(dst: f64[:], src: f64[:], transposed: bool, slice: i64) { + let n = src.len; + var parts = 0; + var lo = 0; + while lo < n { + let hi = min(lo + slice, n); + qput(Job { lo: lo, hi: hi, transposed: transposed, src: src }); + parts++; + lo = hi; + } + for k in parts { + let p = qget(); + copy_into(dst[p.lo..p.lo + p.vals.len], p.vals); + } +} + +fn times_ata(dst: f64[:], src: f64[:], tmp: f64[:], slice: i64) { + product(tmp, src, false, slice); + product(dst, tmp, true, slice); +} + +fn main() { + let a = args(); + var n = 100; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + let workers = max(1, hardware_threads()); + // A multiple of the eight rows times works on at once. + let slice = ((n + workers - 1) / workers + 7) / 8 * 8; + var ids: i64[>..] = []; + for w in workers { ids.push(thread_spawn(worker)); } + var u: f64[>..] = []; + var v: f64[>..] = []; + var tmp: f64[>..] = []; + u.resize(n, 1.0); + v.resize(n, 0.0); + tmp.resize(n, 0.0); + for i in 10 { + times_ata(v, u, tmp, slice); + times_ata(u, v, tmp, slice); + } + for w in workers { qput(Job { lo: 0, hi: -1, transposed: false, src: [] }); } + for id in ids { thread_wait(id); } + var vbv = 0.0; + var vv = 0.0; + for i in n { + vbv += u[i] * v[i]; + vv += v[i] * v[i]; + } + var line: u8[>..] = []; + format_flt(line, sqrt(vbv / vv), 9); + print(line); +} diff --git a/bench/algorithm/spectral-norm/2.goose b/bench/algorithm/spectral-norm/2.goose new file mode 100644 index 000000000..83dd74d18 --- /dev/null +++ b/bench/algorithm/spectral-norm/2.goose @@ -0,0 +1,108 @@ +// A product works on eight rows at once, one row a lane, and on four columns +// at once: as in the C program with intrinsics by Miles, the four terms +// u[j] / A(r, j) of a row are summed over the product of their four +// denominators, so that the four columns take one vector division. Within +// a row the entries of A are made by addition as the column advances, as in +// the C++ program by Bela Pecsek: going from column j - 1 to j adds r + j to +// A(r, j) and r + j + 1 to A(j, r), all of them integers that doubles hold +// exactly. + +import std; +import os; +import math; + +// dst[k] = row lo + k of A src, or with transposed, of A' src. +simd fn times(dst: f64[:], src: f64[:], lo: i64, transposed: bool) { + let cols = src[1..]; + var i = 0; + while i < dst.len { + // a holds the entry of the last column done, and d0 to d3 what the + // next four columns' entries add to it. + var a: f64[8] = [0.0; 8]; + var d0: f64[8] = [0.0; 8]; + var d1: f64[8] = [0.0; 8]; + var d2: f64[8] = [0.0; 8]; + var d3: f64[8] = [0.0; 8]; + var sum: f64[8] = [0.0; 8]; + for k in 8 { + let r = (lo + i + k) as f64; + // A(r, 0) = r (r + 1) / 2 + r + 1, and A(0, r) = r (r + 1) / 2 + 1. + a[k] = r * (r + 1.0) * 0.5 + 1.0; + var step = r + 1.0; + if transposed { step += 1.0; } else { a[k] += r; } + sum[k] = src[0] / a[k]; + d0[k] = step; + d1[k] = 2.0 * step + 1.0; + d2[k] = 3.0 * step + 3.0; + d3[k] = 4.0 * step + 6.0; + } + var j = 0; + while j + 4 <= cols.len { + let s0 = cols[j]; + let s1 = cols[j + 1]; + let s2 = cols[j + 2]; + let s3 = cols[j + 3]; + for k in 8 { + let a0 = a[k] + d0[k]; + let a1 = a[k] + d1[k]; + let a2 = a[k] + d2[k]; + let a3 = a[k] + d3[k]; + a[k] = a3; + d0[k] += 4.0; + d1[k] += 8.0; + d2[k] += 12.0; + d3[k] += 16.0; + // s0 / a0 + s1 / a1 + s2 / a2 + s3 / a3, over a0 a1 a2 a3. + let p01 = a0 * a1; + let p012 = p01 * a2; + let n01 = s0 * a1 + s1 * a0; + let n012 = n01 * a2 + s2 * p01; + sum[k] += (n012 * a3 + s3 * p012) / (p012 * a3); + } + j += 4; + } + while j < cols.len { + let s = cols[j]; + for k in 8 { + a[k] += d0[k]; + d0[k] += 1.0; + sum[k] += s / a[k]; + } + j++; + } + for k in 8 { + if i + k < dst.len { dst[i + k] = sum[k]; } + } + i += 8; + } +} + +fn times_ata(dst: f64[:], src: f64[:], tmp: f64[:]) { + times(tmp, src, 0, false); + times(dst, tmp, 0, true); +} + +fn main() { + let a = args(); + var n = 100; + if a.len > 1 { let v, ok = parse_int(a[1]); if ok { n = v; } } + var u: f64[>..] = []; + var v: f64[>..] = []; + var tmp: f64[>..] = []; + u.resize(n, 1.0); + v.resize(n, 0.0); + tmp.resize(n, 0.0); + for i in 10 { + times_ata(v, u, tmp); + times_ata(u, v, tmp); + } + var vbv = 0.0; + var vv = 0.0; + for i in n { + vbv += u[i] * v[i]; + vv += v[i] * v[i]; + } + var line: u8[>..] = []; + format_flt(line, sqrt(vbv / vv), 9); + print(line); +} diff --git a/bench/bench_goose.yaml b/bench/bench_goose.yaml new file mode 100644 index 000000000..ce3f7ab26 --- /dev/null +++ b/bench/bench_goose.yaml @@ -0,0 +1,82 @@ +lang: goose +problems: + - name: helloworld + source: + - 1.goose + - name: nsieve + source: + - 1.goose + - 2.goose + - name: binarytrees + source: + - 1.goose + - 2-m.goose + - name: merkletrees + source: + - 1.goose + - name: nbody + source: + - 1.goose + - 2.goose + - name: spectral-norm + source: + - 1-m.goose + - 1.goose + - 2-m.goose + - 2.goose + - name: pidigits + source: + - 1.goose + - name: edigits + source: + - 1.goose + - name: mandelbrot + source: + - 1-m.goose + - 1.goose + - name: fannkuch-redux + source: + - 1-m.goose + - 1.goose + - name: fasta + source: + - 1-m.goose + - 1.goose + - name: knucleotide + source: + - 1-m.goose + - 1.goose + - 2.goose + - name: regex-redux + source: + - 1.goose + - 2-m.goose + - name: json-serde + source: + - 1.goose + - name: coro-prime-sieve + source: + - 1.goose + - name: lru + source: + - 1.goose + - 2.goose + - name: secp256k1 + source: + - 1.goose +source_rename_to: app.goose +compiler_version_command: goose version +compiler_version_regex: +runtime_version_parameter: +runtime_version_regex: +environments: + - os: linux + compiler: goose + version: latest + docker: + include: goose + build: goose -O2 --standalone -o app.c app.goose && clang @app.rsp + after_build: + - mv app out + out_dir: out + run_cmd: app diff --git a/bench/include/goose/.gitignore b/bench/include/goose/.gitignore new file mode 100644 index 000000000..24bb004ba --- /dev/null +++ b/bench/include/goose/.gitignore @@ -0,0 +1,2 @@ +!.gitignore +!app.rsp diff --git a/bench/include/goose/app.rsp b/bench/include/goose/app.rsp new file mode 100644 index 000000000..e4e5767ab --- /dev/null +++ b/bench/include/goose/app.rsp @@ -0,0 +1 @@ +-pipe -w -O3 -fomit-frame-pointer -march=broadwell -DGS_SIMD=1 -pthread -o app app.c -lm