PRNGs (PCG, xoshiro)

Deep dive

How it works

In: a seed (expanded with SplitMix64), optionally a stream id or jump count · Out: a stream of uniform 32- or 64-bit words (then integers in [0, n), floats in [0, 1)) · State: 16 bytes (PCG32: state + odd increment) or 32 bytes (xoshiro256**) · Cost: about 1 ns per output; one multiply (PCG) or a few shifts, xors and rotates (xoshiro)

PCG32 (XSH-RR 64/32): state ← state·6364136223846793005 + inc (mod 2⁶⁴). Output from the old state: xorshifted = ((old >> 18) ^ old) >> 27 as u32, rot = old >> 59, result = xorshifted.rotate_right(rot). The top 5 bits choose the rotation, so every output bit depends on the strong high bits of the LCG. Different odd increments give different sequences (streams). Jump-ahead by k steps costs O(log k) via the LCG's affine composition.

xoshiro256**: result = rotl(s₁·5, 7)·9; then t = s₁ << 17; s₂ ^= s₀; s₃ ^= s₁; s₁ ^= s₂; s₀ ^= s₃; s₂ ^= t; s₃ = rotl(s₃, 45). The linear engine has a full period; the ** scrambler removes the linear artifacts from the output bits. jump() applies a precomputed polynomial (the characteristic polynomial power) by xoring together states selected by the jump constant bits.

SplitMix64 (seeding): z = (x += 0x9e3779b97f4a7c15); z = (z ^ (z >> 30))·0xbf58476d1ce4e5b9; z = (z ^ (z >> 27))·0x94d049bb133111eb; return z ^ (z >> 31). Each call gives a well-mixed 64-bit word, so four calls fill a xoshiro state with no zero-state risk in practice.

Bounded integers (Lemire): m = x·n as a 128-bit (or 64-bit for 32-bit x) product; the high half is the candidate. The low half l < n means the value might be biased; then reject if l < (2^w − n) mod n. The modulo is computed only in that rare case.

Code

pub struct Pcg32 { state: u64, inc: u64 }

impl Pcg32 {
    const MUL: u64 = 6364136223846793005;
    /// Same seeding as the reference pcg32_srandom_r(seed, seq).
    pub fn new(seed: u64, seq: u64) -> Self {
        let mut r = Pcg32 { state: 0, inc: (seq << 1) | 1 };
        r.next_u32();
        r.state = r.state.wrapping_add(seed);
        r.next_u32();
        r
    }
    #[inline] pub fn next_u32(&mut self) -> u32 {
        let old = self.state;
        self.state = old.wrapping_mul(Self::MUL).wrapping_add(self.inc);
        let xorshifted = (((old >> 18) ^ old) >> 27) as u32;
        xorshifted.rotate_right((old >> 59) as u32)
    }
    /// Advance by `delta` steps in O(log delta) (Brown, "Random number generation with arbitrary strides").
    pub fn advance(&mut self, mut delta: u64) {
        let (mut acc_mult, mut acc_plus, mut cur_mult, mut cur_plus) = (1u64, 0u64, Self::MUL, self.inc);
        while delta > 0 {
            if delta & 1 == 1 {
                acc_mult = acc_mult.wrapping_mul(cur_mult);
                acc_plus = acc_plus.wrapping_mul(cur_mult).wrapping_add(cur_plus);
            }
            cur_plus = cur_mult.wrapping_add(1).wrapping_mul(cur_plus);
            cur_mult = cur_mult.wrapping_mul(cur_mult);
            delta >>= 1;
        }
        self.state = acc_mult.wrapping_mul(self.state).wrapping_add(acc_plus);
    }
    /// Unbiased integer in [0, n) (Lemire 2019), n > 0.
    pub fn below(&mut self, n: u32) -> u32 {
        let mut m = self.next_u32() as u64 * n as u64;
        if (m as u32) < n {
            let threshold = n.wrapping_neg() % n; // (2^32 - n) mod n
            while (m as u32) < threshold { m = self.next_u32() as u64 * n as u64; }
        }
        (m >> 32) as u32
    }
}

pub fn splitmix64(x: &mut u64) -> u64 {
    *x = x.wrapping_add(0x9e3779b97f4a7c15);
    let mut z = *x;
    z = (z ^ (z >> 30)).wrapping_mul(0xbf58476d1ce4e5b9);
    z = (z ^ (z >> 27)).wrapping_mul(0x94d049bb133111eb);
    z ^ (z >> 31)
}

#[derive(Clone)]
pub struct Xoshiro256ss { s: [u64; 4] }

impl Xoshiro256ss {
    pub fn from_state(s: [u64; 4]) -> Self { assert!(s != [0; 4], "all-zero state is a fixed point"); Xoshiro256ss { s } }
    pub fn seed(seed: u64) -> Self {
        let mut x = seed;
        Xoshiro256ss { s: [splitmix64(&mut x), splitmix64(&mut x), splitmix64(&mut x), splitmix64(&mut x)] }
    }
    #[inline] pub fn next_u64(&mut self) -> u64 {
        let s = &mut self.s;
        let result = s[1].wrapping_mul(5).rotate_left(7).wrapping_mul(9);
        let t = s[1] << 17;
        s[2] ^= s[0];
        s[3] ^= s[1];
        s[1] ^= s[2];
        s[0] ^= s[3];
        s[2] ^= t;
        s[3] = s[3].rotate_left(45);
        result
    }
    /// Equivalent to 2^128 calls of next_u64: gives 2^128 non-overlapping streams.
    pub fn jump(&mut self) {
        const JUMP: [u64; 4] = [0x180ec6d33cfd0aba, 0xd5a61266f0c9392c, 0xa9582618e03fc9aa, 0x39abdc4529b1661c];
        let mut acc = [0u64; 4];
        for &j in &JUMP {
            for b in 0..64 {
                if (j >> b) & 1 == 1 { for i in 0..4 { acc[i] ^= self.s[i]; } }
                self.next_u64();
            }
        }
        self.s = acc;
    }
    /// Uniform f64 in [0, 1) with 53 random bits.
    #[inline] pub fn next_f64(&mut self) -> f64 { (self.next_u64() >> 11) as f64 * (1.0 / (1u64 << 53) as f64) }
    /// Unbiased integer in [0, n), n > 0 (Lemire, 128-bit product).
    pub fn below(&mut self, n: u64) -> u64 {
        let mut m = self.next_u64() as u128 * n as u128;
        if (m as u64) < n {
            let threshold = n.wrapping_neg() % n;
            while (m as u64) < threshold { m = self.next_u64() as u128 * n as u128; }
        }
        (m >> 64) as u64
    }
    /// Fisher-Yates shuffle.
    pub fn shuffle<T>(&mut self, v: &mut [T]) {
        for i in (1..v.len()).rev() { let j = self.below(i as u64 + 1) as usize; v.swap(i, j); }
    }
}

Complexity and constants

• Measured with the code above (Apple M4 Pro, -O, summing outputs in a loop): PCG32 1.0 ns per 32-bit output (the 64-bit multiply is a serial dependency), xoshiro256** 0.76 ns per 64-bit output, SplitMix64 0.50 ns, below(n) 0.84 ns. 3·10⁶ draws of below(3) land 999,444 / 1,000,379 / 1,000,177.
• PCG32 reference outputs for pcg32_srandom_r(42, 54): 0xa15c02b7, 0x7b47f409, 0xba1d3330, 0x83d2f293, 0xbfa4784b, 0xcbed606e (the code above reproduces them).
• xoshiro256** from state [1, 2, 3, 4]: first outputs 11520, 0, 1509978240, 1215971899390074240.
• Lemire rejection happens with probability below n/2³² (or n/2⁶⁴); the division runs only on the slow path.
• Periods: PCG32 2⁶⁴ per stream (exhausting it takes years at 1 ns per output, but use at most √period ≈ 2³² outputs per stream when birthday-type statistics matter); xoshiro256** 2²⁵⁶ − 1.

Efficiency

• Cost: measured with the code above (Apple M4 Pro, -O): PCG32 1.0 ns per 32-bit output (the 64-bit multiply is a serial dependency), xoshiro256** 0.76 ns per 64-bit output (about 10 GB/s), SplitMix64 0.50 ns, below(n) 0.84 ns. State is 16 bytes (PCG32) or 32 bytes (xoshiro256**). Bound: the dependency chain on the state (latency), not memory.
• Levers:
  • Lemire multiply-shift for ranges — the division runs only on a rare slow path (Lemire 2019) — bounded integers, shuffles
  • top 53 bits times 2⁻⁵³ for f64 — one shift and one multiply, unbiased on the 2⁵³ grid — uniform floats
  • several independent generators in SIMD lanes — bulk generation without extra latency — large buffers
  • one jumped or separate-stream generator per thread — no locks, no shared state — parallel runs
  • counter-based Philox or Threefry — a stream per thread index without state (Salmon et al. 2011) — GPUs with millions of threads
• Naive → standard → state of the art:
  • Mersenne Twister · 19937-bit state, fails linear-complexity tests
  • PCG32 or xoshiro256** with Lemire ranges · 0.76-1.0 ns per output, 16-32 bytes of state (measured above)
  • same generators per SIMD lane or thread, or counter-based generators on GPU · throughput scales with lanes

Implementation notes

• Seeding: never seed xoshiro with small raw integers directly (state [1, 0, 0, 0] needs many outputs to look random); expand with SplitMix64. For PCG, the reference seeding above already mixes.
• Streams: PCG gives a stream per increment; xoshiro uses jump() (2¹²⁸ apart) per thread, or long_jump() (2¹⁹²) per process.
• Floats: (x >> 11) as f64 * 2⁻⁵³ is uniform on the 2⁵³ grid in [0, 1); dividing by u64::MAX can return 1.0 and is biased.
• Never use x % n for ranges with large n: bias up to n/2³²; and never (x as f64 * n as f64) as u64 for large n.
• Reproducibility across versions: pin the algorithm in your own code or crate version; library "default RNGs" change (Rust rand's StdRng is explicitly not portable across versions).
• Weak low bits: xoshiro256+ (one +, no ) has low bits of low linear complexity; use the top bits for floats. xoshiro256 has no such caveat.
• Vectorization: run several independent generators in SIMD lanes for bulk generation.

Pitfalls

• Seeding with the time in seconds for parallel jobs started together: identical streams.
• Sharing one generator across threads behind a mutex: slow; or without one: data race and duplicate values. Give each thread its own (jumped) generator.
• Using these for security tokens, password reset links or shuffles in gambling: predictable.
• Expecting identical results from two generators "with the same seed" in different libraries: seeding procedures differ.
• Calling jump() on a copy and forgetting to advance the original: two threads end up on the same stream.

Variants and alternatives

In production

NumPy default_rng (PCG64), Rust rand_pcg and rand_xoshiro (and SmallRng = Xoshiro256++ on 64-bit), Java 17 RandomGenerator (L64X128MixRandom, Xoshiro256PlusPlus), .NET 6+ Random (xoshiro256** for the parameterless constructor), Julia default RNG (Xoshiro256++), Lua 5.4 math.random (xoshiro256**), GNU Fortran random_number (xoshiro256**).

Implementations

Sources

• Melissa E. O'Neill, PCG: A Family of Simple Fast Space-Efficient Statistically Good Algorithms for Random Number Generation, Harvey Mudd College technical report HMC-CS-2014-0905, 2014. pdf
• David Blackman, Sebastiano Vigna, Scrambled Linear Pseudorandom Number Generators, ACM Transactions on Mathematical Software 47(4), 2021. paper · code
• Daniel Lemire, Fast Random Integer Generation in an Interval, ACM Transactions on Modeling and Computer Simulation 29(1), 2019. paper · arXiv
• Guy L. Steele Jr., Doug Lea, Christine H. Flood, Fast Splittable Pseudorandom Number Generators, OOPSLA 2014. paper
• Forrest B. Brown, Random Number Generation with Arbitrary Strides, Transactions of the American Nuclear Society 71, 1994.
• Pierre L'Ecuyer, Richard Simard, TestU01: A C Library for Empirical Testing of Random Number Generators, ACM Transactions on Mathematical Software 33(4), 2007. paper
• John K. Salmon et al., Parallel Random Numbers: As Easy as 1, 2, 3, SC 2011. paper