rand_pcg/
pcg128cm.rs

1// Copyright 2018-2021 Developers of the Rand project.
2// Copyright 2017 Paul Dicker.
3// Copyright 2014-2017, 2019 Melissa O'Neill and PCG Project contributors
4//
5// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
6// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
7// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
8// option. This file may not be copied, modified, or distributed
9// except according to those terms.
10
11//! PCG random number generators
12
13// This is the cheap multiplier used by PCG for 128-bit state.
14const MULTIPLIER: u64 = 15750249268501108917;
15
16use core::fmt;
17use rand_core::{impls, le, RngCore, SeedableRng};
18#[cfg(feature = "serde")]
19use serde::{Deserialize, Serialize};
20
21/// A PCG random number generator (CM DXSM 128/64 (LCG) variant).
22///
23/// Permuted Congruential Generator with 128-bit state, internal Linear
24/// Congruential Generator, and 64-bit output via "double xorshift multiply"
25/// output function.
26///
27/// This is a 128-bit LCG with explicitly chosen stream with the PCG-DXSM
28/// output function. This corresponds to `pcg_engines::cm_setseq_dxsm_128_64`
29/// from pcg_cpp and `PCG64DXSM` from NumPy.
30///
31/// Despite the name, this implementation uses 32 bytes (256 bit) space
32/// comprising 128 bits of state and 128 bits stream selector. These are both
33/// set by `SeedableRng`, using a 256-bit seed.
34///
35/// Note that while two generators with different stream parameter may be
36/// closely correlated, this is [mitigated][upgrading-pcg64] by the DXSM output function.
37///
38/// [upgrading-pcg64]: https://numpy.org/doc/stable/reference/random/upgrading-pcg64.html
39#[derive(Clone, PartialEq, Eq)]
40#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
41pub struct Lcg128CmDxsm64 {
42    state: u128,
43    increment: u128,
44}
45
46/// [`Lcg128CmDxsm64`] is also known as `PCG64DXSM`.
47pub type Pcg64Dxsm = Lcg128CmDxsm64;
48
49impl Lcg128CmDxsm64 {
50    /// Multi-step advance functions (jump-ahead, jump-back)
51    ///
52    /// The method used here is based on Brown, "Random Number Generation
53    /// with Arbitrary Stride,", Transactions of the American Nuclear
54    /// Society (Nov. 1994).  The algorithm is very similar to fast
55    /// exponentiation.
56    ///
57    /// Even though delta is an unsigned integer, we can pass a
58    /// signed integer to go backwards, it just goes "the long way round".
59    ///
60    /// Using this function is equivalent to calling `next_64()` `delta`
61    /// number of times.
62    #[inline]
63    pub fn advance(&mut self, delta: u128) {
64        let mut acc_mult: u128 = 1;
65        let mut acc_plus: u128 = 0;
66        let mut cur_mult = MULTIPLIER as u128;
67        let mut cur_plus = self.increment;
68        let mut mdelta = delta;
69
70        while mdelta > 0 {
71            if (mdelta & 1) != 0 {
72                acc_mult = acc_mult.wrapping_mul(cur_mult);
73                acc_plus = acc_plus.wrapping_mul(cur_mult).wrapping_add(cur_plus);
74            }
75            cur_plus = cur_mult.wrapping_add(1).wrapping_mul(cur_plus);
76            cur_mult = cur_mult.wrapping_mul(cur_mult);
77            mdelta /= 2;
78        }
79        self.state = acc_mult.wrapping_mul(self.state).wrapping_add(acc_plus);
80    }
81
82    /// Construct an instance compatible with PCG seed and stream.
83    ///
84    /// Note that the highest bit of the `stream` parameter is discarded
85    /// to simplify upholding internal invariants.
86    ///
87    /// Note that while two generators with different stream parameter may be
88    /// closely correlated, this is [mitigated][upgrading-pcg64] by the DXSM output function.
89    ///
90    /// PCG specifies the following default values for both parameters:
91    ///
92    /// - `state = 0xcafef00dd15ea5e5`
93    /// - `stream = 0xa02bdbf7bb3c0a7ac28fa16a64abf96`
94    ///
95    /// [upgrading-pcg64]: https://numpy.org/doc/stable/reference/random/upgrading-pcg64.html
96    pub fn new(state: u128, stream: u128) -> Self {
97        // The increment must be odd, hence we discard one bit:
98        let increment = (stream << 1) | 1;
99        Self::from_state_incr(state, increment)
100    }
101
102    #[inline]
103    fn from_state_incr(state: u128, increment: u128) -> Self {
104        let mut pcg = Self { state, increment };
105        // Move away from initial value:
106        pcg.state = pcg.state.wrapping_add(pcg.increment);
107        pcg.step();
108        pcg
109    }
110
111    #[inline(always)]
112    fn step(&mut self) {
113        // prepare the LCG for the next round
114        self.state = self
115            .state
116            .wrapping_mul(MULTIPLIER as u128)
117            .wrapping_add(self.increment);
118    }
119}
120
121// Custom Debug implementation that does not expose the internal state
122impl fmt::Debug for Lcg128CmDxsm64 {
123    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
124        write!(f, "Lcg128CmDxsm64 {{}}")
125    }
126}
127
128impl SeedableRng for Lcg128CmDxsm64 {
129    type Seed = [u8; 32];
130
131    /// We use a single 255-bit seed to initialise the state and select a stream.
132    /// One `seed` bit (lowest bit of `seed[8]`) is ignored.
133    fn from_seed(seed: Self::Seed) -> Self {
134        let mut seed_u64 = [0u64; 4];
135        le::read_u64_into(&seed, &mut seed_u64);
136        let state = u128::from(seed_u64[0]) | (u128::from(seed_u64[1]) << 64);
137        let incr = u128::from(seed_u64[2]) | (u128::from(seed_u64[3]) << 64);
138
139        // The increment must be odd, hence we discard one bit:
140        Self::from_state_incr(state, incr | 1)
141    }
142}
143
144impl RngCore for Lcg128CmDxsm64 {
145    #[inline]
146    fn next_u32(&mut self) -> u32 {
147        self.next_u64() as u32
148    }
149
150    #[inline]
151    fn next_u64(&mut self) -> u64 {
152        let res = output_dxsm(self.state);
153        self.step();
154        res
155    }
156
157    #[inline]
158    fn fill_bytes(&mut self, dest: &mut [u8]) {
159        impls::fill_bytes_via_next(self, dest)
160    }
161}
162
163#[inline(always)]
164fn output_dxsm(state: u128) -> u64 {
165    // See https://github.com/imneme/pcg-cpp/blob/ffd522e7188bef30a00c74dc7eb9de5faff90092/include/pcg_random.hpp#L1016
166    // for a short discussion of the construction and its original implementation.
167    let mut hi = (state >> 64) as u64;
168    let mut lo = state as u64;
169
170    lo |= 1;
171    hi ^= hi >> 32;
172    hi = hi.wrapping_mul(MULTIPLIER);
173    hi ^= hi >> 48;
174    hi = hi.wrapping_mul(lo);
175
176    hi
177}
rand_pcg/pcg128cm.rs

rand_pcg/
pcg128cm.rs
