feat(random): numpy.random with NumPy's numbers for the same seed

[tannevaled]

What

A new subpackage, random: numpy.random on go-ndarray arrays, giving NumPy's numbers for the same seed, so a notebook migrated from NumPy keeps its results.

rng := random.DefaultRNG(42)          // np.random.default_rng(42)
x, _ := rng.StandardNormal(3, 4)      // same 12 values as NumPy
random.Seed(0); u, _ := random.Rand(3) // np.random.seed(0); np.random.rand(3)

• SeedSequence (pool, generate_state uint32/uint64, spawn, spawn keys, pool size).
• Bit generators PCG64 (default), PCG64DXSM, MT19937, Philox, SFC64, with Jumped/Advance/RandomRaw/Spawn/state.
• Generator: Random (float64/float32), Integers (all 9 integer/bool dtypes, endpoint; Lemire with the buffered 8/16-bit draws), Uniform, ziggurat StandardNormal/Normal and StandardExponential/Exponential (float64 and float32), StandardGamma/Gamma, Beta, ChiSquare, Binomial, Poisson, Geometric, Lognormal, Laplace, Logistic, Choice (replace, p, axis, shuffle), Shuffle, Permutation, Permuted, Bytes, Dirichlet, Multinomial, plus F, noncentral χ²/F, t, Cauchy, Gumbel, Weibull, Pareto, power, Rayleigh, Wald, triangular, von Mises, Zipf, logseries, negative binomial.
• Legacy RandomState (MT19937 legacy seeding, polar Gaussian with its cached value, masked randint, NumPy 1.16 binomial) and the module functions Seed, Rand, Randn, Randint, RandomSample, Normal, Uniform, Choice, Shuffle, Permutation, Binomial, Poisson.

Base: #42 is merged; this targets main (rebased on #43).

Verification

• Fixture from NumPy 2.5.3 (x86-64 build), random/testdata/cases.py, embedded: SeedSequence (10 entropies), raw streams/jumps/advance/spawn of all 5 generators × 6 seeds (0, 1, 42, 2⁶⁴−1, [1,2,3], 2¹²⁸+12345), legacy MT19937 seeding, 2,492 calls / 81,806 values, and 50 long runs by SHA-256 (4.49 M values; ziggurat tails and wedges included).
• Bit-exact: integers, raw output, every arithmetic-only sampler (normal, exponential, binomial, Poisson, geometric, choice, shuffle, legacy rand/randint/binomial/poisson …).
• libm: log/log1p/exp/expm1/pow are implemented here correctly rounded (Ziv fast path ~2⁻⁶⁸ + double-double fallback; 0 misrounded on 63,000 mpmath-checked arguments), deterministic on every arch. Residual: 0.03–0.5% of values one ULP from NumPy, each a NumPy libm misrounding (checked with mpmath); with Go's math it was 3–10%. Bounded per sampler in the test (libmULP), never globally.
• Passes unchanged on darwin/arm64, darwin/amd64 (Rosetta), and linux 386, arm, riscv64, ppc64le, s390x (big-endian; Bytes = NumPy's little-endian bytes) under qemu.
• Negative controls (PCG64 multiplier, SeedSequence MIX_MULT_L/MULT_A, a ziggurat entry, masked↔Lemire, MT seeding position, legacy BTPE signs) all fail the suite; one does not (legacy BTPE constant 13680→13860: a ~2·10⁻⁹ term that flips nothing in 10⁶ draws) — documented.
• CI gate: go test -race -coverpkg=all → 100.0%; go vet (native, linux/386, js/wasm, windows), GOOS=linux GOARCH=s390x go test -c, gofmt -l . empty.

Finding: NumPy's aarch64 build disagrees with its own x86-64 build (FMA contraction, different libm): 20.8% of uniform values, 14% of legacy randn. This package follows the x86-64 (unfused) arithmetic on every architecture.

Performance (Apple M4 Max, loaded host, 7 interleaved rounds, medians, ms per 10⁶)

	Go	NumPy	Go/NumPy
random	3.51	3.13	1.12
random(float32)	1.89	1.64	1.16
integers(0, 1000)	2.31	1.70	1.36
standard_normal	4.12	3.68	1.12
legacy randn	20.36	9.65	2.11

Slower everywhere; details and causes in docs/random.md.

Not done

Broadcast (array) parameters and bounds, out=, hypergeometric, multivariate_normal, multivariate_hypergeometric, the RandomState samplers beyond those listed (their 1.16 versions were not ported), pickling/state dicts as Python objects.

🤖 Generated with Claude Code