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