Merge pull request #101 from bashtage/ziggurat-exponential

ENH: Ziggurat method for exponential generator

Author: bashtage

ci/performance.py

@@ -25,7 +25,7 @@
 functions = {'randint': {'low': 2 ** 31, 'dtype': 'uint32'},
              'random_sample': {},
              'random_raw': {'output': False},
-             'standard_exponential': {},
+             'standard_exponential': {'method': 'zig'},
              'standard_gamma': {'shape': 2.4},
              'standard_normal': {'method': 'zig'},
              'multinomial': {'n': 20, 'pvals': [1.0 / 6.0] * np.ones(6)},
@@ -41,7 +41,7 @@ def timer(prng: str, fn: str, args: dict):
         # Differences with NumPy
         if fn in ('random_raw', 'complex_normal'):
             return np.nan
-        if fn == 'standard_normal':
+        if fn in ('standard_normal','standard_exponential'):
             args = {k: v for k, v in args.items() if k != 'method'}
     elif prng == 'mt19937' and fn == 'random_raw':  # To make comparable
         args['size'] = 2 * args['size']
@@ -125,4 +125,4 @@ def timer(prng: str, fn: str, args: dict):
 lines.insert(1, '    \n')
 lines.insert(1, '    :widths: 14,14,14,14,14,14,14,14\n')
 lines.insert(0, '.. csv-table::\n')
-print(''.join(lines))
+print(''.join(lines))

doc/source/change-log.rst

@@ -5,7 +5,22 @@ Change Log
 
 Since Version 1.13
 ------------------
-* Add SIMD-oriented Fast Mersenne Twister (`SFMT <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/>`) generator.
+* Add Ziggurat generator for standard exponential
+  (:meth:`~randomstate.prng.mt19937.standard_exponential`) for both floats and
+  doubles
+
+.. ipython:: python
+
+   import numpy as np
+   import randomstate as rs
+   rs.seed(23456)
+   rs.standard_exponential(3, method='zig') # New method
+
+   rs.standard_exponential(3, method='inv') # Old method
+
+
+* Add SIMD-oriented Fast Mersenne Twister
+  (`SFMT <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/>`_) generator.
 * Add complex normal (:meth:`~randomstate.prng.mt19937.complex_normal`)
 
 Version 1.13

doc/source/index.rst

@@ -10,14 +10,16 @@ What's New or Different
   ``/dev/urandom`` on Unix).
 * Simulate from the complex normal distribution
   (:meth:`~randomstate.prng.mt19937.complex_normal`)
-* The normal generator supports a 256-step Ziggurat method which is 2-6 times
-  faster than NumPy's ``standard_normal``.  This generator can be accessed
-  by passing the keyword argument ``method='zig'``.
+* The both the normal and exponential generators support 256-step Ziggurat
+  methods which are 2-6 times faster than NumPy's default implementation in
+  ``standard_normal`` or ``standard_exponential``.  This Ziggurat generator can
+  be accessed by passing the keyword argument ``method='zig'``.
 
 .. ipython:: python
 
-  from randomstate import standard_normal
+  from randomstate import standard_normal, standard_exponential
   standard_normal(100000, method='zig')
+  standard_exponential(100000, method='zig')
 
 * Optional ``dtype`` argument that accepts ``np.float32`` or ``np.float64``
   to produce either single or double prevision uniform random variables for
@@ -45,7 +47,7 @@ What's New or Different
   * Standard Gammas (:meth:`~randomstate.prng.mt19937.standard_gamma`)
   * Standard Exponentials (:meth:`~randomstate.prng.mt19937.standard_exponential`)
 
-  This allows mulththreading to fill large arrays in chunks using suitable
+  This allows multithreading to fill large arrays in chunks using suitable
   PRNGs in parallel.
 
 .. ipython:: python

papers/ziggurat.pdf

@@ -1557,3 +1557,99 @@ void random_bounded_bool_fill(aug_state *state, npy_bool off, npy_bool rng, npy_
         out[i] = (buf & 0x00000001) != 0;
     }
 }
+
+static inline double standard_exponential_zig_double(aug_state* state);
+
+static double standard_exponential_zig_double_unlikely(aug_state* state, uint8_t idx, double x)
+{
+    if (idx == 0)
+    {
+        return ziggurat_exp_r - log(random_double(state));
+    }
+    else if ((fe_double[idx - 1] - fe_double[idx]) * random_double(state) + fe_double[idx] < exp(-x))
+    {
+        return x;
+    }
+    else
+    {
+        return standard_exponential_zig_double(state);
+    }
+}
+
+static inline double standard_exponential_zig_double(aug_state* state)
+{
+    uint64_t ri;
+    uint8_t idx;
+    double x;
+    ri = random_uint64(state);
+    ri >>= 3;
+    idx = ri & 0xFF;
+    ri >>= 8;
+    x = ri * we_double[idx];
+    if (ri < ke_double[idx]){
+        return x;  // 98.9% of the time we return here 1st try
+    }
+    return standard_exponential_zig_double_unlikely(state, idx, x);
+}
+
+double random_standard_exponential_zig_double(aug_state* state)
+{
+    return standard_exponential_zig_double(state);
+}
+
+void random_standard_exponential_zig_double_fill(aug_state* state, npy_intp count, double *out)
+{
+    npy_intp i;
+    for (i=0; i < count; i++) {
+        out[i] = standard_exponential_zig_double(state);
+    }
+}
+
+static inline float standard_exponential_zig_float(aug_state* state);
+
+static double standard_exponential_zig_float_unlikely(aug_state* state, uint8_t idx, float x)
+{
+    if (idx == 0)
+    {
+        // TODO: Replace with float ?
+        return ziggurat_exp_r - logf(random_float(state));
+    }
+    else if ((fe_float[idx - 1] - fe_float[idx]) * random_float(state) + fe_float[idx] < expf(-x))
+    {
+        return x;
+    }
+    else
+    {
+        return standard_exponential_zig_float(state);
+    }
+}
+
+static inline float standard_exponential_zig_float(aug_state* state)
+{
+    uint32_t ri;
+    uint8_t idx;
+    float x;
+    ri = random_uint32(state);
+    ri >>= 1;
+    idx = ri & 0xFF;
+    ri >>= 8;
+    x = ri * we_float[idx];
+    if (ri < ke_float[idx]){
+        return x;  // 98.9% of the time we return here 1st try
+    }
+    return standard_exponential_zig_float_unlikely(state, idx, x);
+}
+
+
+float random_standard_exponential_zig_float(aug_state* state)
+{
+    return standard_exponential_zig_float(state);
+}
+
+void random_standard_exponential_zig_float_fill(aug_state* state, npy_intp count, float *out)
+{
+    npy_intp i;
+    for (i=0; i < count; i++) {
+        out[i] = standard_exponential_zig_float(state);
+    }
+}

randomstate/distributions.c

@@ -162,4 +162,12 @@ extern void random_gauss_zig_julia_fill(aug_state *state, npy_intp count, double
 
 extern void random_gauss_zig_float_fill(aug_state *state, npy_intp count, float *out);
 
-extern void random_gauss_zig_double_fill(aug_state *state, npy_intp count, double *out);
+extern void random_gauss_zig_double_fill(aug_state *state, npy_intp count, double *out);
+
+extern double random_standard_exponential_zig(aug_state* state);
+
+extern void random_standard_exponential_zig_double_fill(aug_state* state, npy_intp count, double *out);
+
+extern float random_standard_exponential_zig_float(aug_state* state);
+
+extern void random_standard_exponential_zig_float_fill(aug_state* state, npy_intp count, float *out);

randomstate/distributions.h

@@ -80,6 +80,7 @@ cdef extern from "distributions.h":
     cdef double random_gauss_zig(aug_state* state) nogil
     cdef double random_gauss_zig_julia(aug_state* state) nogil
     cdef double random_standard_exponential(aug_state* state) nogil
+    cdef double random_standard_exponential_ziggurat(aug_state* state) nogil
     cdef double random_standard_cauchy(aug_state* state) nogil
 
     cdef double random_exponential(aug_state *state, double scale) nogil
@@ -124,12 +125,14 @@ cdef extern from "distributions.h":
 
     cdef void random_gauss_zig_float_fill(aug_state *state, intptr_t count, float *out) nogil
     cdef void random_uniform_fill_float(aug_state *state, intptr_t cnt, double *out) nogil
+    cdef void random_standard_exponential_zig_float_fill(aug_state* state, intptr_t count, float *out) nogil
     cdef void random_standard_exponential_fill_float(aug_state* state, intptr_t count, float *out) nogil
     cdef void random_gauss_fill_float(aug_state* state, intptr_t count, float *out) nogil
 
     cdef void random_gauss_zig_double_fill(aug_state* state, intptr_t count, double *out) nogil
     cdef void random_uniform_fill_double(aug_state *state, intptr_t cnt, double *out) nogil
     cdef void random_standard_exponential_fill_double(aug_state* state, intptr_t count, double *out) nogil
+    cdef void random_standard_exponential_zig_double_fill(aug_state* state, intptr_t count, double *out) nogil
     cdef void random_gauss_fill(aug_state* state, intptr_t count, double *out) nogil
     cdef void random_gauss_zig_julia_fill(aug_state* state, intptr_t count, double *out) nogil
 
@@ -2014,9 +2017,9 @@ cdef class RandomState:
                     0.0, '', CONS_NONE,
                     None)
 
-    def standard_exponential(self, size=None, dtype=np.float64, out=None):
+    def standard_exponential(self, size=None, dtype=np.float64, method=u'inv', out=None):
         """
-        standard_exponential(size=None, dtype=np.float64, out=None)
+        standard_exponential(size=None, dtype=np.float64, method='inv', out=None)
 
         Draw samples from the standard exponential distribution.
 
@@ -2033,6 +2036,9 @@ cdef class RandomState:
             Desired dtype of the result. All dtypes are determined by their
             name, either 'float64' or 'float32'. The default value is
             'float64'.
+        method : str, optional
+            Either 'inv' or 'zig'. 'inv' uses the default inverse CDF method.
+            'zig' uses the much faster Ziggurat method of Marsaglia and Tsang.
         out : ndarray, optional
             Alternative output array in which to place the result. If size is not None,
             it must have the same shape as the provided size and must match the type of
@@ -2048,17 +2054,28 @@ cdef class RandomState:
         Output a 3x8000 array:
 
         >>> n = np.random.standard_exponential((3, 8000))
-
         """
+        if method != u'zig' and method != u'inv':
+            raise ValueError("method must be either 'bm' or 'zig'")
         key = np.dtype(dtype).name
         if key == 'float64':
-            return double_fill(&self.rng_state,
-                               &random_standard_exponential_fill_double,
-                               size, self.lock, out)
+            if method == 'zig':
+                return double_fill(&self.rng_state,
+                                   &random_standard_exponential_zig_double_fill,
+                                   size, self.lock, out)
+            else:
+                return double_fill(&self.rng_state,
+                                   &random_standard_exponential_fill_double,
+                                   size, self.lock, out)
         elif key == 'float32':
-            return float_fill(&self.rng_state,
-                              &random_standard_exponential_fill_float,
-                              size, self.lock, out)
+            if method == 'zig':
+                return float_fill(&self.rng_state,
+                                   &random_standard_exponential_zig_float_fill,
+                                   size, self.lock, out)
+            else:
+                return float_fill(&self.rng_state,
+                                  &random_standard_exponential_fill_float,
+                                  size, self.lock, out)
         else:
             raise TypeError('Unsupported dtype "%s" for standard_exponential'
                             % key)

randomstate/randomstate.pyx

@@ -3,7 +3,7 @@
 
 extern inline uint64_t xoroshiro128plus_next(xoroshiro128plus_state* state);
 
-static inline uint64_t rotl(const uint64_t x, int k);
+inline uint64_t rotl(const uint64_t x, int k);
 
 void xoroshiro128plus_jump(xoroshiro128plus_state* state) {
 	static const uint64_t JUMP[] = { 0xbeac0467eba5facb, 0xd86b048b86aa9922 };

randomstate/src/xoroshiro128plus/xoroshiro128plus.c

@@ -20,7 +20,7 @@ void xoroshiro128plus_seed_by_array(xoroshiro128plus_state* state, uint64_t *see
 
 void xoroshiro128plus_init_state(xoroshiro128plus_state* state, uint64_t seed, uint64_t inc);
 
-static inline uint64_t rotl(const uint64_t x, int k) {
+inline uint64_t rotl(const uint64_t x, int k) {
 	return (x << k) | (x >> (64 - k));
 }