random — numpy.random, the same numbers

Mirrored from docs/random.md at v0.27.0.

Package random is numpy.random over go-ndarray arrays. Its goal is migration: a notebook moved from NumPy to Go keeps its numbers for the same seed.

rng := random.DefaultRNG(42)          // np.random.default_rng(42)
x, _ := rng.StandardNormal(3, 4)      // rng.standard_normal((3, 4))
k, _ := rng.Integers(0, 10, 5)        // rng.integers(0, 10, 5)

random.Seed(0)                        // np.random.seed(0)
u, _ := random.Rand(3)                // np.random.rand(3)

Example, ExampleSeed and ExampleGenerator_ChoiceN in example_test.go print NumPy’s own numbers.

What is implemented

The size comes last; no size gives a 0-d array holding what NumPy returns as a scalar (size=None), which draws exactly the same values. ChoiceOptions' zero value is NumPy’s defaults (NoReplace, NoShuffle are the negations).

Not implemented: array-valued (broadcast) distribution parameters and integer bounds, out=, hypergeometric, multivariate_normal, multivariate_hypergeometric, the RandomState samplers not listed above (beta, f, pareto …: NumPy 1.16’s versions differ from Generator’s and were not ported), tomaxint, random_integers, pickling, and float16.

Verification

random/testdata/cases.py runs NumPy 2.5.3 and records:

• SeedSequence: the pool, 9 uint32 and 5 uint64 state words, and the pools and spawn keys of 4 children, for 10 entropies (0, 42, 2⁶⁴−1, [1,2,3], 2¹²⁸+12345, a 10-word list, [], pool size 8, spawn keys with a 2⁴⁰ entry, a key after a 6-word entropy);
• bit generators: for each of the five and six seeds (0, 1, 42, 2⁶⁴−1, [1,2,3], 2¹²⁸+12345), the seeded state, 100 raw outputs (700 for MT19937, which crosses a regeneration), jumped(1), jumped(3), advance(2¹⁰⁰+7), advance(-3), two spawned children; Philox from an explicit counter and key and across a counter carry; legacy MT19937 seeding from 6 integer and array seeds, and its jump;
• 2,222 scripts, 2,492 calls, 81,806 values: every Generator function on PCG64 with every seed, every dtype and bound case of integers (9 dtypes, open and closed, full ranges, ranges of 2³²−1, 2³², 3·2³⁰, 3·2⁶², empty ranges); a subset on the four other bit generators; an interleaved script per generator and seed (the buffered upper half of a 64-bit draw outlives a call); RandomState on 6 seeds (including a 700-word array);
• 50 long runs as SHA-256 digests (4.49 million values): random, standard_normal and standard_exponential (both dtypes, 200,000 each, which reaches the ziggurat tails and wedges), integers (Lemire with rejections, 8- and 16-bit buffering, bool), binomial (inversion and BTPE), poisson (both methods), choice with and without replacement, and RandomState’s rand, randint, binomial, poisson, uniform.

The fixture is embedded, so CI’s emulated targets check it too. It passes, unchanged, on darwin/arm64, darwin/amd64 (Rosetta), and linux 386, arm, riscv64, ppc64le and s390x (big-endian: Bytes gives NumPy’s little-endian bytes) under qemu.

Integers, raw output, and every sampler that is plain arithmetic match NumPy to the bit — including the ziggurat normal and exponential (tails included, over 200,000 draws), binomial, Poisson, geometric, choice, shuffle and the legacy randint/rand/binomial/poisson.

Negative controls

Each mutation below was made on purpose and the suite run:

Where the numbers can differ: the C library

Some samplers call log, exp, log1p, expm1, pow, cos or acos. NumPy takes them from the C library; Go’s math package does not round them the same way. Measured against the libm of the x86-64 NumPy used for the fixture (Apple’s), Go’s math differs in 6.5% (log) to 71% (pow) of random arguments. Used as-is, that made 3% of np.random.randn values and 10% of lognormal values one ULP off.

So random rounds log, log1p, exp, expm1 and pow correctly (libm.go): a fast evaluation to about 2⁻⁶⁸ with Ziv’s rounding test, and a double-double fallback (about 2⁻¹⁰⁰) for the one argument in ~2¹⁵ that needs it. Checked against mpmath at 300 bits on 63,000 arguments (0 misrounded; the committed libm.json keeps 8,400), and the fast paths against the fallback on 400,000 more. The code converts every product explicitly, so Go cannot fuse it into an FMA, and the result is the same on every architecture.

The C libraries round these functions correctly almost always — Apple’s x86-64 libm misrounds 0.06% (log) to 0.3% (exp) of arguments; pow(x, 0.5) is misrounded for 0.19% (it is not sqrt). What remains, on 10⁶ draws against x86-64 NumPy, seed 0:

¹ Large ULP counts are absolute errors of about one ULP of the inner value, measured relative to a result near zero (1 + 3z for z ≈ −1/3; log(-log(U)) for U near 1/e; an angle wrapped to near 0).

Each of the residual differences checked by hand (with mpmath) is a value NumPy’s libm rounded incorrectly. vonmises keeps Go’s cos and acos: Apple’s cos itself misrounds 17% of arguments in [0, π], so rounding correctly would not bring it closer. No difference has ever changed an accept/reject decision — every integer-valued sampler matches exactly — and the test bounds each float sampler’s ULP distance individually (libmULP in numpy_test.go).

NumPy does not agree with itself

The same fixture generated by NumPy 2.5.3’s aarch64 macOS build differs from its x86-64 build: the C compiler fuses a*b + c into one FMA there (clang’s default -ffp-contract=on), and Apple’s aarch64 libm rounds differently.

This package follows the C source operation by operation without fusing, which is what NumPy’s x86-64 builds (Linux, Windows, Intel Macs) compute, and it gives the same numbers on every architecture. A notebook run on an Apple silicon Mac or an aarch64 Linux box will see NumPy’s last bits move in those samplers; that is NumPy, not a migration error.

Deliberate differences

Performance

Apple M4 Max, 16 cores, shared and loaded (load average 10.7–14.4 during the runs). Seven interleaved rounds: the Go benchmark (go test -bench, 20 iterations) then NumPy 2.5.3 aarch64 (median of 20 calls); medians of the rounds, in milliseconds for 10⁶ values:

Go is slower everywhere: 12–36% for the Generator, 2.1× for the legacy randn. The Generator gap is the 128-bit LCG step (3.3 ns per output, a dependent multiply chain) plus zeroing the result (ZerosOf; NumPy’s np.empty does not); the default PCG64 already gets direct, inlinable calls for these four functions. randn pays for a correctly rounded log on every pair (6 ns against 3.4 ns for math.Log) and for MT19937 behind an interface. The samplers that call pow/exp per draw cost about what math.Pow/math.Exp would (19.5 vs 18.4 ns; 6.3 vs 3.9 ns).

Regenerating

uv run --python cpython-3.12-macos-x86_64 --with numpy==2.5.3 python random/testdata/cases.py > random/testdata/cases.json
uv run --with mpmath python random/testdata/libm.py > random/testdata/libm.json
uv run --with mpmath python random/testdata/libm_tables.py > random/zz_libm_tables.go
python3 random/testdata/zig2go.py <numpy>/numpy/random/src/distributions/ziggurat_constants.h random/zz_ziggurat.go

random/testdata/quant.py and TestQuantifyLibm (run with RANDOM_QUANT_DIR) reproduce the 10⁶-draw table.