Every function is a graph.

Absolute value, built from Relu. Exact: one of the two terms is always zero.

Elementwise square.

Elementwise cube.

Linear interpolation from a to b by a scalar t.

Leaky ReLU with slope alpha for negative inputs, built from Relu.

SiLU (also called swish): x times its own sigmoid.

GELU, tanh approximation, as used in GPT-style networks.

Softsign: a smooth, bounded alternative to tanh. Calls abs.

Clamp every element into [lo, hi] with exact Maximum and Minimum.

ReLU capped at 6, as in mobile networks.

Piecewise-linear sigmoid: relu6(x + 3) / 6. Calls relu6.

Sign of each element (−1, 0 or 1), from two comparison masks. Exact.

Softplus in its overflow-free form. Mathematically identical to log(1 + eˣ), which the reference uses.

ELU: the identity for positive inputs, α(eˣ − 1) otherwise.

Mish: x · tanh(softplus(x)). Calls softplus.

An angle wrapped into [0, 2π): the remainder after dividing by 2π, which takes the divisor's sign, so negative angles come out positive.

The day of the week of a date given as days since 1970-01-01, Monday = 0: that day was a Thursday, so (days + 3) mod 7. Dates before 1970 work too, because the remainder is never negative.

1 for a leap year of the proleptic Gregorian calendar, else 0: divisible by 4, except centuries, except every fourth century.

Days since 1970-01-01 for a date of the proleptic Gregorian calendar, by Howard Hinnant's algorithm: the year is counted from March, so the leap day falls last, and every division is a floor.

The year, month and day of a date given as days since 1970-01-01, by Howard Hinnant's algorithm: floor divisions by the 400-year era, the year of the era, and a year that starts in March.

One step of Euclid's algorithm on pairs stored as two rows: (a, b) becomes (b, a mod b); a finished pair, with b = 0, stays as it is. gcd repeats it.

Whether any pair still has a non-zero second number, for pairs stored as two rows of non-negative integers: 1 or 0. The condition gcd's loop tests.

The greatest common divisor of each pair of positive integers, by Euclid's algorithm: a While loop of euclid_step until every remainder is zero.

The elementwise product of two complex vectors, each stored as a [2, n] array of real parts and imaginary parts: (a + bi)(c + di) = (ac − bd) + (ad + bc)i.

log(1 − e^(−x)) for x > 0, accurate everywhere (Mächler's rule): log(−expm1(−x)) below log 2, log1p(−e^(−x)) above. Each naive form loses its digits on one side.

Inner product of two vectors.

Squared Euclidean norm. Calls dot.

Matrix times vector.

Row vector times matrix.

Outer product: every product xᵢ·yⱼ, by broadcasting. Exact up to one rounding per entry.

Gram matrix of the columns of X.

Bilinear form xᵀ·A·y. Calls matvec, then dot.

Linear combination of two vectors with scalar weights.

Orthogonal projection of x onto the line spanned by u. Calls dot twice.

The part of x orthogonal to u: x minus its projection. Calls proj.

Squared Frobenius norm: the sum of every squared entry.

Euclidean norm. Calls norm_sq.

Scale a non-zero vector to unit length. Calls norm.

Cosine of the angle between two non-zero vectors. Calls dot and norm.

Frobenius norm of a matrix. Calls frobenius_sq.

One Jacobi step for A·x = b, with A split into its diagonal d and the rest R: solve each row for its own unknown, using the others' current values. The body that jacobi_solve runs.

Whether x still misses A·x = b, with A given as its off-diagonal part R and diagonal d: ‖b − R·x − d⊙x‖² above 10⁻²⁰·‖b‖² (a relative residual of 10⁻¹⁰), as 1 or 0. The condition jacobi_solve loops on.

Solve A·x = b by Jacobi iteration for a diagonally dominant A, from x = 0 until the relative residual is below 10⁻¹⁰. The diagonal comes from an Iota mask; the loop is a While.

The solution of A·x = b, by Gaussian elimination with partial pivoting (the Solve primitive).

The inverse of A: Solve against the identity, which is built from two Iotas.

The Cholesky factor of a symmetric positive-definite A: the lower-triangular L with L·Lᵀ = A.

Solve A·x = b for a symmetric positive-definite A the way it is done in practice: factor A = L·Lᵀ, then one forward and one back substitution.

The log-determinant of a symmetric positive-definite A, from its Cholesky factor: twice the sum of the logs of L's diagonal. No determinant is ever formed, so nothing overflows.

The eigenvalues of a symmetric matrix, ascending.

The unit eigenvectors of a symmetric matrix, as columns in ascending order of eigenvalue, each signed so its largest component is positive.

Ridge regression: the coefficients β minimizing ‖X·β − y‖² + λ‖β‖², from the regularized normal equations. λ > 0 keeps them solvable for any X.

The orthonormal factor Q of A = Q·R, for a tall matrix of full column rank: its columns are an orthonormal basis of A's column space.

The triangular factor R of A = Q·R, upper triangular with a positive diagonal.

Least squares: the x minimizing ‖A·x − b‖ for a tall A of full column rank, through QR: R·x = Qᵀ·b. Unlike the normal equations, it does not square A's condition number.

The orthogonal projection onto A's column space, Q·Qᵀ from the QR factorization.

The singular values of A, largest first.

The spectral norm of A: its largest singular value, the most A can stretch a unit vector.

The condition number of A in the 2-norm, largest over smallest singular value: how much A can amplify a relative error.

The nuclear norm of A: the sum of its singular values, the convex stand-in for rank in low-rank recovery.

The Moore–Penrose pseudo-inverse of a tall A of full column rank, from its SVD: V·diag(1/σ)·Uᵀ. It maps b to the least-squares solution.

The orthogonal factor of the polar decomposition A = W·H: the orthogonal matrix nearest to A, U·Vᵀ from the SVD.

The best rank-1 approximation of A (Eckart–Young): σ₁·u₁·v₁ᵀ, from the leading singular triplet.

Squared Euclidean distance between two points. Calls norm_sq.

Squared distances between every point of X and every point of Y, with one matrix product.

Centroid (mean point) of a set of points.

Shift a set of points so that its centroid is the origin.

Euclidean distance between two points. Calls sqdist.

Distances between every point of X and every point of Y. Calls pairwise_sqdist, and clamps at zero before the square root.

Lengths of the segments of a polyline through the points P, in order.

Total length of the polyline through the points P. Calls segment_lengths.

Arc length along the polyline at every point, starting from 0: the running sum of the segment lengths. Calls segment_lengths.

Signed area of a polygon by the shoelace formula: positive when the vertices run counter-clockwise. The next vertex of each one comes from joining P[1:] and P[:1].

One-nearest-neighbour classification: the label of the point closest to q. ArgMin gives the position; a mask built from Iota reads the label there.

k-nearest-neighbour regression: the mean value of the k points closest to q (ties by index). The ranks of the distances become a mask, so k can be an input.

Population variance, two-pass: the mean first, then the mean squared deviation.

Population covariance of two paired samples.

Weighted mean with positive weights. Calls dot.

Population standard deviation. Calls variance.

Sample variance, dividing by n − 1. The count comes from the Size primitive. Calls variance.

Standard scores: each element's distance from the mean, in standard deviations. Calls std.

Population covariance matrix of the columns of X. Calls center; the row count comes from Size.

Pearson correlation of two paired samples. Calls covariance and std.

log Σ eˣ, computed stably by shifting by the maximum first.

Rescale to [0, 1] by the minimum and the maximum.

The q-quantile with linear interpolation (NumPy's default): sort, then weigh each sorted value by a tent 1 − |i − h| at position h = (n − 1)·q.

The median: the 0.5-quantile, the mean of the two middle values when n is even. Calls quantile.

Interquartile range: the 0.75-quantile minus the 0.25-quantile. Calls quantile twice.

Median absolute deviation, a robust spread: the median distance from the median. Calls median and abs.

Mean without the smallest and the largest value, which makes it resistant to one outlier at each end.

Maximum drawdown: the largest fall from a running peak to a later value, as used for price series.

Counts per bin for increasing bin edges, with NumPy's rules: half-open bins, the last one closed, values outside the edges not counted. Each value's bin comes from comparisons, becomes an integer slot with Cast, and is counted with ScatterAdd.

Add each value into its segment: x[k] goes to base[seg[k]], and segments that repeat accumulate. A ScatterAdd.

Reorder the values by the stable order of the keys: ArgSort the keys, Cast the positions to indices, Gather the values.

The rank of each value, 0 for the smallest; ties are ranked in input order. The argsort of the argsort.

The log-density of a multivariate normal N(μ, S) at x, computed the stable way: one Cholesky factorization gives both the log-determinant and, by a triangular solve, the quadratic form.

The Mahalanobis distance of x from μ under covariance S: the length of x − μ after whitening by S's Cholesky factor.

Principal component analysis, as the share of variance along each principal axis, largest first: the eigenvalues of the covariance matrix over their sum.

The principal axes of a data set, as columns, largest variance first: the right singular vectors of the centred data, each signed so its largest component is positive.

The bin each value falls into, for bins of equal width starting at lo: ⌊(x − lo)/width⌋, as an integer index. Values below lo get negative bins; nothing is clipped.

The normal distribution's cumulative probability, Φ((x − μ)/σ), written with erfc so the lower tail keeps its relative accuracy where 1 + erf would cancel to 0.

The logarithm of the binomial coefficient C(k + m, k), from log-gamma: log Γ(k+m+1) − log Γ(k+1) − log Γ(m+1). The coefficient itself would overflow long before its logarithm does.

The log-probability of k events under a Poisson distribution with mean λ: k·log λ − λ − log k!, with log k! as log Γ(k + 1).

The log-density of the gamma distribution with shape a and rate b: a·log b + (a − 1)·log x − b·x − log Γ(a).

The log-density of the beta distribution on (0, 1): (a − 1)·log x + (b − 1)·log(1 − x) − log B(a, b), with log(1 − x) as log1p(−x) and log B from log-gamma.

The log-density of Student's t distribution with ν degrees of freedom: log Γ((ν+1)/2) − log Γ(ν/2) − ½·log(νπ) − ((ν+1)/2)·log1p(x²/ν).

The normal distribution's quantile: the x with Φ((x − μ)/σ) = p, as μ + σ·Φ⁻¹(p).

One standard normal number for each element of x, by inverse-transform sampling: Φ⁻¹ of the Wichmann–Hill uniform draws. Reproducible from its state.

The chi-square distribution's cumulative probability with k degrees of freedom: P(k/2, x/2), the regularized lower incomplete gamma function.

The probability of at most k events under a Poisson distribution with mean λ: 1 − P(k + 1, λ).

The gamma distribution's cumulative probability with shape a and rate b: P(a, b·x).

The beta distribution's cumulative probability: I_x(a, b), the regularized incomplete beta function.

Student's t distribution's cumulative probability with ν degrees of freedom, through the incomplete beta function: ½·I_{ν/(ν+x²)}(ν/2, ½) on the left, one minus that on the right.

The probability of at most k successes in k + m independent trials with success probability p (m ≥ 1): I_{1−p}(m, k + 1).

The chi-square distribution's quantile with k degrees of freedom: the x with F(x; k) = p, as 2·P⁻¹(k/2, p).

The gamma distribution's quantile with shape a and rate b: P⁻¹(a, p)/b.

The beta distribution's quantile: the x with I_x(a, b) = p, the inverse incomplete beta function.

Student's t distribution's quantile with ν degrees of freedom, through the inverse incomplete beta function: with z = I⁻¹_{2·min(p, 1−p)}(ν/2, ½), x = ±√(ν(1 − z)/z), negative below the median.

First differences: each element minus the one before it. One element shorter than x.

Running total: element i is the sum of x₀ … xᵢ.

Running mean: element i is the mean of x₀ … xᵢ. The counts 1, 2, … come from Iota.

1² + 2² + … + n², with n the length of x (its values are not used): Iota counts 1 … n, each count is squared, and the squares are added. This is the intent the family home follows through EML-U, EML-P and EML-NOVA; for x of length 100 it is 338350.

Trapezoid rule with unit spacing: the area under the samples y.

Trapezoid rule over sample points x (any spacing, in order). Calls diff.

Moving average over a window of three, at every position where the window fits.

One-dimensional convolution with a three-tap kernel w, valid positions only — cross-correlation, as deep-learning libraries define it.

One step of an exponential moving average: blend the new value into the running one. The body that ema scans.

Exponential moving average of a series, started at its first value: a Scan of ema_step.

One step of a discounted return: this reward plus γ times the return after it. The body that discounted_returns scans, backwards.

Discounted returns of a reward sequence, as in reinforcement learning: a reverse Scan of discount_step.

One step of Horner's rule: multiply by t, add the next coefficient. The body that horner scans.

Evaluate a polynomial at t by Horner's rule, coefficients from the highest degree down: a Scan of horner_step, keeping the last value.

One forward-Euler step of the linear system x′ = A·x over a time step dt. The body that euler_linear scans.

Integrate x′ = A·x from x₀ by forward Euler over the time steps dts; returns the state after every step. A Scan of euler_step.

Reorder a vector by a permutation: element i of the result is x[perm[i]]. A Gather.

The inverse of a permutation: inv[perm[i]] = i. Each position is scattered to where it came from, then Cast back to integers.

One step of Newton's method for √a: average the guess with a divided by it. The body that sqrt_newton runs.

Whether a square-root guess still misses: |s² − a| above 10⁻¹² of a, as 1 or 0. The condition sqrt_newton loops on.

√a by Newton's method, starting from (a + 1)/2 and stopping once |s² − a| ≤ 10⁻¹²·a: a While over newton_sqrt_step until sqrt_unconverged says stop.

One bisection step for x³ = a: halve the bracket [lo, hi], keeping the half where x³ − a changes sign. The body that cbrt_bisect runs.

Whether a bisection bracket is still too wide: hi − lo above 10⁻¹²·(|a| + 1), as 1 or 0. The condition cbrt_bisect loops on.

∛a by bisection: start from the bracket ±(|a| + 1), which always contains the root, and halve it until it is narrower than 10⁻¹²·(|a| + 1). A While whose state is the bracket itself.

One step of a polynomial rolling hash: multiply by 131, add the next code, reduce modulo the prime 1 000 000 007. Every value stays below 2⁵³, so the arithmetic is exact. The body poly_hash scans.

A polynomial rolling hash of a sequence of byte codes: a Scan of hash_step from 0, keeping the last value. Equal sequences hash equally; it is not a cryptographic hash.

The discrete Fourier transform of a complex vector, unnormalized: yₖ = Σⱼ zⱼ·e^(−2πi·jk/n).

The inverse discrete Fourier transform: zⱼ = (1/n)·Σₖ yₖ·e^(2πi·jk/n), so ifft(fft(z)) = z.

The power spectrum of a real signal: |yₖ|² for its discrete Fourier transform y, the energy at each frequency.

The circular convolution of two real signals through the Fourier transform: transform both, multiply, transform back, keep the real part. Takes n log n operations instead of n².

The circular autocorrelation of a real signal, rₖ = Σⱼ xⱼ·x₍ⱼ₊ₖ₎ mod n, by the Wiener–Khinchin theorem: the inverse transform of the power spectrum.

One step of the Wichmann–Hill generator (AS 183, 1982): three small multiplicative congruential generators side by side, each sᵢ ↦ aᵢ·sᵢ mod mᵢ. It is the body uniform_draws scans; the element it is handed only sets how many steps there are.

One uniform number in [0, 1) for each element of x, from the Wichmann–Hill generator started at state: a Scan of wichmann_hill_step, then the fractional part of s₁/m₁ + s₂/m₂ + s₃/m₃. The same state always gives the same numbers, on every backend. Its period is about 7·10¹²; it is not for cryptography.

A reproducible random permutation of x: draw one uniform number per element, then sort the elements by their numbers.

Affine layer: a batch of rows times a weight matrix, plus a bias.

Two-layer perceptron with a ReLU between the layers. Calls linear twice.

Scaled dot-product attention for one head.

Log of the softmax, computed as x − logsumexp(x) without forming the softmax. Calls logsumexp.

Layer normalisation of one feature vector, with scale γ, shift β and a small ε. Calls variance.

Scaled dot-product attention with the standard 1/√d scale computed inside the graph (Size, then Sqrt). Calls attention.

One step of an Elman recurrent network: the next hidden state from the current one and an input. The body that rnn scans.

Run an Elman recurrent network over a sequence of inputs; returns every hidden state. A Scan of rnn_cell, differentiable end to end.

Classification accuracy: the share of rows whose largest logit is at the position of the largest label entry (one-hot labels).

Causal (decoder) self-attention: row i attends only to rows j ≤ i. The mask compares positions from Iota (j > i is the future); future scores are replaced by the row minimum before the stable softmax, so nothing can overflow, and their weights are set to zero after it.

Embedding lookup: the rows of the table E at the integer ids, in order (repeats allowed).

The mean of the embeddings of a bag of ids: one vector for a whole set of tokens. Calls embedding_lookup.

The posterior mean of a Gaussian process with a squared-exponential kernel of length scale ℓ and noise variance σ², at new points Xs. The training system is solved by Cholesky (calls cholesky_solve).

Mean squared error.

Mean absolute error. Calls abs.

Hinge loss for labels ±1 and real-valued scores.

L2 weight penalty. Calls frobenius_sq.

Cross-entropy between a target distribution and the softmax of logits. Calls log_softmax.

Huber loss: quadratic for small errors, linear beyond δ. Calls abs.

Binary cross-entropy for probabilities strictly between 0 and 1.

Log-cosh loss, in a form that never overflows. Calls abs.

Negative log-likelihood with integer class labels: minus the mean of each row's log-probability at its label.

Softmax cross-entropy over rows of logits with integer labels, the way classifiers are trained: a stable log-softmax per row, then nll_labels.

Binary cross-entropy taken from logits: max(z, 0) − z·t + log1p(e^(−|z|)). It never forms the probability, so it stays finite and accurate where sigmoid would round to 0 or 1.