Standard library · Numbers

num/col

Imported as import "num/col" as col, its names are then col.…. Every signature below is the one the checker infers.

Column aggregates, vectorized by hand.

A plain for over a buffer is vectorized by the compiler on its own, and for anything elementwise — a projection, a filter, a predicate count — that is the end of it: the generated code is what C's would be. Reductions over Float are the exception, and the reason is worth stating, because it is not a shortcoming of the loop.

Floating-point addition is not associative: (a + b) + c and a + (b + c) round differently, so they are different numbers. Vectorizing a sum means splitting it into independent partial sums and adding those at the end, which is exactly that reassociation — so a compiler may only do it if it has been told the program will accept a different answer. Rill does not tell it that. + on floats is loosened for contraction and nothing else, and the whole language keeps that promise, so LLVM has to keep the additions in the order they were written: it loads a vector's worth, extracts each lane, and adds them one at a time into a single accumulator. Every add waits on the one before it, and a loop that should stream at memory speed runs at the latency of fadd.

What that costs, on an M4, over a range that fits in L1: about six times. The answer is not to loosen the language but to write the accumulators out — several of them, so the adds are independent, in Float2 so each one holds two lanes. That is the same trade lib/num/grid.rill makes for the matrix product, and it lands in the same place: the kernels below run at 100 GB/s and better, matching a C loop compiled with -ffast-math, while a program that wants the strictly ordered sum still has it by writing the loop.

The order these add in is the order eight partial sums happen to reach the end, not the order the elements are in, so a sum here and a sum written as a loop may differ in the last bits. That is the whole of what is being traded, and it is traded per call rather than per build.

Every function takes a half-open range, lo up to but not hi, because a column arrives in chunks and a scan wants one chunk at a time without slicing a copy out of it. #v is the whole of it.

Reductions over Int are not here: integer addition is associative, so the compiler vectorizes those loops itself and a plain for is already as fast as this. The same goes for counting a predicate, which compiles to a vector compare and a subtract with no help from anyone.

prices = buf_float(4)
prices[0] = 1.5
prices[1] = 2.5
prices[2] = 3.0
prices[3] = 5.0
col.sum(prices, 0, #prices)          # => 12
col.mean(prices, 0, #prices)         # => Some(3)
col.mean(prices, 0, 0)               # => None

Functions

fn sum(v: Buf(Float), lo: Int, hi: Int) -> Float

Eight lanes at a time in four independent Float2 accumulators, then the elements past the last block of eight, added into the result the same way a scalar loop would.

v = buf_float(3)
v[0] = 0.5
v[1] = 0.25
v[2] = 0.25
col.sum(v, 0, 3)   # => 1
col.sum(v, 1, 3)   # => 0.5

fn sum_sq(v: Buf(Float), lo: Int, hi: Int) -> Float

The sum of squares, which is what a variance wants: one fma2 per lane pair, so the multiply and the add are one instruction and round once.

v = buf_float(2)
v[0] = 3.0
v[1] = 4.0
col.sum_sq(v, 0, 2)   # => 25

fn dot(x: Buf(Float), y: Buf(Float), lo: Int, hi: Int) -> Float

The dot product of two columns over the same range — a weighted sum, and what a covariance is made of.

x = buf_float(2)
y = buf_float(2)
x[0] = 1.0
x[1] = 2.0
y[0] = 3.0
y[1] = 4.0
col.dot(x, y, 0, 2)   # => 11

fn mean(v: Buf(Float), lo: Int, hi: Int) -> Option(Float)

None over an empty range, because the mean of nothing is not a number and an Option says so where a 0.0 would lie.

v = buf_float(2)
v[0] = 1.0
v[1] = 2.0
col.mean(v, 0, 2)   # => Some(1.5)
col.mean(v, 2, 2)   # => None

fn variance(v: Buf(Float), lo: Int, hi: Int) -> Option(Float)

The population variance, in one pass: E[x²] - E[x]². One pass is what a scan wants, and the price is the usual one — a column whose values are large and whose spread is small loses precision in the subtraction. None over an empty range.

v = buf_float(4)
v[0] = 2.0
v[1] = 4.0
v[2] = 4.0
v[3] = 6.0
col.variance(v, 0, 4)   # => Some(2)

fn min_of(v: Buf(Float), lo: Int, hi: Int) -> Option(Float)

v = buf_float(3)
v[0] = 2.0
v[1] = -1.0
v[2] = 5.0
col.min_of(v, 0, 3)   # => Some(-1)

fn max_of(v: Buf(Float), lo: Int, hi: Int) -> Option(Float)

v = buf_float(3)
v[0] = 2.0
v[1] = -1.0
v[2] = 5.0
col.max_of(v, 0, 3)   # => Some(5)
col.max_of(v, 0, 0)   # => None