Merge pull request #59 from bashtage/ziggurat-float

ENH: Initial work on 32-bit ziggurat

Author: bashtage

README.md

@@ -3,8 +3,8 @@
 [![Travis Build Status](https://travis-ci.org/bashtage/ng-numpy-randomstate.svg?branch=master)](https://travis-ci.org/bashtage/ng-numpy-randomstate) 
 [![Appveyor Build Status](https://ci.appveyor.com/api/projects/status/odc5c4ukhru5xicl/branch/master?svg=true)](https://ci.appveyor.com/project/bashtage/ng-numpy-randomstate/branch/master)
 
-This is a library and generic interface for alternative random generators 
-in Python and Numpy. 
+This is a library and generic interface for alternative random 
+generators in Python and Numpy. 
 
 ## Features
 
@@ -30,22 +30,19 @@ import randomstate as rnd
 w = rnd.standard_normal(10000, method='zig')
 ```
 
-* Preliminary support for 32-bit floating randoms for core generators. 
-  Currently only uniforms (`random_sample`), exponentials 
-  (`standard_exponential`) and normals (`standard_normal` but only 
-  using Box-Muller, so `method='bm'` is required) have been implemented. 
-  Ultimately support should be avialable for:
-  
-    * Uniforms
-    * Exponentials
-    * Standard Gammas (via `standard_gamma`)
-    * Normals (currently only implemented using Box-Muller transformation)
+* Support for 32-bit floating randoms for core generators. 
+  Currently supported:
+
+    * Uniforms (`random_sample`)
+    * Exponentials (`standard_exponential`) 
+    * Normals (`standard_normal`, both Box-Muller and Ziggurat)
+    * Standard Gammas (via `standard_gamma`) 
 
   **WARNING**: The 32-bit generators are **experimental** and subjust 
   to change.
 
-  **Note**: There are no plans to extend the alternative precision generation to 
-  all random number types.
+  **Note**: There are _no_ plans to extend the alternative precision 
+  generation to all random number types.
 
 
 
@@ -59,11 +56,14 @@ The RNGs include:
 
 * [MT19937](https://github.com/numpy/numpy/blob/master/numpy/random/mtrand/),
  the NumPy rng
-* [dSFMT](http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/) a SSE2-aware 
-version of the MT19937 generator that is especially fast at generating doubles
-* [xorshift128+](http://xorshift.di.unimi.it/), [xoroshiro128+](http://xoroshiro.di.unimi.it/)
-[xorshift1024*](http://xorshift.di.unimi.it/)
-* [PCG32](http://www.pcg-random.org/) and [PCG64](http:w//www.pcg-random.org/)
+* [dSFMT](http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/) a 
+  SSE2-aware version of the MT19937 generator that is especially fast at 
+  generating doubles
+* [xorshift128+](http://xorshift.di.unimi.it/), 
+  [xoroshiro128+](http://xoroshiro.di.unimi.it/) and
+  [xorshift1024*](http://xorshift.di.unimi.it/)
+* [PCG32](http://www.pcg-random.org/) and 
+  [PCG64](http:w//www.pcg-random.org/)
 * [MRG32K3A](http://simul.iro.umontreal.ca/rng)
 * A multiplicative lagged fibonacci generator (LFG(63, 1279, 861, *))
 
@@ -75,8 +75,8 @@ support an additional `method` keyword argument which can be `bm` or
 `zig` where `bm` corresponds to the current method using the Box-Muller 
 transformation and `zig` uses the much faster (100%+) ziggurat method.
 * `random_sample` can produce either single precision (`np.float32`) or
-double precision (`np.float64`, the default) using an the optional keyword
-argument `dtype`.
+double precision (`np.float64`, the default) using an the optional 
+keyword argument `dtype`.
 
 ### New Functions
 
@@ -93,9 +93,9 @@ the RNG._
 
 ## Status
 
-* Complete drop-in replacement for `numpy.random.RandomState`. The `mt19937` 
-generator is identical to `numpy.random.RandomState`, and will produce an 
-identical sequence of random numbers for a given seed.   
+* Complete drop-in replacement for `numpy.random.RandomState`. The 
+`mt19937` generator is identical to `numpy.random.RandomState`, and 
+will produce an identical sequence of random numbers for a given seed.   
 * Builds and passes all tests on:
   * Linux 32/64 bit, Python 2.7, 3.4, 3.5 (should work on 2.6 and 3.3)
   * PC-BSD (FreeBSD) 64-bit, Python 2.7
@@ -133,13 +133,16 @@ Testing requires nose (1.3+).
 **Note:** it might work with other versions but only tested with these 
 versions. 
 
+## Development and Testing
+
 All development has been on 64-bit Linux, and it is regularly tested on 
-Travis-CI. The library is occasionally tested on Linux 32-bit,  OSX 10.10, 
-PC-BSD 10.2 (should also work on Free BSD) and Windows (Python 2.7/3.5, 
-both 32 and 64-bit).
+Travis-CI. The library is occasionally tested on Linux 32-bit,  
+OSX 10.10, PC-BSD 10.2 (should also work on Free BSD) and Windows 
+(Python 2.7/3.5, both 32 and 64-bit).
 
-Basic tests are in place for all RNGs. The MT19937 is tested against NumPy's 
-implementation for identical results. It also passes NumPy's test suite.
+Basic tests are in place for all RNGs. The MT19937 is tested against 
+NumPy's implementation for identical results. It also passes NumPy's 
+test suite.
 
 ## Installing
 
@@ -179,8 +182,8 @@ rs = randomstate.prng.mt19937.RandomState()
 rs.random_sample(100)
 ```
 
-Like NumPy, `randomstate` also exposes a single instance of the `mt19937` 
-generator directly at the module level so that commands like
+Like NumPy, `randomstate` also exposes a single instance of the 
+`mt19937` generator directly at the module level so that commands like
 
 ```python
 import randomstate
@@ -194,7 +197,8 @@ will work.
 Standard NCSA, plus sub licenses for components.
 
 ## Performance
-Performance is promising, and even the mt19937 seems to be faster than NumPy's mt19937. 
+Performance is promising, and even the mt19937 seems to be faster than 
+NumPy's mt19937. 
 
 ```
 Speed-up relative to NumPy (Uniform Doubles)

doc/make.bat

@@ -44,6 +44,7 @@ if "%1" == "help" (
 if "%1" == "clean" (
 	for /d %%i in (%BUILDDIR%\*) do rmdir /q /s %%i
 	del /q /s %BUILDDIR%\*
+	rmdir /q /s source\generated
 	goto end
 )
 

doc/source/change-log.rst

@@ -2,6 +2,20 @@
 
 Change Log
 ==========
+
+Head (not released)
+-------------------
+* Extended 32-bit generation to
+
+  * Uniforms (:meth:`~randomstate.entropy.random_sample` and :meth:`~randomstate.entropy.rand`)
+  * Normals (:meth:`~randomstate.entropy.standard_normal` and :meth:`~randomstate.entropy.randn`)
+  * Standard Gammas (:meth:`~randomstate.entropy.standard_gamma`)
+  * Standard Exponentials (:meth:`~randomstate.entropy.standard_exponential`)
+
+  using the ``dtype`` keyword.
+* Removed ``random_uintegers`` since these are special cases of
+  :meth:`~randomstate.entropy.randint`
+
 Version 1.11.2
 --------------
 * Added keyword argument `dtype` to `random_sample` which allows for single

doc/source/index.rst

@@ -11,16 +11,34 @@ What's New or Different
 * 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'``.
-* ``random_sample`` accepts the optional keyword argument ``dtype`` which
-  accepts ``np.float32`` or ``np.float64`` to produce either single or
-  double prevision uniform random variables
-* For changes since the previous release, see the :ref:`change-log`
 
-.. code-block:: python
+.. ipython:: python
 
   from randomstate import standard_normal
   standard_normal(100000, method='zig')
 
+* Optional ``dtype`` argument that accepts ``np.float32`` or ``np.float64``
+  to produce either single or double prevision uniform random variables for
+  select distributions
+
+    * Uniforms (:meth:`~randomstate.entropy.random_sample` and :meth:`~randomstate.entropy.rand`)
+    * Normals (:meth:`~randomstate.entropy.standard_normal` and :meth:`~randomstate.entropy.randn`)
+    * Standard Gammas (:meth:`~randomstate.entropy.standard_gamma`)
+    * Standard Exponentials (:meth:`~randomstate.entropy.standard_exponential`)
+
+.. ipython:: python
+
+  import numpy as np
+  import randomstate as rs
+  rs.seed(0)
+  rs.random_sample(3, dtype=np.float64)
+  rs.seed(0)
+  rs.random_sample(3, dtype=np.float32)
+
+
+
+* For changes since the previous release, see the :ref:`change-log`
+
 Parallel Generation
 -------------------
 

papers/generate_consts.py

@@ -0,0 +1,83 @@
+import numpy as np
+from numpy import log, sqrt, exp
+
+pdf = lambda x: exp(-0.5 * x * x)
+
+def inverse_pdf(p):
+    return sqrt(-2 * log(p))
+
+
+C = 256
+r = 3.6541528853610087963519472518
+Vr = 0.0049286732339746519
+x = np.zeros(C)
+x[0] = r
+for i in range(1, C - 1):
+    x[i] = inverse_pdf(pdf(x[i - 1]) + Vr / x[i - 1])
+
+x = x[::-1]
+
+# ki
+ki_final = {}
+scales = (2 ** 52, 2 ** 23)
+names = ('double', 'float')
+for name, scale in zip(names, scales):
+    ki = np.zeros(C, np.uint64)
+    prob_ratio = x[:-1] / x[1:]
+    ki[1:] = np.round(prob_ratio * scale).astype(np.uint64)
+    ki[0] = np.uint64(np.round(((x.max() * pdf(x.max())) / Vr) * scale))
+    digits = 10 if name == 'float' else 18
+    out = ["{0:#0{1}X}".format((ki[i]), digits) for i in range(C)]
+    ki_final[name] = out
+
+# wi
+wi_final = {}
+scales = (2 ** 52, 2 ** 23)
+names = ('double', 'float')
+for name, scale in zip(names, scales):
+    wi = x.copy()
+    wi[0] = Vr / pdf(x.max())
+    wi_final[name] = wi / scale
+
+# fi
+fi_final = {'double': pdf(x), 'float': pdf(x)}
+
+constants = {'ziggurat_nor_r': 3.6541528853610087963519472518,
+             'ziggurat_nor_inv_r': 0.27366123732975827203338247596,
+             'ziggurat_exp_r': 7.6971174701310497140446280481}
+
+type_map = {'uint64': 'uint64_t', 'uint32': 'uint32_t', 'double': 'double',
+            'float': 'float'}
+extra_text = {'uint64': 'ULL', 'uint32': 'UL', 'double': '', 'float': 'f'}
+
+
+def write(a, name, dtype):
+    ctype = type_map[dtype]
+    out = 'static const ' + ctype + ' ' + name + '[] = { \n'
+    format_str = '{0: .20e}' if dtype in ('double', 'float') else '{0}'
+    formatted = [format_str.format(a[i]) + extra_text[dtype] for i in range(len(a))]
+    lines = len(formatted) // 4
+    for i in range(lines):
+        temp = '    ' + ', '.join(formatted[4 * i:4 * i + 4])
+        if i < (lines - 1):
+            temp += ',\n'
+        else:
+            temp += '\n'
+        out += temp
+    out += '};'
+    return out
+
+
+with open('./ziggurat_constants.h', 'w') as f:
+    f.write(write(ki_final['double'], 'ki_double', 'uint64'))
+    f.write('\n\n')
+    f.write(write(wi_final['double'], 'wi_double', 'double'))
+    f.write('\n\n')
+    f.write(write(fi_final['double'], 'fi_double', 'double'))
+
+    f.write('\n\n\n\n')
+    f.write(write(ki_final['float'], 'ki_float', 'uint32'))
+    f.write('\n\n')
+    f.write(write(wi_final['float'], 'wi_float', 'float'))
+    f.write('\n\n')
+    f.write(write(fi_final['float'], 'fi_float', 'float'))