ENH: Add access to raw PRNG values

Sheppard, Kevin · Sheppard, Kevin · commit 479676a902d1 · 2016-02-17T17:07:25.000-05:00
Implement method that allows direct access to the values produced by the PRNG
diff --git a/README.md b/README.md
@@ -58,7 +58,10 @@ an additional `method` keyword argument which can be `bm` or `zig` where
 
 * `random_entropy` - Read from the system entropy provider, which is commonly 
 used in cryptographic applications
-* `random_uintegers` - unsigned integers `[0, 2**64-1]` 
+* `random_uintegers` - unsigned integers `[0, 2**64-1]`
+* `random_raw` - Direct access to the values produced by the underlying PRNG.  
+The range of the values returned depends on the specifics of the PRNG 
+implementation.
 * `jump` - Jumps RNGs that support it.  `jump` moves the state a great 
 distance. _Only available if supported by the RNG._
 * `advance` - Advanced the core RNG 'as-if' a number of draws were made, 
diff --git a/README.rst b/README.rst
@@ -66,6 +66,10 @@ New Functions
 -  ``random_entropy`` - Read from the system entropy provider, which is
    commonly used in cryptographic applications
 -  ``random_uintegers`` - unsigned integers ``[0, 2**64-1]``
+-  ``random_raw`` - Direct access to the values produced by the
+   underlying PRNG.
+   The range of the values returned depends on the specifics of the PRNG
+   implementation.
 -  ``jump`` - Jumps RNGs that support it. ``jump`` moves the state a
    great distance. *Only available if supported by the RNG.*
 -  ``advance`` - Advanced the core RNG 'as-if' a number of draws were
diff --git a/doc/source/dsfmt.rst b/doc/source/dsfmt.rst
@@ -37,6 +37,8 @@ Simple random data
    sample
    choice
    bytes
+   random_uintegers
+   random_raw
 
 Permutations
 ============
diff --git a/doc/source/mlfg.rst b/doc/source/mlfg.rst
@@ -31,6 +31,8 @@ Simple random data
    sample
    choice
    bytes
+   random_uintegers
+   random_raw
 
 Permutations
 ============
diff --git a/doc/source/mrg32k3a.rst b/doc/source/mrg32k3a.rst
@@ -38,6 +38,8 @@ Simple random data
    sample
    choice
    bytes
+   random_uintegers
+   random_raw
 
 Permutations
 ============
diff --git a/doc/source/mt19937.rst b/doc/source/mt19937.rst
@@ -32,6 +32,8 @@ Simple random data
    sample
    choice
    bytes
+   random_uintegers
+   random_raw
 
 Permutations
 ============
diff --git a/doc/source/pcg32.rst b/doc/source/pcg32.rst
@@ -38,6 +38,8 @@ Simple random data
    sample
    choice
    bytes
+   random_uintegers
+   random_raw
 
 Permutations
 ============
diff --git a/doc/source/pcg64.rst b/doc/source/pcg64.rst
@@ -38,6 +38,8 @@ Simple random data
    sample
    choice
    bytes
+   random_uintegers
+   random_raw
 
 Permutations
 ============
diff --git a/doc/source/xorshift1024.rst b/doc/source/xorshift1024.rst
@@ -38,6 +38,8 @@ Simple random data
    sample
    choice
    bytes
+   random_uintegers
+   random_raw
 
 Permutations
 ============
diff --git a/doc/source/xorshift128.rst b/doc/source/xorshift128.rst
@@ -38,6 +38,8 @@ Simple random data
    sample
    choice
    bytes
+   random_uintegers
+   random_raw
 
 Permutations
 ============
diff --git a/randomstate/interface/dSFMT/dSFMT-shim.c b/randomstate/interface/dSFMT/dSFMT-shim.c
@@ -7,6 +7,8 @@ extern inline uint64_t random_uint64(aug_state* state);
 
 extern inline double random_double(aug_state* state);
 
+extern inline uint64_t random_raw(aug_state* state);
+
 void reset_buffer(aug_state* state)
 {
     int i = 0;
diff --git a/randomstate/interface/dSFMT/dSFMT-shim.h b/randomstate/interface/dSFMT/dSFMT-shim.h
@@ -63,6 +63,14 @@ static inline double random_double(aug_state* state)
     // return dsfmt_genrand_close1_open2(state->rng) - 1.0;
 }
 
+
+static inline uint64_t random_raw(aug_state* state)
+{
+    double d = random_double_from_buffer(state);
+    return (uint64_t *)&d;
+}
+
+
 extern void entropy_init(aug_state* state);
 
 extern void set_seed_by_array(aug_state* state, uint32_t init_key[], int key_length);
diff --git a/randomstate/interface/mlfg-1279-861/mlfg-1279-861-shim.c b/randomstate/interface/mlfg-1279-861/mlfg-1279-861-shim.c
@@ -6,6 +6,8 @@ extern inline uint64_t random_uint64(aug_state* state);
 
 extern inline double random_double(aug_state* state);
 
+extern inline uint64_t random_raw(aug_state* state);
+
 void set_seed(aug_state* state, uint64_t val)
 {
     mlfg_seed(state->rng, val);
diff --git a/randomstate/interface/mlfg-1279-861/mlfg-1279-861-shim.h b/randomstate/interface/mlfg-1279-861/mlfg-1279-861-shim.h
@@ -32,6 +32,10 @@ inline uint64_t random_uint64(aug_state* state)
     return out;
 }
 
+inline uint64_t random_raw(aug_state* state)
+{
+    return  mlfg_next(state->rng) >> 1;
+}
 
 inline double random_double(aug_state* state)
 {
diff --git a/randomstate/interface/mrg32k3a/mrg32k3a-shim.c b/randomstate/interface/mrg32k3a/mrg32k3a-shim.c
@@ -6,6 +6,8 @@ extern inline uint64_t random_uint64(aug_state* state);
 
 extern inline double random_double(aug_state* state);
 
+extern inline uint64_t random_raw(aug_state* state);
+
 void entropy_init(aug_state* state)
 {
     uint32_t buf[6] = { 0 };
diff --git a/randomstate/interface/mrg32k3a/mrg32k3a-shim.h b/randomstate/interface/mrg32k3a/mrg32k3a-shim.h
@@ -31,6 +31,11 @@ inline uint64_t random_uint64(aug_state* state)
     return (((uint64_t) mrg32k3a_random(state->rng) << 32) | mrg32k3a_random(state->rng));
 }
 
+inline uint64_t random_raw(aug_state* state)
+{
+    return (uint64_t)random_uint32(state);
+}
+
 inline double random_double(aug_state* state)
 {
     int32_t a = random_uint32(state) >> 5, b = random_uint32(state) >> 6;
diff --git a/randomstate/interface/pcg-32/pcg-32-shim.c b/randomstate/interface/pcg-32/pcg-32-shim.c
@@ -4,10 +4,13 @@ extern inline uint32_t random_uint32(aug_state* state);
 
 extern inline uint64_t random_uint64(aug_state* state);
 
+extern inline double random_double(aug_state* state);
+
+extern inline uint64_t random_raw(aug_state* state);
+
 extern inline void set_seed(aug_state* state, uint64_t seed, uint64_t inc);
 
 extern inline void advance_state(aug_state* state, uint64_t delta);
 
 extern inline void entropy_init(aug_state* state);
 
-extern inline double random_double(aug_state* state);
diff --git a/randomstate/interface/pcg-32/pcg-32-shim.h b/randomstate/interface/pcg-32/pcg-32-shim.h
@@ -35,6 +35,11 @@ inline void advance_state(aug_state* state, uint64_t delta)
     pcg32_advance_r(state->rng, delta);
 }
 
+inline uint64_t random_raw(aug_state* state)
+{
+    return (uint64_t)random_uint32(state);
+}
+
 inline void entropy_init(aug_state* state)
 {
     uint64_t seeds[2];
diff --git a/randomstate/interface/pcg-64/pcg-64-shim.c b/randomstate/interface/pcg-64/pcg-64-shim.c
@@ -5,10 +5,13 @@ extern inline uint32_t random_uint32(aug_state* state);
 
 extern inline uint64_t random_uint64(aug_state* state);
 
+extern inline double random_double(aug_state* state);
+
+extern inline uint64_t random_raw(aug_state* state);
+
 extern inline void advance_state(aug_state* state, pcg128_t delta);
 
 extern inline void entropy_init(aug_state* state);
 
 extern inline void set_seed(aug_state* state, pcg128_t seed, pcg128_t inc);
 
-extern inline double random_double(aug_state* state);
diff --git a/randomstate/interface/pcg-64/pcg-64-shim.h b/randomstate/interface/pcg-64/pcg-64-shim.h
@@ -38,6 +38,11 @@ inline uint64_t random_uint64(aug_state* state)
     return pcg64_random_r(state->rng);
 }
 
+inline uint64_t random_raw(aug_state* state)
+{
+    return random_uint64(state);
+}
+
 inline void set_seed(aug_state* state, pcg128_t seed, pcg128_t inc)
 {
     pcg64_srandom_r(state->rng, seed, inc);
diff --git a/randomstate/interface/random-kit/random-kit-shim.c b/randomstate/interface/random-kit/random-kit-shim.c
@@ -4,6 +4,8 @@ extern inline uint32_t random_uint32(aug_state* state);
 
 extern inline uint64_t random_uint64(aug_state* state);
 
+extern inline uint64_t random_raw(aug_state* state);
+
 extern inline double random_double(aug_state* state);
 
 void set_seed_by_array(aug_state* state, uint32_t *init_key, int key_length)
diff --git a/randomstate/interface/random-kit/random-kit-shim.h b/randomstate/interface/random-kit/random-kit-shim.h
@@ -35,6 +35,11 @@ inline double random_double(aug_state* state)
     return (a * 67108864.0 + b) / 9007199254740992.0;
 }
 
+inline uint64_t random_raw(aug_state* state)
+{
+    return (uint64_t)random_uint32(state);
+}
+
 extern void entropy_init(aug_state* state);
 
 extern void set_seed_by_array(aug_state* state, uint32_t *init_key, int key_length);
diff --git a/randomstate/interface/xorshift1024/xorshift1024-shim.c b/randomstate/interface/xorshift1024/xorshift1024-shim.c
@@ -6,6 +6,8 @@ extern inline uint64_t random_uint64(aug_state* state);
 
 extern inline double random_double(aug_state* state);
 
+extern inline uint64_t random_raw(aug_state* state);
+
 void set_seed(aug_state* state, uint64_t seed)
 {
     xorshift1024_seed(state->rng, seed);
diff --git a/randomstate/interface/xorshift1024/xorshift1024-shim.h b/randomstate/interface/xorshift1024/xorshift1024-shim.h
@@ -39,6 +39,11 @@ inline uint64_t random_uint64(aug_state* state)
     return xorshift1024_next(state->rng);
 }
 
+inline uint64_t random_raw(aug_state* state)
+{
+    return random_uint64(state);
+}
+
 inline double random_double(aug_state* state)
 {
     uint64_t rn;
diff --git a/randomstate/interface/xorshift128/xorshift128-shim.c b/randomstate/interface/xorshift128/xorshift128-shim.c
@@ -6,7 +6,9 @@ extern inline uint64_t random_uint64(aug_state* state);
 
 extern inline double random_double(aug_state* state);
 
-inline void set_seed(aug_state* state, uint64_t seed)
+extern inline uint64_t random_raw(aug_state* state);
+
+void set_seed(aug_state* state, uint64_t seed)
 {
     xorshift128_seed(state->rng, seed);
 }
diff --git a/randomstate/interface/xorshift128/xorshift128-shim.h b/randomstate/interface/xorshift128/xorshift128-shim.h
@@ -39,6 +39,11 @@ inline uint64_t random_uint64(aug_state* state)
     return xorshift128_next(state->rng);
 }
 
+inline uint64_t random_raw(aug_state* state)
+{
+    return random_uint64(state);
+}
+
 inline double random_double(aug_state* state)
 {
     uint64_t rn;
diff --git a/randomstate/randomstate.pyx b/randomstate/randomstate.pyx
@@ -58,6 +58,7 @@ cdef extern from "distributions.h":
 
     cdef uint64_t random_uint64(aug_state* state) nogil
     cdef uint32_t random_uint32(aug_state* state) nogil
+    cdef uint64_t random_raw(aug_state* state) nogil
 
     cdef long random_positive_int(aug_state* state) nogil
     cdef unsigned long random_uint(aug_state* state) nogil
@@ -585,20 +586,64 @@ cdef class RandomState:
                 return random_uint64(&self.rng_state)
             n = compute_numel(size)
             randoms64 = np.empty(n, np.uint64)
-            for i in range(n):
-                randoms64[i] = random_uint64(&self.rng_state)
+            with self.lock, nogil:
+                for i in range(n):
+                    randoms64[i] = random_uint64(&self.rng_state)
             return np.asarray(randoms64).reshape(size)
         elif bits == 32:
             if size is None:
                 return random_uint32(&self.rng_state)
             n = compute_numel(size)
             randoms32 = np.empty(n, np.uint32)
-            for i in range(n):
-                randoms32[i] = random_uint32(&self.rng_state)
+            with self.lock, nogil:
+                for i in range(n):
+                    randoms32[i] = random_uint32(&self.rng_state)
             return np.asarray(randoms32).reshape(size)
         else:
             raise ValueError('Unknown value of bits.  Must be either 32 or 64.')
 
+    def random_raw(self, size=None):
+        """
+        random_raw(self, size=None)
+
+        Return randoms as generated by the underlying PRNG
+
+        Parameters
+        ----------
+        size : int or tuple of ints, optional
+            Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
+            ``m * n * k`` samples are drawn.  Default is None, in which case a
+            single value is returned.
+
+        Returns
+        -------
+        out : uint or ndarray
+            Drawn samples.
+
+        Notes
+        -----
+        This method directly exposes the the raw underlying pseudo-random
+        number generator. All values are returned as unsigned 64-bit
+        values irrespective of the number of bits produced by the PRNG.
+
+        See the class docstring for the number of bits returned.
+        """
+        cdef np.ndarray randoms
+        cdef uint64_t *randoms_data
+        cdef Py_ssize_t i, n
+
+        if size is None:
+            return random_raw(&self.rng_state)
+
+        randoms = <np.ndarray>np.empty(size, np.uint64)
+        randoms_data = <uint64_t*>np.PyArray_DATA(randoms)
+        n = np.PyArray_SIZE(randoms)
+
+        with self.lock, nogil:
+            for i in range(n):
+                randoms_data[i] = random_raw(&self.rng_state)
+        return randoms
+
     # Pickling support:
     def __getstate__(self):
         return self.get_state()
diff --git a/randomstate/tests/test_smoke.py b/randomstate/tests/test_smoke.py
@@ -120,6 +120,10 @@ def test_jump(self):
     def test_random_uintegers(self):
         assert_(len(self.rs.random_uintegers(10)) == 10)
 
+    def test_random_uintegers(self):
+        assert_(len(self.rs.random_raw(10)) == 10)
+        assert_(self.rs.random_raw((10,10)).shape == (10,10))
+
     def test_uniform(self):
         r = self.rs.uniform(-1.0, 0.0, size=10)
         assert_(len(r) == 10)
