Merge pull request #13 from bashtage/mlfg-64

ENH: Update MLFG to 64 bits

Author: bashtage

README.md

@@ -35,10 +35,9 @@ identical sequence of random numbers for a given seed.
 ## Plans
 It is essentially complete.  There are a few rough edges that need to be smoothed.
 
-  * Document core RNG classes
   * Pickling support
-  * Verify entropy based initialization is missing for some RNGs
-  * Multiple stream support for MLFG and MRG32K3A
+  * Verify entropy based initialization for some RNGs
+  * Stream support for MLFG and MRG32K3A
   * Creation of additional streams from a RandomState where supported (i.e. 
   a `next_stream()` method)
 
@@ -134,36 +133,36 @@ Performance is promising.  Some early numbers:
 ```
 Time to produce 1,000,000 Uniforms
 ************************************************************
-numpy-random-random_sample                 11.44 ms
-randomstate-mlfg_1279_861-random_sample     7.18 ms
-randomstate-mrg32k3a-random_sample         35.70 ms
-randomstate-mt19937-random_sample           9.26 ms
-randomstate-pcg32-random_sample             6.75 ms
-randomstate-pcg64-random_sample             5.79 ms
-randomstate-xorshift1024-random_sample      5.73 ms
-randomstate-xorshift128-random_sample       5.38 ms
+numpy-random-random_sample                 13.68 ms
+randomstate-mlfg_1279_861-random_sample     6.64 ms
+randomstate-mrg32k3a-random_sample         37.87 ms
+randomstate-mt19937-random_sample          13.33 ms
+randomstate-pcg32-random_sample            10.20 ms
+randomstate-pcg64-random_sample             7.83 ms
+randomstate-xorshift1024-random_sample      6.20 ms
+randomstate-xorshift128-random_sample       5.49 ms
 
 Uniforms per second
 ************************************************************
-numpy-random-random_sample                  87.39 million
-randomstate-mlfg_1279_861-random_sample    139.27 million
-randomstate-mrg32k3a-random_sample          28.01 million
-randomstate-mt19937-random_sample          107.94 million
-randomstate-pcg32-random_sample            148.16 million
-randomstate-pcg64-random_sample            172.63 million
-randomstate-xorshift1024-random_sample     174.63 million
-randomstate-xorshift128-random_sample      185.83 million
-dtype: object
+numpy-random-random_sample                  73.11 million
+randomstate-mlfg_1279_861-random_sample    150.71 million
+randomstate-mrg32k3a-random_sample          26.41 million
+randomstate-mt19937-random_sample           75.03 million
+randomstate-pcg32-random_sample             98.00 million
+randomstate-pcg64-random_sample            127.77 million
+randomstate-xorshift1024-random_sample     161.39 million
+randomstate-xorshift128-random_sample      182.29 million
 
 Speed-up relative to NumPy
 ************************************************************
-randomstate-mlfg_1279_861-random_sample     59.4%
-randomstate-mrg32k3a-random_sample         -67.9%
-randomstate-mt19937-random_sample           23.5%
-randomstate-pcg32-random_sample             69.5%
-randomstate-pcg64-random_sample             97.5%
-randomstate-xorshift1024-random_sample      99.8%
-randomstate-xorshift128-random_sample      112.6%
+randomstate-mlfg_1279_861-random_sample    106.1%
+randomstate-mrg32k3a-random_sample         -63.9%
+randomstate-mt19937-random_sample            2.6%
+randomstate-pcg32-random_sample             34.0%
+randomstate-pcg64-random_sample             74.8%
+randomstate-xorshift1024-random_sample     120.7%
+randomstate-xorshift128-random_sample      149.3%
+
 --------------------------------------------------------------------------------
 ```
 

randomstate/__init__.py

@@ -1,29 +1,9 @@
-#import randomstate.mlfg_1279_861 as mlfg_1279_861
-#import randomstate.mrg32k3a as mrg32k3a
-import randomstate.mt19937 as __mt19937
-#import randomstate.pcg32 as pcg32
-#import randomstate.xorshift1024 as xorshift1024
-#import randomstate.xorshift128 as xorshift128
+from __future__ import division, absolute_import, print_function
 
-#try:
-#    import randomstate.pcg64
-#except ImportError:
-#    pass
-
-__all__ = ['standard_normal', 'standard_gamma', 'get_state', 'set_state', 'seed', 'bytes',
-           'standard_exponential', 'standard_cauchy', 'standard_t', 'rand', 'randn']
-
-__rs = __mt19937.RandomState()
-RandomState = __mt19937.RandomState
-standard_normal = __rs.standard_normal
-standard_gamma = __rs.standard_gamma
-standard_exponential = __rs.standard_exponential
-standard_cauchy = __rs.standard_cauchy
-standard_t = __rs.standard_t
-random_sample = __rs.random_sample
-get_state = __rs.get_state
-set_state = __rs.set_state
-seed = __rs.seed
-bytes = __rs.bytes
-randn = __rs.randn
-rand = __rs.rand
+from randomstate.mt19937 import *
+import randomstate.mlfg_1279_861 as mlfg_1279_861
+import randomstate.mrg32k3a as mrg32k3a
+import randomstate.pcg32 as pcg32
+import randomstate.pcg64 as pcg64
+import randomstate.xorshift1024 as xorshift1024
+import randomstate.xorshift128 as xorshift128

randomstate/setup-all.py

@@ -1,3 +1,4 @@
+
 #ifdef _WIN32
 #include "../../src/common/stdint.h"
 #define inline __inline
@@ -9,10 +10,6 @@
 #include "../../src/common/binomial.h"
 #include "../../src/mlfg-1279-861/mlfg-1279-861.h"
 
-#define BITMASK_UPPER 0xffff0000
-#define BITMAST_UPPER22 0xfffffc00
-#define BITMAST_UPPER21 0xfffff800
-
 typedef struct s_aug_state {
     mlfg_state *rng;
     binomial_t *binomial;
@@ -25,24 +22,13 @@ typedef struct s_aug_state {
 
 inline uint32_t random_uint32(aug_state* state)
 {
-    // Two are needed since there is only 31 bits in each
-#if defined(_WIN32) && !defined(_WIN64)
-   uint32_t a, b;
-    a = mlfg_next(state->rng) & BITMASK_UPPER;
-    b = mlfg_next(state->rng) >> 16;
-    return a | b;
-#else
-    return (mlfg_next(state->rng) & BITMASK_UPPER) | (mlfg_next(state->rng) >> 16);
-#endif
+    return (uint32_t)(mlfg_next(state->rng) >> 32);
 }
 
 inline uint64_t random_uint64(aug_state* state)
 {
-
-    uint64_t out = ((uint64_t)(mlfg_next(state->rng)) & BITMAST_UPPER22) << 32;
-    out |= ((uint64_t)(mlfg_next(state->rng)) & BITMAST_UPPER21) << 10;
-    out |= ((mlfg_next(state->rng)) & BITMAST_UPPER21) >> 11;
-    // Three are needed since there is only 31 bits in each
+    uint64_t out  = mlfg_next(state->rng) & 0xffffffff00000000ULL;
+    out |= mlfg_next(state->rng) >> 32;
     return out;
 }
 
@@ -60,7 +46,10 @@ inline void entropy_init(aug_state* state)
 
 inline double random_double(aug_state* state)
 {
-    // These differ by 1 bit since there are only 31 bits in each
-    int32_t a = mlfg_next(state->rng) >> 5, b = mlfg_next(state->rng) >> 6;
+    uint64_t rn;
+    int32_t a, b;
+    rn = mlfg_next(state->rng);
+    a = rn >> 37;
+    b = (rn & 0xFFFFFFFFULL) >> 6;
     return (a * 67108864.0 + b) / 9007199254740992.0;
-}
+}