ENH: Add advance and jump to PCG64

Add advance and jump Clean other PRNGs
numpy · seberg · May 28, 2019 · Feb 20, 2018 · Feb 22, 2018 · Feb 23, 2018
commit 6b81fa0a1ebdd769dcd3ce82ef98dcbb0fa17767
diff --git a/_randomgen/TODO.md b/_randomgen/TODO.md
@@ -1,28 +1,28 @@
 # TODO
 2. Add dSFMT
-5. Augment state to have binomial structure
 6. Port over 0 parameter distributions
    * standard exponential ziggurat float   
    * standard normal ziggurat
    * standard normal ziggurat float
 7. Remove SplitMix64 as an external generator
 8. Restore ability to use `out` in core distributions
-10. Seed/Inc for PCG64
-11. Advance/Jump for PCG64
 12. Key/Counter for ThreeFry
+13. Simplify state
 
 ## Done
 1. Add PCG64
 3. Add xorshift2014
 4. Augment state to include has_gauss and gauss
+5. Augment state to have binomial structure
 6. Port over 0 parameter distributions
    * standard exponential ziggurat
    * standard exponential float
    * standard normal
    * standard normal float
    * standard gamma - Not implement: This is a 1 param
    * standard gamma float - Not implement: This is a 1 param
-
 9. Add correct carry for ThreeFry to allow full set of counters.  Important when implemeting jump
+10. Seed/Inc for PCG64
+11. Advance/Jump for PCG64
 0. NOT IMPLEMENTABLE due to limits on inheritance in Cython: Use inheritance to simplify CorePRNG structure. The natural base is 
    xoroshiro128.
diff --git a/_randomgen/core_prng/common.pxd b/_randomgen/core_prng/common.pxd
@@ -5,6 +5,27 @@ cdef extern from "src/distributions/distributions.h":
 
     ctypedef float (*random_float_0)(void *st) nogil
 
+    cdef struct s_binomial_t:
+        int has_binomial;
+        double psave;
+        long nsave;
+        double r;
+        double q;
+        double fm;
+        long m;
+        double p1;
+        double xm;
+        double xl;
+        double xr;
+        double c;
+        double laml;
+        double lamr;
+        double p2;
+        double p3;
+        double p4;
+
+    ctypedef s_binomial_t binomial_t
+
     cdef struct prng:
         void *state
         uint64_t (*next_uint64)(void *st)
@@ -14,6 +35,7 @@ cdef extern from "src/distributions/distributions.h":
         double gauss
         int has_gauss_f
         float gauss_f
+        binomial_t *binomial
 
     ctypedef prng prng_t
 

diff --git a/_randomgen/core_prng/generator.pyx b/_randomgen/core_prng/generator.pyx
@@ -2,6 +2,7 @@ import numpy as np
 cimport numpy as np
 from cpython.pycapsule cimport PyCapsule_IsValid, PyCapsule_GetPointer
 from common cimport *
+from libc.stdlib cimport malloc, free
 
 cimport numpy as np
 import numpy as np

diff --git a/_randomgen/core_prng/mt19937.pyx b/_randomgen/core_prng/mt19937.pyx
@@ -56,6 +56,7 @@ cdef class MT19937:
     def __init__(self, seed=None):
         self.rng_state = <mt19937_state *>malloc(sizeof(mt19937_state))
         self._prng = <prng_t *>malloc(sizeof(prng_t))
+        self._prng.binomial = <binomial_t *>malloc(sizeof(binomial_t))
         self.seed(seed)
 
         self._prng.state = <void *>self.rng_state
@@ -68,6 +69,7 @@ cdef class MT19937:
 
     def __dealloc__(self):
         free(self.rng_state)
+        free(self._prng.binomial)
         free(self._prng)
 
     # Pickling support:
@@ -82,10 +84,15 @@ cdef class MT19937:
                 (self.state['prng'],),
                 self.state)
 
-    def __random_integer(self):
+    def __random_integer(self, bits=64):
         """
         64-bit Random Integers from the PRNG
 
+        Parameters
+        ----------
+        bits : {32, 64}
+            Number of random bits to return
+
         Returns
         -------
         rv : int
@@ -95,7 +102,12 @@ cdef class MT19937:
         -----
         Testing only
         """
-        return mt19937_next64(self.rng_state)
+        if bits == 64:
+            return self._prng.next_uint64(self._prng.state)
+        elif bits == 32:
+            return self._prng.next_uint32(self._prng.state)
+        else:
+            raise ValueError('bits must be 32 or 64')
 
     def seed(self, seed=None):
         """

diff --git a/_randomgen/core_prng/pcg64.pyx b/_randomgen/core_prng/pcg64.pyx
@@ -33,6 +33,7 @@ cdef extern from "src/pcg64/pcg64.h":
     uint64_t pcg64_next64(pcg64_state *state)  nogil
     uint64_t pcg64_next32(pcg64_state *state)  nogil
     void pcg64_jump(pcg64_state  *state)
+    void pcg64_advance(pcg64_state *state, uint64_t *step)
 
 
 cdef uint64_t pcg64_uint64(void* st):# nogil:
@@ -62,20 +63,18 @@ cdef class PCG64:
     cdef prng_t *_prng
     cdef public object _prng_capsule
 
-    def __init__(self, seed=None):
+    def __init__(self, seed=None, inc=1):
         self.rng_state = <pcg64_state *>malloc(sizeof(pcg64_state))
         self.rng_state.pcg_state = <pcg64_random_t *>malloc(sizeof(pcg64_random_t))
         self._prng = <prng_t *>malloc(sizeof(prng_t))
-        self.seed(seed)
+        self._prng.binomial = <binomial_t *>malloc(sizeof(binomial_t))
+        self.seed(seed, inc)
 
         self._prng.state = <void *>self.rng_state
         self._prng.next_uint64 = &pcg64_uint64
         self._prng.next_uint32 = &pcg64_uint32
         self._prng.next_double = &pcg64_double
 
-        self.rng_state.pcg_state.inc.high = 0
-        self.rng_state.pcg_state.inc.low = 1
-
         cdef const char *name = "CorePRNG"
         self._prng_capsule = PyCapsule_New(<void *>self._prng, name, NULL)
 
@@ -93,6 +92,7 @@ cdef class PCG64:
 
     def __dealloc__(self):
         free(self.rng_state)
+        free(self._prng.binomial)
         free(self._prng)
 
     def _reset_state_variables(self):
@@ -124,7 +124,7 @@ cdef class PCG64:
         else:
             raise ValueError('bits must be 32 or 64')
 
-    def seed(self, seed=None):
+    def seed(self, seed=None, inc=1):
         """
         seed(seed=None, stream=None)
 
@@ -138,6 +138,8 @@ cdef class PCG64:
         ----------
         seed : int, optional
             Seed for ``RandomState``.
+        inc : int, optional
+            Increment to use for PCG stream
 
         Raises
         ------
@@ -156,6 +158,8 @@ cdef class PCG64:
             state = entropy.seed_by_array(seed, 2)
         self.rng_state.pcg_state.state.high = <uint64_t>int(state[0])
         self.rng_state.pcg_state.state.low = <uint64_t>int(state[1])
+        self.rng_state.pcg_state.inc.high = inc // 2**64
+        self.rng_state.pcg_state.inc.low = inc % 2**64
         self._reset_state_variables()
 
     @property
@@ -186,3 +190,13 @@ cdef class PCG64:
         self.rng_state.pcg_state.inc.low = value['state']['inc'] % 2 ** 64
         self.rng_state.has_uint32 = value['has_uint32']
         self.rng_state.uinteger = value['uinteger']
+
+    def advance(self, step):
+        cdef np.ndarray delta = np.empty(2,dtype=np.uint64)
+        delta[0] = step // 2**64
+        delta[1] = step % 2**64
+        pcg64_advance(self.rng_state, <uint64_t *>delta.data)
+        return self
+
+    def jump(self):
+        return self.advance(2**64)
diff --git a/_randomgen/core_prng/splitmix64.pyx b/_randomgen/core_prng/splitmix64.pyx
@@ -53,6 +53,7 @@ cdef class SplitMix64:
     def __init__(self, seed=None):
         self.rng_state = <splitmix64_state *>malloc(sizeof(splitmix64_state))
         self._prng = <prng_t *>malloc(sizeof(prng_t))
+        self._prng.binomial = <binomial_t *>malloc(sizeof(binomial_t))
         self.seed(seed)
 
         self._prng.state = <void *>self.rng_state
@@ -65,6 +66,7 @@ cdef class SplitMix64:
 
     def __dealloc__(self):
         free(self.rng_state)
+        free(self._prng.binomial)
         free(self._prng)
 
     # Pickling support:
@@ -84,10 +86,15 @@ cdef class SplitMix64:
         self.rng_state.has_uint32 = 0
         self.rng_state.uinteger = 0
 
-    def __random_integer(self):
+    def __random_integer(self, bits=64):
         """
         64-bit Random Integers from the PRNG
 
+        Parameters
+        ----------
+        bits : {32, 64}
+            Number of random bits to return
+
         Returns
         -------
         rv : int
@@ -97,7 +104,12 @@ cdef class SplitMix64:
         -----
         Testing only
         """
-        return splitmix64_next64(self.rng_state)
+        if bits == 64:
+            return self._prng.next_uint64(self._prng.state)
+        elif bits == 32:
+            return self._prng.next_uint32(self._prng.state)
+        else:
+            raise ValueError('bits must be 32 or 64')
 
     def seed(self, seed=None):
         """

diff --git a/_randomgen/core_prng/src/distributions/binomial.h b/_randomgen/core_prng/src/distributions/binomial.h
diff --git a/_randomgen/core_prng/src/distributions/distributions.h b/_randomgen/core_prng/src/distributions/distributions.h
@@ -20,6 +20,27 @@ typedef int bool;
 typedef double (*random_double_0)(void *st);
 typedef float (*random_float_0)(void *st);
 
+typedef struct s_binomial_t
+{
+    int has_binomial; /* !=0: following parameters initialized for binomial */
+    double psave;
+    long nsave;
+    double r;
+    double q;
+    double fm;
+    long m;
+    double p1;
+    double xm;
+    double xl;
+    double xr;
+    double c;
+    double laml;
+    double lamr;
+    double p2;
+    double p3;
+    double p4;
+} binomial_t;
+
 typedef struct prng {
   void *state;
   uint64_t (*next_uint64)(void *st);
@@ -29,6 +50,7 @@ typedef struct prng {
   double gauss;
   int has_gauss_f;
   float gauss_f;
+  binomial_t *binomial;
 } prng_t;
 
 float random_float(prng_t *prng_state);

diff --git a/_randomgen/core_prng/src/pcg64/pcg64.c b/_randomgen/core_prng/src/pcg64/pcg64.c
@@ -89,4 +89,17 @@ pcg128_t pcg_advance_lcg_128(pcg128_t state, pcg128_t delta, pcg128_t cur_mult,
 #endif
 
 extern inline uint64_t pcg64_next64(pcg64_state *state);
-extern inline uint32_t pcg64_next32(pcg64_state *state);
+extern inline uint32_t pcg64_next32(pcg64_state *state);
+
+#if __SIZEOF_INT128__ && !defined(PCG_FORCE_EMULATED_128BIT_MATH)
+extern void pcg64_advance(pcg64_state *state, pcg128_t step) {
+    pcg64_advance_r(state->pcg_state, step);
+}
+#else
+extern void pcg64_advance(pcg64_state *state, uint64_t *step) {
+    pcg128_t delta;
+    delta.high = step[0];
+    delta.low = step[1];
+    pcg64_advance_r(state->pcg_state, delta);
+}
+#endif
diff --git a/_randomgen/core_prng/src/pcg64/pcg64.h b/_randomgen/core_prng/src/pcg64/pcg64.h
@@ -234,4 +234,10 @@ static inline uint32_t pcg64_next32(pcg64_state *state) {
   state->has_uint32 = 1;
   state->uinteger = (uint32_t)(next & 0xffffffff);
   return (uint32_t)(next >> 32);
-}
+}
+
+#if __SIZEOF_INT128__ && !defined(PCG_FORCE_EMULATED_128BIT_MATH)
+void pcg64_advance(pcg64_state *state, pcg128_t step);
+#else
+void pcg64_advance(pcg64_state *state, uint64_t *step);
+#endif
diff --git a/_randomgen/core_prng/threefry.pyx b/_randomgen/core_prng/threefry.pyx
@@ -68,6 +68,7 @@ cdef class ThreeFry:
         self.rng_state.ctr = <threefry4x64_ctr_t *>malloc(sizeof(threefry4x64_ctr_t))
         self.rng_state.key = <threefry4x64_key_t *>malloc(sizeof(threefry4x64_key_t))
         self._prng = <prng_t *>malloc(sizeof(prng_t))
+        self._prng.binomial = <binomial_t *>malloc(sizeof(binomial_t))
         self.seed(seed)
 
         self._prng.state = <void *>self.rng_state
@@ -94,6 +95,7 @@ cdef class ThreeFry:
         free(self.rng_state.ctr)
         free(self.rng_state.key)
         free(self.rng_state)
+        free(self._prng.binomial)
         free(self._prng)
 
     def _reset_state_variables(self):
@@ -122,9 +124,9 @@ cdef class ThreeFry:
         Testing only
         """
         if bits == 64:
-            return threefry_next64(self.rng_state)
+            return self._prng.next_uint64(self._prng.state)
         elif bits == 32:
-            return threefry_next32(self.rng_state)
+            return self._prng.next_uint32(self._prng.state)
         else:
             raise ValueError('bits must be 32 or 64')
 
@@ -200,6 +202,7 @@ cdef class ThreeFry:
     def jump(self):
         """Jump the state as-if 2**128 draws have been made"""
         threefry_jump(self.rng_state)
+        return self
 
     def advance(self, step):
         """Advance the state as-if a specific number of draws have been made"""

diff --git a/_randomgen/core_prng/xoroshiro128.pyx b/_randomgen/core_prng/xoroshiro128.pyx
@@ -55,6 +55,7 @@ cdef class Xoroshiro128:
     def __init__(self, seed=None):
         self.rng_state = <xoroshiro128_state *>malloc(sizeof(xoroshiro128_state))
         self._prng = <prng_t *>malloc(sizeof(prng_t))
+        self._prng.binomial = <binomial_t *>malloc(sizeof(binomial_t))
         self.seed(seed)
 
         self._prng.state = <void *>self.rng_state
@@ -79,6 +80,7 @@ cdef class Xoroshiro128:
 
     def __dealloc__(self):
         free(self.rng_state)
+        free(self._prng.binomial)
         free(self._prng)
 
     def _reset_state_variables(self):
@@ -146,6 +148,7 @@ cdef class Xoroshiro128:
 
     def jump(self):
         xoroshiro128_jump(self.rng_state)
+        return self
 
     @property
     def state(self):