Core Maths & Linear Algebra

The tables below are generated from the interpreter's own builtin dispatch table (the fn call_builtin match statement in engine/crates/qu-interp/src/lib.rs), not written speculatively ahead of the implementation — every function listed here runs today. Functions covering regression, clustering, and probability *models* (polyfit, ridge, kmeans, pca, normpdf, mvnpdf, ...) are documented in Statistics & Machine Learning instead, to avoid splitting that workflow across two chapters.

Example data

Names used by the examples below, defined here so each one runs.

x = linspace(0, 1, 12)
a = randn(3, 4, seed = 6)
b = randn(3, 4, seed = 7)
v = [1, 2, 3, 4, 5, 6]
n = 4
r = 3
c = 4
df = table(price = [10.0, 12.5, 9.75, 14.0], volume = [100, 80, 140, 60])

Elementwise Math

All of these accept a scalar, a Vec, or a Mat/Signal and apply elementwise; a handful (abs, angle/arg/phase, real/re, imag/ im, conj) are complex-aware and branch on whether the input is a real or complex value/array.

Overloads: abs

Case: scalar

x a real number. Returns a Num, the standard real absolute value.

print(abs(-5))   # 5

5

Case: vector

x a Vec. Returns a Vec of the same length, abs applied elementwise.

print(abs([-2, -0.5, 0.5, 2]))   # [2, 0.5, 0.5, 2]

[2, 0.5, 0.5, 2]

Case: complex

x a Complex value (the i/j-suffix literal). Returns a real Num — the magnitude sqrt(re^2 + im^2), not either component alone — the case that most often surprises newcomers.

print(abs(3 + 4j))   # 5 — magnitude, not the real part

5

Overloads: angle

Case: real

x a signed real number (scalar, Vec, or Mat). Degrades to 0 for a non-negative value or π for a negative one — there is no real angle to compute, so this is a convention, not a measurement.

print("angle(5)={angle(5)}  angle(-5)={angle(-5)}")

angle(5)=0  angle(-5)=3.141593

Case: complex

x a Complex value. Returns the true angle from the positive real axis, atan2(im, re), in radians.

z = 3 + 4j
print(angle(z))   # 0.927295

0.927295

Overloads: re

Case: real

x a real scalar, Vec, or Mat. Identity — passed through unchanged.

print(re(5))   # 5

5

Case: complex

x a Complex value. Returns just the real component, discarding the imaginary part.

print(re(3 + 4j))   # 3

3

Overloads: im

Case: real

x a real scalar, Vec, or Mat. Always 0, since a real value has no imaginary part.

print(im(5))   # 0

0

Case: complex

x a Complex value. Returns just the imaginary component.

print(im(3 + 4j))   # 4

4

Overloads: conj

Case: real

x a real scalar, Vec, or Mat. Identity — passed through unchanged, since conjugation is a no-op on a value with no imaginary part.

print(conj(5))   # 5

5

Case: complex

x a Complex value. Negates the imaginary part, leaving the real part alone.

print(conj(3 + 4j))   # 3 - 4j

3 - 4j

sign's exact-zero case

sign has its own exact-zero edge case worth seeing separately: it returns +1, never 0, per Rust's f64::signum.

print(sign([-2, -0.5, 0.5, 2]))   # [-1, -1, 1, 1] — sign(0) would be +1, never 0

[-1, -1, 1, 1]

The remaining functions in this table — the rest of the trig and inverse-trig family, the hyperbolic functions, the log family, the complex-component accessors (arg/phase/re/im/conj), and ceil — all follow the exact same elementwise pattern already shown above: one function in, one value of the same shape out, with only the domain notes already given in the table above to keep in mind. The example below runs each of them once on a scalar (or, for the complex ones, a single complex literal) so their actual return values are visible side by side.

x = 0.5
print("tan={tan(x)}  asin={asin(x)}  acos={acos(x)}  atan={atan(x)}")
print("sinh={sinh(x)}  cosh={cosh(x)}  tanh={tanh(x)}")
print("ln(e^2)={ln(exp(2))}  log(e^2)={log(exp(2))}  log2(8)={log2(8)}")

z = 3 + 4j
print("arg={arg(z)}  phase={phase(z)}  re={re(z)}  imag={imag(z)}  im={im(z)}  conj={conj(z)}")

print(ceil(2.1))   # 3

tan=0.546302  asin=0.523599  acos=1.047198  atan=0.463648
sinh=0.521095  cosh=1.127626  tanh=0.462117
ln(e^2)=2  log(e^2)=2  log2(8)=3
arg=0.927295  phase=0.927295  re=3  imag=4  im=4  conj=3 - 4j
3

Reductions & Descriptive Statistics

Basic order-statistic and spread reductions. corr, cov, mode, skewness, and kurtosis are documented in Statistics & Machine Learning. Several of the pure reductions below (sum, prod, mean, std, var, max, min, argmax, argmin, median, quantile) additionally accept an axis= keyword when x is a Mat: axis=0 collapses each column to one number (a length-cols Vec), axis=1 collapses each row (a length-rows Vec), instead of flattening the whole matrix to a single scalar.

mean, min, max, sum, std, var, median, quantile, argmin, and argmax also accept on_invalid= — R's na.rm= for a NaN or Infinity in the input, extended to cover Inf too: on_invalid="propagate" (the default, unchanged from before this keyword existed) lets a NaN/Inf flow through untouched; on_invalid="ignore" drops every NaN/Inf element before reducing (for mean/std/var this also shrinks the n the average/variance divides by, since it is computed over the already- filtered values; for argmin/argmax the returned index still refers to the position in the ORIGINAL, unfiltered x); on_invalid="error" raises instead of silently returning a NaN/Inf result, naming what was found and at what index. On a Mat under axis=, "ignore" drops NaN/Inf independently within each row/column, not across the whole matrix.

Overloads: sum

Case: real

x a real scalar, Vec, or Mat/Signal. Returns a real scalar Num.

print(sum([1, 2, 3]))   # 6

6

Case: complex

x a CVec or CMat. Returns a Complex scalar instead of a Num — the one reduction in this table that changes its own return *type* based on its input, rather than just its shape.

print(sum([1 + 1j, 2 + 2j]))   # 3 + 3j

3 + 3j

Overloads: max

print(max([1, 5, 3]))      # 5 -- one-argument reduction

print(max([1, 5, 3], 4))   # [4, 5, 4] -- two-argument elementwise clip

5
[4, 5, 4]

Overloads: min

print(min([1, 5, 3]))      # 1 -- one-argument reduction

print(min([1, 5, 3], 4))   # [1, 4, 3] -- two-argument elementwise clip

1
[1, 4, 3]

The next example runs the everyday descriptive-statistics toolkit — mean, std, iqr, and percentile — over a batch of 1000 standard-normal draws from randn, the kind of quick numeric summary a script reaches for before deciding whether a distribution looks the way it should. percentile(x, 95) is included specifically to show the 0-100 scale, as opposed to quantile's 0-1 scale for the same computation.

samples = randn(1000, 1)
print("mean={mean(samples)}  std={std(samples)}  iqr={iqr(samples)}")
print("p95={percentile(samples, 95)}")

mean=0.012612  std=1.005695  iqr=1.356399
p95=1.671269

The second batch below covers the functions that don't fit neatly into the first example: prod (multiplicative reduction), the whole arg* index family (argmedian, argquantile, argmean — each returning a *position* into the data rather than a value), clamp's three-argument clipping, and the var/std relationship (var is always std squared, which the last line checks directly), plus lse as the numerically stable form of log(sum(exp(x))).

v = [5.0, 3.0, 8.0, 1.0, 9.0]
print("prod={prod(v)}")
print("argmedian={argmedian(v)}  argquantile(0.9)={argquantile(v, 0.9)}  argmean={argmean(v)}")
print("clamp={clamp([-5, 0, 5, 10], 0, 5)}")
print("var={var(v)}  std^2={std(v)^2}")
print("lse={lse(v)}")

prod=1080
argmedian=0  argquantile(0.9)=4  argmean=0
clamp=[0, 0, 5, 5]
var=11.2  std^2=11.2
lse=9.328591

Shape, Construction & Introspection

Overloads: zeros

print(zeros(3))                              # [0, 0, 0]

print(zeros(2, 3))                           # [0, 0, 0; 0, 0, 0]

filt = butter(4, "low", 100, 1000)
print(zeros(filt))                           # [-1, -1, -1, -1] -- a CVec, not a fresh zero-fill

[0, 0, 0]
[0, 0, 0; 0, 0, 0]
[-1, -1, -1, -1]

Overloads: reshape

v = [1, 2, 3, 4, 5, 6]
print(reshape(v, 2, 3))        # [1, 3, 5; 2, 4, 6] -- vector source, real relayout

[1, 3, 5; 2, 4, 6]

reshape(A, 3, 2) on the general 2×3 matrix A = [1, 2, 3; 4, 5, 6] from the earlier example, by contrast, does not relayout the data — it hits the broadcast check and errors:

qu: runtime error: shapes (2x3) and (3x2) do not broadcast (each axis must match or be 1)

The example below builds a small matrix and reads its shape back with size/numel, constructs identity matrices with both spellings (eye/identity), and shows the two cat orientations side by side — dim=1 stacks the two vectors as new rows, dim=2 lays them end to end as one longer row — the same distinction vstack/hstack make under more MATLAB-flavored names.

A = [1, 2, 3; 4, 5, 6]
print("size={size(A)}  numel={numel(A)}")
print(eye(3))
print(identity(2))
p = [1, 2]
q = [3, 4]
print(cat(1, p, q))   # vertical: stacks rows

print(cat(2, p, q))   # horizontal: concatenates

size=[2, 3]  numel=6
[1, 0, 0; 0, 1, 0; 0, 0, 1]
[1, 0; 0, 1]
[1, 2; 3, 4]
[1, 2, 3, 4]

-> reshape operator

expr -> (rows, cols) is sugar for reshape(expr, rows, cols) — rewritten to that exact call at parse time (no separate evaluation path, so results are always identical to calling reshape directly). Scoped to the 2-argument form only (matching reshape's own signature); there is no 1-argument -> n shorthand. Because reshape itself only truly relays out data for a vector-like source (see the reshape row above), chaining two reshapes only works when the intermediate shape keeps one dimension equal to 1 — reshaping a genuinely two-dimensional intermediate result into yet another two-dimensional shape hits the same broadcast-only limitation and errors.

Precedence: -> sits *above* |> (pipe) but *below* as/??/the ternary ? : — it takes the entire preceding pipe/arithmetic expression as its value operand, not just the nearest operand:

v = [1, 2, 3, 4, 5, 6]
m = v -> (2, 3)                  # same as reshape(v, 2, 3)

print(m)
mp = df.price -> (n, 1)          # field access resolves first: (df.price) -> (n, 1)

print(mp)
s = a + b -> (r, c)              # reshapes the SUM: (a + b) -> (r, c)

print(s)
r2 = x |> normalize -> (r, c)    # reshapes the PIPED RESULT: (x |> normalize) -> (r, c)

print(r2)
chained = v -> (6, 1) -> (2, 3)  # left-associative: reshape twice, via an intermediate column shape

print(chained)

[1, 3, 5; 2, 4, 6]
[10; 12.5; 9.75; 14]
[3.512478, 2.311214, -0.06689, 1.375279; -3.006714, 1.368403, -1.154908, -1.612839; -1.413089, 0.622247, -0.150671, -0.127341]
[0, 0, 0, 0; 0, 0, 0, 0; 0, 0, 0, 0]
[1, 3, 5; 2, 4, 6]

Rearranging Arrays

The example below runs the identity operation val on a plain number (to show it really is just a pass-through today), transposes a small matrix, and reverses a vector with mirror — the three "rearrange, don't recompute" operations in this table that are easiest to check by eye.

A = [1, 2, 3; 4, 5, 6]
x = 5
print("val={val(x)}")
print(transpose(A))
print("mirror={mirror([1, 2, 3])}")

val=5
[1, 4; 2, 5; 3, 6]
mirror=[3, 2, 1]

Overloads: flip

Case: Str

A a Str. Every spelling (flip, flipud, fliplr, mirror) reverses character order the same way, regardless of dim — a string has no row/ column orientation to pick between.

print(flip("hello"))   # olleh

olleh

Case: Vec

A a plain Vec, which likewise has no row/column orientation. Both directions just reverse element order, so dim makes no difference here either.

print(flip([1, 2, 3]))   # [3, 2, 1]

[3, 2, 1]

Case: Mat

A an r×c Mat (or Image). Here dim genuinely matters: 1 reverses row order (top-to-bottom), 2 reverses column order (left-to-right).

M = [1, 2; 3, 4]
print(flip(M, 1))   # [3, 4; 1, 2] -- rows reversed

[3, 4; 1, 2]

Selection

Overloads: where

mask = [true, false, true]
print(where(mask))          # [0, 2] -- indices, 1-argument form

print(where(mask, 1, 0))    # [1, 0, 1] -- elementwise pick, 3-argument form

[0, 2]
[1, 0, 1]

Overloads: find

mask = [true, false, true]
print(find(mask))                 # [0, 2] -- same as where(mask)


is_even(n) := n mod 2 == 0
xs = [1, 3, 4, 5, 6]
print(find(xs, "is_even"))        # 4 -- first element, not an index

[0, 2]
4

The example below builds a small vector, zeroes out its negative entries with a boolean-mask assignment, then uses where to find which (now non-negative) entries exceed a threshold and pulls just those values back out with x[idx] — the round trip of "get the indices, then use them to index back into the data" that where exists for.

x = [-2, 0.5, 4, -1, 5, 2]
x[x < 0] = 0
print(x)              # [0, 0.5, 4, 0, 5, 2]

idx := where x > 3
print(idx)             # [2, 4] -- 0-based positions of the entries above 3

print(x[idx])          # [4, 5] -- the entries themselves

plot(x[idx], "x", color="red")

[0, 0.5, 4, 0, 5, 2]
[2, 4]
[4, 5]
· plot: figure 1 (2 points)

find's two-argument form searches by a named predicate instead of a mask — useful when the condition isn't a simple elementwise comparison but a small function of its own. The example below defines is_even as an ordinary one-argument Qu function, then hands its *name* (as a string) to find to get back the first element of xs that satisfies it.

is_even(n) := n mod 2 == 0
xs = [1, 3, 4, 5, 6]
print(find(xs, "is_even"))   # 4 -- first element for which is_even is true

4

Matrix Algebra & Decompositions

Decompositions that return more than one related result (svd, qr, lu, eig) come back as a Value::Model handle with named fields — the same convention the fitted-model protocol uses (see Statistics & Machine Learning) — rather than a tuple, since Qu has no multi-return values.

Unlike pinv/inv/det/svd, matmul does not accept a CMat — calling it on a complex matrix errors with "cannot use complex matrix as a matrix" rather than producing a complex product; there is no complex overload to reach for here.

One convention worth knowing before the example: a bracket literal like [1, 1, 0, 0] is a column vector (4×1) in Qu, the opposite of MATLAB's row-vector default for the same syntax — reshape with -> (1, n) to get a row. The example below builds the row/column pair, multiplies them one way to get the 1×1 inner product (the dot product), then the other way to get the 4×4 outer product — the same progression MATLAB's own mtimes documentation uses, and the same two numbers it gets — before moving on to two same-shape square matrices multiplied via both names matmul/mtimes are known by.

row = [1, 1, 0, 0] -> (1, 4)
col = [1; 2; 3; 4]
print(matmul(row, col))    # 3 — inner product (dot product) of the two

print(matmul(col, row))    # 4x4 outer product: [1,1,0,0; 2,2,0,0; 3,3,0,0; 4,4,0,0]


A = [1, 2; 3, 4]
B = [5, 6; 7, 8]
print(matmul(A, B))
print(mtimes(A, B))   # matmul and mtimes are the same builtin under two names

3
[1, 1, 0, 0; 2, 2, 0, 0; 3, 3, 0, 0; 4, 4, 0, 0]
[19, 22; 43, 50]
[19, 22; 43, 50]

gpu_matmul(A, B) has the same signature and result as matmul, but only runs in a build compiled with --features gpu (and dispatches to a real GPU device at runtime); a normal build errors clearly if called.

The next example demonstrates two of the heavier decompositions: chol, which factors a symmetric positive-definite matrix into a lower-triangular square root (useful for generating correlated random samples or solving normal equations efficiently), and svd, whose singular values summarize how well- or ill-conditioned a matrix is. Multiplying the Cholesky factor by its own conjugate transpose should reproduce the original matrix exactly, which is the cheapest possible sanity check on the factorization.

A = [4, 1; 1, 3]
L = chol(A)
print(L * L.ctranspose())   # reproduces A


r = svd(A)
print(r.s)   # singular values

[4, 1; 1, 3]
[4.618034, 2.381966]

Seeing what a matrix does

cond(A) is a number, and a number is hard to have intuition about. It is the ratio of the longest to the shortest axis of the ellipse A turns the unit circle into, so it can just be drawn:

A = [4, 1, 1, 3] as matrix(2, 2)
print("cond(A) = {round(cond(A), 3)}")

angles = linspace(0, 2 * pi, 200)
unit_x = cos(angles)
unit_y = sin(angles)
out_x = zeros(200)
out_y = zeros(200)
for i in 0 to 199
    v = A * [unit_x[i], unit_y[i]]
    out_x[i] = v[0]
    out_y[i] = v[1]
end for

plot(unit_x, unit_y, color = "#94a3b8", label = "unit circle")
plot(out_x, out_y, color = "#4169e1", width = 2, label = "A applied to it")
legend()
title("what a 2x2 matrix does to a circle")

cond(A) = 1.939
· plot: figure 2 (200 points)
· plot: figure 3 (200 points)

A well-conditioned matrix gives a nearly round ellipse. As the condition number grows the ellipse flattens, and a solve loses digits along the short axis — which is what "ill-conditioned" means, drawn.

Overloads: svd

Case: real

A a real r×c Mat. Returns a Model with fields u (r×k Mat), s (length-k Vec of singular values), and vt (k×c Mat), where k = min(r,c) — the full factorization, ready to reassemble A.

A = [4, 1; 1, 3]
r = svd(A)
print(r.s)   # [4.618034, 2.381966]

[4.618034, 2.381966]

Case: complex

A a CMat. Returns a *differently-shaped* Model — only values (the singular values), condition, and rank; there is no u/vt factor pair to reassemble A from in this branch.

C = [1 + 1j, 2; 3, 4 - 1j]
rc = svd(C)
print("values={rc.values}  condition={round(rc.condition, 3)}  rank={rc.rank}")

values=[5.628889, 0.561794]  condition=10.019  rank=2

More functions

The example below runs the small-matrix toolkit that doesn't fit neatly into the decomposition table above: diag/trace/cond read three different summary numbers off the same 2×2 matrix, and nnls fits a non-negative least-squares solution to a tiny system, showing both the fitted x (each entry constrained >= 0) and the leftover residual.

A = [4.0, 1.0; 1.0, 3.0]
print("diag={diag(A)}  trace={trace(A)}  cond={cond(A)}")

b = [1.0, 2.0]
r = nnls(A, b)
print("nnls x={r.x}  residual={r.residual}")

diag=[4, 3]  trace=7  cond=1.938749
nnls x=[0.090909, 0.636364]  residual=4.965068e-16

This next example covers three unrelated but frequently-needed shape utilities together: meshgrid builds the pair of coordinate matrices a contour or surface plot needs from two plain axis vectors, ones_like/ zeros_like stamp out a filler array that automatically matches another array's shape, and row_mean/row_sum collapse a matrix down to one number per row.

g = meshgrid([1, 2, 3], [1, 2])
print(g.x)   # x-coordinate at every grid point

print(g.y)   # y-coordinate at every grid point


M = [1, 2, 3; 4, 5, 6]
print(ones_like(M))
print(zeros_like(M))
print("row_mean={row_mean(M)}  row_sum={row_sum(M)}")

[1, 2, 3; 1, 2, 3]
[1, 1, 1; 2, 2, 2]
[1, 1, 1; 1, 1, 1]
[0, 0, 0; 0, 0, 0]
row_mean=[2; 5]  row_sum=[6; 15]

Finally, the length-unit helpers and the fifo occupancy checks — two unrelated conveniences that happen to round out this chapter's "more functions" table. mm/cm/inch/pt convert a physical length into the figure's internal coordinate units (used anywhere a plotting call wants an exact size, like legend(fontsize = pt(7))); is_empty/ is_full let a script check a fifo ring buffer's occupancy before deciding whether to push into it.

print("mm(10)={mm(10)}  cm(1)={cm(1)}  inch(1)={inch(1)}")   # e.g. figure_size(mm(180), mm(120))


f = fifo(2)
print("is_empty={is_empty(f)}")
push(f, 1)
push(f, 2)
print("is_full={is_full(f)}")   # at capacity -- the next push would evict the oldest

mm(10)=100  cm(1)=100  inch(1)=254
is_empty=true
is_full=true

Overloads: diag

Case: matrix

M an r×c Mat. Reads the main diagonal off M and returns it as a length-min(r,c) Vec — extraction.

M = [1, 2; 3, 4]
print(diag(M))   # [1, 4]

[1, 4]

Case: vector

v a length-n Vec (or scalar/other coercible value). Builds a new n×n Mat with v's entries on the diagonal and zeros elsewhere — construction, the opposite direction from the matrix case.

print(diag([1, 2, 3]))   # [1, 0, 0; 0, 2, 0; 0, 0, 3]

[1, 0, 0; 0, 2, 0; 0, 0, 3]