Add Lcg128CmDxsm64 generator compatible with NumPy's PCG64DXSM

benches/generators.rs

@@ -21,7 +21,7 @@ use rand::prelude::*;
 use rand::rngs::adapter::ReseedingRng;
 use rand::rngs::{mock::StepRng, OsRng};
 use rand_chacha::{ChaCha12Rng, ChaCha20Core, ChaCha20Rng, ChaCha8Rng};
-use rand_pcg::{Pcg32, Pcg64, Pcg64Mcg};
+use rand_pcg::{Pcg32, Pcg64, Pcg64Mcg, Pcg64Dxsm};
 
 macro_rules! gen_bytes {
     ($fnn:ident, $gen:expr) => {
@@ -44,6 +44,7 @@ gen_bytes!(gen_bytes_step, StepRng::new(0, 1));
 gen_bytes!(gen_bytes_pcg32, Pcg32::from_entropy());
 gen_bytes!(gen_bytes_pcg64, Pcg64::from_entropy());
 gen_bytes!(gen_bytes_pcg64mcg, Pcg64Mcg::from_entropy());
+gen_bytes!(gen_bytes_pcg64dxsm, Pcg64Dxsm::from_entropy());
 gen_bytes!(gen_bytes_chacha8, ChaCha8Rng::from_entropy());
 gen_bytes!(gen_bytes_chacha12, ChaCha12Rng::from_entropy());
 gen_bytes!(gen_bytes_chacha20, ChaCha20Rng::from_entropy());
@@ -73,6 +74,7 @@ gen_uint!(gen_u32_step, u32, StepRng::new(0, 1));
 gen_uint!(gen_u32_pcg32, u32, Pcg32::from_entropy());
 gen_uint!(gen_u32_pcg64, u32, Pcg64::from_entropy());
 gen_uint!(gen_u32_pcg64mcg, u32, Pcg64Mcg::from_entropy());
+gen_uint!(gen_u32_pcg64dxsm, u32, Pcg64Dxsm::from_entropy());
 gen_uint!(gen_u32_chacha8, u32, ChaCha8Rng::from_entropy());
 gen_uint!(gen_u32_chacha12, u32, ChaCha12Rng::from_entropy());
 gen_uint!(gen_u32_chacha20, u32, ChaCha20Rng::from_entropy());
@@ -85,6 +87,7 @@ gen_uint!(gen_u64_step, u64, StepRng::new(0, 1));
 gen_uint!(gen_u64_pcg32, u64, Pcg32::from_entropy());
 gen_uint!(gen_u64_pcg64, u64, Pcg64::from_entropy());
 gen_uint!(gen_u64_pcg64mcg, u64, Pcg64Mcg::from_entropy());
+gen_uint!(gen_u64_pcg64dxsm, u64, Pcg64Dxsm::from_entropy());
 gen_uint!(gen_u64_chacha8, u64, ChaCha8Rng::from_entropy());
 gen_uint!(gen_u64_chacha12, u64, ChaCha12Rng::from_entropy());
 gen_uint!(gen_u64_chacha20, u64, ChaCha20Rng::from_entropy());
@@ -109,6 +112,7 @@ macro_rules! init_gen {
 init_gen!(init_pcg32, Pcg32);
 init_gen!(init_pcg64, Pcg64);
 init_gen!(init_pcg64mcg, Pcg64Mcg);
+init_gen!(init_pcg64dxsm, Pcg64Dxsm);
 init_gen!(init_chacha, ChaCha20Rng);
 
 const RESEEDING_BYTES_LEN: usize = 1024 * 1024;

rand_pcg/CHANGELOG.md

@@ -4,6 +4,9 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](http://keepachangelog.com/en/1.0.0/)
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [Unreleased]
+- Add `Lcg128CmDxsm64` generator compatible with NumPy's `PCG64DXSM` (#1202)
+
 ## [0.3.1] - 2021-06-15
 - Add `advance` methods to RNGs (#1111)
 - Document dependencies between streams (#1122)

rand_pcg/src/lib.rs

@@ -38,7 +38,9 @@
 #![no_std]
 
 mod pcg128;
+mod pcg128cm;
 mod pcg64;
 
 pub use self::pcg128::{Lcg128Xsl64, Mcg128Xsl64, Pcg64, Pcg64Mcg};
+pub use self::pcg128cm::{Lcg128CmDxsm64, Pcg64Dxsm};
 pub use self::pcg64::{Lcg64Xsh32, Pcg32};

rand_pcg/src/pcg128.rs

@@ -14,7 +14,7 @@
 const MULTIPLIER: u128 = 0x2360_ED05_1FC6_5DA4_4385_DF64_9FCC_F645;
 
 use core::fmt;
-use rand_core::{le, Error, RngCore, SeedableRng};
+use rand_core::{impls, le, Error, RngCore, SeedableRng};
 #[cfg(feature = "serde1")] use serde::{Deserialize, Serialize};
 
 /// A PCG random number generator (XSL RR 128/64 (LCG) variant).
@@ -77,6 +77,9 @@ impl Lcg128Xsl64 {
 
     /// Construct an instance compatible with PCG seed and stream.
     ///
+    /// Note that the highest bit of the `stream` parameter is discarded
+    /// to simplify upholding internal invariants.
+    ///
     /// Note that two generators with different stream parameters may be closely
     /// correlated.
     ///
@@ -116,11 +119,11 @@ impl fmt::Debug for Lcg128Xsl64 {
     }
 }
 
-/// We use a single 255-bit seed to initialise the state and select a stream.
-/// One `seed` bit (lowest bit of `seed[8]`) is ignored.
 impl SeedableRng for Lcg128Xsl64 {
     type Seed = [u8; 32];
 
+    /// We use a single 255-bit seed to initialise the state and select a stream.
+    /// One `seed` bit (lowest bit of `seed[8]`) is ignored.
     fn from_seed(seed: Self::Seed) -> Self {
         let mut seed_u64 = [0u64; 4];
         le::read_u64_into(&seed, &mut seed_u64);
@@ -146,7 +149,7 @@ impl RngCore for Lcg128Xsl64 {
 
     #[inline]
     fn fill_bytes(&mut self, dest: &mut [u8]) {
-        fill_bytes_impl(self, dest)
+        impls::fill_bytes_via_next(self, dest)
     }
 
     #[inline]
@@ -237,8 +240,7 @@ impl SeedableRng for Mcg128Xsl64 {
         // Read as if a little-endian u128 value:
         let mut seed_u64 = [0u64; 2];
         le::read_u64_into(&seed, &mut seed_u64);
-        let state = u128::from(seed_u64[0])  |
-                    u128::from(seed_u64[1]) << 64;
+        let state = u128::from(seed_u64[0]) | u128::from(seed_u64[1]) << 64;
         Mcg128Xsl64::new(state)
     }
 }
@@ -257,7 +259,7 @@ impl RngCore for Mcg128Xsl64 {
 
     #[inline]
     fn fill_bytes(&mut self, dest: &mut [u8]) {
-        fill_bytes_impl(self, dest)
+        impls::fill_bytes_via_next(self, dest)
     }
 
     #[inline]
@@ -278,19 +280,3 @@ fn output_xsl_rr(state: u128) -> u64 {
     let xsl = ((state >> XSHIFT) as u64) ^ (state as u64);
     xsl.rotate_right(rot)
 }
-
-#[inline(always)]
-fn fill_bytes_impl<R: RngCore + ?Sized>(rng: &mut R, dest: &mut [u8]) {
-    let mut left = dest;
-    while left.len() >= 8 {
-        let (l, r) = { left }.split_at_mut(8);
-        left = r;
-        let chunk: [u8; 8] = rng.next_u64().to_le_bytes();
-        l.copy_from_slice(&chunk);
-    }
-    let n = left.len();
-    if n > 0 {
-        let chunk: [u8; 8] = rng.next_u64().to_le_bytes();
-        left.copy_from_slice(&chunk[..n]);
-    }
-}

rand_pcg/src/pcg128cm.rs

@@ -0,0 +1,182 @@
+// Copyright 2018-2021 Developers of the Rand project.
+// Copyright 2017 Paul Dicker.
+// Copyright 2014-2017, 2019 Melissa O'Neill and PCG Project contributors
+//
+// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
+// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
+// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
+// option. This file may not be copied, modified, or distributed
+// except according to those terms.
+
+//! PCG random number generators
+
+// This is the cheap multiplier used by PCG for 128-bit state.
+const MULTIPLIER: u64 = 15750249268501108917;
+
+use core::fmt;
+use rand_core::{impls, le, Error, RngCore, SeedableRng};
+#[cfg(feature = "serde1")] use serde::{Deserialize, Serialize};
+
+/// A PCG random number generator (CM DXSM 128/64 (LCG) variant).
+///
+/// Permuted Congruential Generator with 128-bit state, internal Linear
+/// Congruential Generator, and 64-bit output via "double xorshift multiply"
+/// output function.
+///
+/// This is a 128-bit LCG with explicitly chosen stream with the PCG-DXSM
+/// output function. This corresponds to `pcg_engines::cm_setseq_dxsm_128_64`
+/// from pcg_cpp and `PCG64DXSM` from NumPy.
+///
+/// Despite the name, this implementation uses 32 bytes (256 bit) space
+/// comprising 128 bits of state and 128 bits stream selector. These are both
+/// set by `SeedableRng`, using a 256-bit seed.
+///
+/// Note that while two generators with different stream parameter may be
+/// closely correlated, this is [mitigated][upgrading-pcg64] by the DXSM output function.
+///
+/// [upgrading-pcg64]: https://numpy.org/doc/stable/reference/random/upgrading-pcg64.html
+#[derive(Clone, PartialEq, Eq)]
+#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
+pub struct Lcg128CmDxsm64 {
+    state: u128,
+    increment: u128,
+}
+
+/// [`Lcg128CmDxsm64`] is also known as `PCG64DXSM`.
+pub type Pcg64Dxsm = Lcg128CmDxsm64;
+
+impl Lcg128CmDxsm64 {
+    /// Multi-step advance functions (jump-ahead, jump-back)
+    ///
+    /// The method used here is based on Brown, "Random Number Generation
+    /// with Arbitrary Stride,", Transactions of the American Nuclear
+    /// Society (Nov. 1994).  The algorithm is very similar to fast
+    /// exponentiation.
+    ///
+    /// Even though delta is an unsigned integer, we can pass a
+    /// signed integer to go backwards, it just goes "the long way round".
+    ///
+    /// Using this function is equivalent to calling `next_64()` `delta`
+    /// number of times.
+    #[inline]
+    pub fn advance(&mut self, delta: u128) {
+        let mut acc_mult: u128 = 1;
+        let mut acc_plus: u128 = 0;
+        let mut cur_mult = MULTIPLIER as u128;
+        let mut cur_plus = self.increment;
+        let mut mdelta = delta;
+
+        while mdelta > 0 {
+            if (mdelta & 1) != 0 {
+                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);
+            mdelta /= 2;
+        }
+        self.state = acc_mult.wrapping_mul(self.state).wrapping_add(acc_plus);
+    }
+
+    /// Construct an instance compatible with PCG seed and stream.
+    ///
+    /// Note that the highest bit of the `stream` parameter is discarded
+    /// to simplify upholding internal invariants.
+    ///
+    /// Note that while two generators with different stream parameter may be
+    /// closely correlated, this is [mitigated][upgrading-pcg64] by the DXSM output function.
+    ///
+    /// PCG specifies the following default values for both parameters:
+    ///
+    /// - `state = 0xcafef00dd15ea5e5`
+    /// - `stream = 0xa02bdbf7bb3c0a7ac28fa16a64abf96`
+    ///
+    /// [upgrading-pcg64]: https://numpy.org/doc/stable/reference/random/upgrading-pcg64.html
+    pub fn new(state: u128, stream: u128) -> Self {
+        // The increment must be odd, hence we discard one bit:
+        let increment = (stream << 1) | 1;
+        Self::from_state_incr(state, increment)
+    }
+
+    #[inline]
+    fn from_state_incr(state: u128, increment: u128) -> Self {
+        let mut pcg = Self { state, increment };
+        // Move away from initial value:
+        pcg.state = pcg.state.wrapping_add(pcg.increment);
+        pcg.step();
+        pcg
+    }
+
+    #[inline(always)]
+    fn step(&mut self) {
+        // prepare the LCG for the next round
+        self.state = self
+            .state
+            .wrapping_mul(MULTIPLIER as u128)
+            .wrapping_add(self.increment);
+    }
+}
+
+// Custom Debug implementation that does not expose the internal state
+impl fmt::Debug for Lcg128CmDxsm64 {
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+        write!(f, "Lcg128CmDxsm64 {{}}")
+    }
+}
+
+impl SeedableRng for Lcg128CmDxsm64 {
+    type Seed = [u8; 32];
+
+    /// We use a single 255-bit seed to initialise the state and select a stream.
+    /// One `seed` bit (lowest bit of `seed[8]`) is ignored.
+    fn from_seed(seed: Self::Seed) -> Self {
+        let mut seed_u64 = [0u64; 4];
+        le::read_u64_into(&seed, &mut seed_u64);
+        let state = u128::from(seed_u64[0]) | (u128::from(seed_u64[1]) << 64);
+        let incr = u128::from(seed_u64[2]) | (u128::from(seed_u64[3]) << 64);
+
+        // The increment must be odd, hence we discard one bit:
+        Self::from_state_incr(state, incr | 1)
+    }
+}
+
+impl RngCore for Lcg128CmDxsm64 {
+    #[inline]
+    fn next_u32(&mut self) -> u32 {
+        self.next_u64() as u32
+    }
+
+    #[inline]
+    fn next_u64(&mut self) -> u64 {
+        let val = output_dxsm(self.state);
+        self.step();
+        val
+    }
+
+    #[inline]
+    fn fill_bytes(&mut self, dest: &mut [u8]) {
+        impls::fill_bytes_via_next(self, dest)
+    }
+
+    #[inline]
+    fn try_fill_bytes(&mut self, dest: &mut [u8]) -> Result<(), Error> {
+        self.fill_bytes(dest);
+        Ok(())
+    }
+}
+
+#[inline(always)]
+fn output_dxsm(state: u128) -> u64 {
+    // See https://github.com/imneme/pcg-cpp/blob/ffd522e7188bef30a00c74dc7eb9de5faff90092/include/pcg_random.hpp#L1016
+    // for a short discussion of the construction and its original implementation.
+    let mut hi = (state >> 64) as u64;
+    let mut lo = state as u64;
+
+    lo |= 1;
+    hi ^= hi >> 32;
+    hi = hi.wrapping_mul(MULTIPLIER);
+    hi ^= hi >> 48;
+    hi = hi.wrapping_mul(lo);
+
+    hi
+}

rand_pcg/src/pcg64.rs

@@ -77,6 +77,9 @@ impl Lcg64Xsh32 {
 
     /// Construct an instance compatible with PCG seed and stream.
     ///
+    /// Note that the highest bit of the `stream` parameter is discarded
+    /// to simplify upholding internal invariants.
+    ///
     /// Note that two generators with different stream parameters may be closely
     /// correlated.
     ///
@@ -117,11 +120,11 @@ impl fmt::Debug for Lcg64Xsh32 {
     }
 }
 
-/// We use a single 127-bit seed to initialise the state and select a stream.
-/// One `seed` bit (lowest bit of `seed[8]`) is ignored.
 impl SeedableRng for Lcg64Xsh32 {
     type Seed = [u8; 16];
 
+    /// We use a single 127-bit seed to initialise the state and select a stream.
+    /// One `seed` bit (lowest bit of `seed[8]`) is ignored.
     fn from_seed(seed: Self::Seed) -> Self {
         let mut seed_u64 = [0u64; 2];
         le::read_u64_into(&seed, &mut seed_u64);