Random numbers — rml internals

In short: Generator encrypts a 128-bit counter with a 64-bit key (the seed) through ten Philox rounds and hands out the four 32-bit lanes of each block as words. Skipping ahead is O(1) and so is saving or restoring state. On top of the words sit exact conversions to float and double, unbiased integers by rejection, and Box-Muller normals. distributions.h adds sampled tensors, Fisher-Yates permutations, and normal, Bernoulli and categorical distributions with log-probability, entropy and KL divergence in nats.

The headers link reproducibility, probability and distributions and entropy and KL divergence on /ml-models. Every sampler in rml takes a random::Generator&, so a fixed seed makes a whole run reproducible (who draws from it).

Why a counter-based generator

A conventional generator such as std::mt19937 carries a large state and steps it once per draw, so draw n exists only after draws 0 through n − 1. Philox (Salmon et al., "Parallel Random Numbers: As Easy as 1, 2, 3", SC 2011) is instead a keyed bijection on 128-bit counters (philox.h L4–12). Block i of the stream is Philox(counter = i, key = seed), computed directly with nothing carried over from block i − 1. Three things follow:

Any element of a stream can be produced on its own, so Skip costs nothing.
Work split across threads can draw the same numbers a serial run would.
A generator's entire state is (seed, stream, position): three uint64s, trivially saved in a checkpoint.

The word layout from philox.h L25–27. A different stream with the same seed gives an independent sequence.

Philox4x32-10

Philox4x32 (philox.cc L33–42) runs ten rounds. Each round (L24–29) multiplies two of the four lanes by fixed odd constants, keeps both halves of each 64-bit product, and mixes them with the other two lanes and the key:

Output lane	Value
0	hi(M1 · c2) ⊕ c1 ⊕ k0
1	lo(M1 · c2)
2	hi(M0 · c0) ⊕ c3 ⊕ k1
3	lo(M0 · c0)

M0 = 0xD2511F53 and M1 = 0xCD9E8D57 (L12–16). Between rounds the key is bumped by two Weyl constants, 0x9E3779B9 (the golden ratio) and 0xBB67AE85 (√3 − 1), each as a 32-bit fraction. The high half of each product depends on every bit of its input lane, which is where the mixing comes from. The implementation is checked against three known-answer vectors from the Random123 reference distribution.

Generator

Generator holds the seed, the stream, the number of words consumed, and the current four-word block. Refill (philox.cc L46–52) recomputes the block for position / 4. Everything else is arithmetic on position:

Method	What it returns	Words used	Source
`NextU32`	The next lane, refilling when the position is a multiple of 4.	1	L54
`NextU64`	Two words, the first as the high half.	2	L59
`NextFloat`	The top 24 bits × 2⁻²⁴: uniform in [0, 1), and every value is an exact float.	1	L64
`NextDouble`	The top 53 bits of a u64 × 2⁻⁵³: uniform in [0, 1) at full double precision.	2	L66
`UniformInt(n)`	Uniform in [0, n) with no modulo bias (below). n = 0 aborts.	2 or more	L68
`NextNormal`	A standard normal by Box-Muller (below).	4	L78
`Skip(words)`	Advances the position in O(1), refilling only if it lands mid-block.	—	L85
`state()`, `set_state`	The (seed, stream, position) triple, and restoring it.	—	L90

Unbiased integers

Taking x % n of a uniform 64-bit x would make the smallest residues slightly more likely whenever n does not divide 2⁶⁴. UniformInt rejects the bottom 2⁶⁴ mod n values of the range, computed in unsigned arithmetic as -n % n, and redraws. What remains is a whole multiple of n, so every residue is equally likely. For n = 3 the threshold is 1: only x = 0 is rejected, a chance of 1 in 2⁶⁴, so a redraw almost never happens.

Normals

Box-Muller turns two uniforms into a standard normal: z = √(−2 ln u₁) · cos(2π u₂). NextNormal takes u₁ = 1 − NextDouble(), which lies in (0, 1], so the logarithm is always finite. It computes in double and rounds to float. The paired sine value is not kept, so each normal costs two doubles (four words), and a draw never depends on whether an earlier call left a cached value behind.

Distributions

distributions.h fills tensors from a Generator in row-major element order, so a fixed seed gives a fixed tensor. Entropy and KL divergence are in nats.

API	Behavior	Source
`Uniform(shape, gen, low, high)`	low + (high − low) · `NextFloat()` per element. Requires low ≤ high.	L47
`Normal(shape, gen, mean, stddev)`	mean + stddev · `NextNormal()`. Requires stddev ≥ 0.	L55
`Permutation(n, gen)`	Fisher-Yates: from the end, swap each slot i with a uniform j in [0, i]. Every one of the n! orders is equally likely.	L63
`NormalDistribution`	Scalar mean and stddev (stddev > 0). `LogProb` = −z²/2 − ln σ − ln(2π)/2, entropy = ln(σ√(2πe)).	L76
`Bernoulli`	One probability per element, checked to be in [0, 1]. `FromLogits` applies a sigmoid. Samples are 0 or 1 as f32 (`NextFloat() < p`). `LogProb` uses `log1p(−p)` for outcome 0 and rejects outcomes other than 0 and 1.	L112
`Categorical`	Over the last axis; leading axes are a batch. The constructor normalizes each row (non-negative, nonzero sum). `FromLogits` is a softmax. `Sample` returns i32 indices with the batch shape, and `LogProb` takes them back.	L170
`KlDivergence(p, q)`	For each pair of the same family. Zero exactly when p = q, and not symmetric.	.h L85

Categorical sampling

Categorical::Sample (L192–213) inverts the CDF. It draws u = NextDouble() and walks the row, adding probabilities, until u falls below the running sum. Zero-probability categories are skipped, so they can never be chosen. If rounding leaves u above the final sum, the last category with nonzero probability is taken. For the row [0.25, 0.5, 0, 0.25], u = 0.6 passes 0.25, stops at 0.75, and returns index 1.

Entropy and KL edge cases

Two helpers (L23–31) carry the limits that keep entropy and KL finite where they should be. x log x is taken as 0 at x = 0, so a category with probability 0 adds nothing to the entropy. x log(x / y) is 0 when x = 0, and infinity when x > 0 but y = 0. That is, KL(p ‖ q) is infinite when p puts mass where q has none. Both sums are accumulated in double.

Family	KL(p ‖ q)	Source
Normal	ln(σ_q / σ_p) + (σ_p² + (μ_p − μ_q)²) / 2σ_q² − 1/2	L102
Bernoulli	p ln(p/q) + (1 − p) ln((1 − p)/(1 − q)), per element	L156
Categorical	Σₖ pₖ ln(pₖ / qₖ) over the last axis	L247

Who draws from it

No module keeps a hidden global generator. Every caller takes a random::Generator& or owns one, so the seed is always visible at the call site:

Caller	What it samples	Source
Weight init	Xavier and Kaiming, uniform and normal; the `Linear`, `Conv2d`, `Embedding`, attention and recurrent constructors take a generator.	nn/init.h L36
Dropout	Owns a generator seeded at construction and samples a `Bernoulli(1 − p)` keep mask.	nn/dropout.cc L15
Data splits	Train/test, stratified and k-fold splits.	data/split.h L25
Variational losses	`ReparameterizedSample`: z = mean + exp(log_var / 2) · ε with ε ~ N(0, I).	losses/variational.h L30
Classical ML	k-means++ seeding, and feature subsampling when a tree grows (`GrowTree` draws only then).	classical/kmeans.h L64
LLM sampler	Owns a generator for temperature, top-k, top-p and min-p sampling.	llm/sampler.h L80

Those modules are covered on Neural network layers, Data and metrics, Classical ML and Generation.

Tests

File	Tests	What they pin down
`philox_test.cc`	8	Three Random123 known-answer vectors; words equal to the counter blocks for a given seed and stream; the same seed and stream repeating, and a different seed or stream not; `Skip` and state restore matching sequential draws; floats in [0, 1) with mean 0.5; `UniformInt` filling bins evenly and aborting on 0; normal mean, variance and fourth moment (3, a Gaussian's kurtosis).
`distributions_test.cc`	13	Uniform range, moments and seeding; normal moments; permutations hold every index once and all 6 orders of 3 elements appear about equally often in 36,000 draws; normal log-probability, entropy and KL against the closed form and a Monte Carlo estimate; Bernoulli samples, log-probability (−∞ for an impossible outcome), entropy, logits and KL; categorical normalization, sampling in proportion, log-probability, entropy, KL and softmax logits; invalid inputs aborting.

Randomized tests elsewhere in rml use fixed seeds of this generator, so their failures reproduce. How these run under sanitizers is on Testing.