Merge pull request #1111 from markwkm/advance

Add an advance function to PCG generators

rand_pcg/src/pcg128.rs

@@ -40,6 +40,38 @@ pub struct Lcg128Xsl64 {
 pub type Pcg64 = Lcg128Xsl64;
 
 impl Lcg128Xsl64 {
+    /// 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;
+        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 PCG specifies default values for both parameters:
@@ -140,6 +172,38 @@ pub struct Mcg128Xsl64 {
 pub type Pcg64Mcg = Mcg128Xsl64;
 
 impl Mcg128Xsl64 {
+    /// 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;
+        let mut cur_plus: u128 = 0;
+        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.
     ///
     /// Note that PCG specifies a default value for the parameter:

rand_pcg/src/pcg64.rs

@@ -40,6 +40,38 @@ pub struct Lcg64Xsh32 {
 pub type Pcg32 = Lcg64Xsh32;
 
 impl Lcg64Xsh32 {
+    /// 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_32()` `delta`
+    /// number of times.
+    #[inline]
+    pub fn advance(&mut self, delta: u64) {
+        let mut acc_mult: u64 = 1;
+        let mut acc_plus: u64 = 0;
+        let mut cur_mult = MULTIPLIER;
+        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 PCG specifies default values for both parameters:

rand_pcg/tests/lcg128xsl64.rs

@@ -1,6 +1,19 @@
 use rand_core::{RngCore, SeedableRng};
 use rand_pcg::{Lcg128Xsl64, Pcg64};
 
+#[test]
+fn test_lcg128xsl64_advancing() {
+    for seed in 0..20 {
+        let mut rng1 = Lcg128Xsl64::seed_from_u64(seed);
+        let mut rng2 = rng1.clone();
+        for _ in 0..20 {
+            rng1.next_u64();
+        }
+        rng2.advance(20);
+        assert_eq!(rng1, rng2);
+    }
+}
+
 #[test]
 fn test_lcg128xsl64_construction() {
     // Test that various construction techniques produce a working RNG.

rand_pcg/tests/lcg64xsh32.rs

@@ -1,6 +1,19 @@
 use rand_core::{RngCore, SeedableRng};
 use rand_pcg::{Lcg64Xsh32, Pcg32};
 
+#[test]
+fn test_lcg64xsh32_advancing() {
+    for seed in 0..20 {
+        let mut rng1 = Lcg64Xsh32::seed_from_u64(seed);
+        let mut rng2 = rng1.clone();
+        for _ in 0..20 {
+            rng1.next_u32();
+        }
+        rng2.advance(20);
+        assert_eq!(rng1, rng2);
+    }
+}
+
 #[test]
 fn test_lcg64xsh32_construction() {
     // Test that various construction techniques produce a working RNG.

rand_pcg/tests/mcg128xsl64.rs

@@ -1,6 +1,19 @@
 use rand_core::{RngCore, SeedableRng};
 use rand_pcg::{Mcg128Xsl64, Pcg64Mcg};
 
+#[test]
+fn test_mcg128xsl64_advancing() {
+    for seed in 0..20 {
+        let mut rng1 = Mcg128Xsl64::seed_from_u64(seed);
+        let mut rng2 = rng1.clone();
+        for _ in 0..20 {
+            rng1.next_u64();
+        }
+        rng2.advance(20);
+        assert_eq!(rng1, rng2);
+    }
+}
+
 #[test]
 fn test_mcg128xsl64_construction() {
     // Test that various construction techniques produce a working RNG.