ENH: Improve double generation for 64-bit generators

Improve random double generation when using 64 bit generators
by avoiding a division.

Author: bashtage

randomstate/interface/mlfg-1279-861/mlfg-1279-861-shim.h

@@ -40,11 +40,8 @@ inline uint64_t random_raw_values(aug_state* state)
 inline double random_double(aug_state* state)
 {
     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;
+    return (rn >> 11) * (1.0 / 9007199254740992.0);
 }
 
 extern void set_seed(aug_state* state, uint64_t seed);

randomstate/interface/pcg-64/pcg-64-shim.h

@@ -63,9 +63,6 @@ inline void entropy_init(aug_state* state)
 inline double random_double(aug_state* state)
 {
     uint64_t rn;
-    int32_t a, b;
     rn = random_uint64(state);
-    a = rn >> 37;
-    b = (rn & 0xFFFFFFFFLL) >> 6;
-    return (a * 67108864.0 + b) / 9007199254740992.0;
+    return (rn >> 11) * (1.0 / 9007199254740992.0);
 }

randomstate/interface/xoroshiro128plus/xoroshiro128plus-shim.h

@@ -47,11 +47,8 @@ inline uint64_t random_raw_values(aug_state* state)
 inline double random_double(aug_state* state)
 {
     uint64_t rn;
-    int32_t a, b;
     rn = random_uint64(state);
-    a = rn >> 37;
-    b = (rn & 0xFFFFFFFFLL) >> 6;
-    return (a * 67108864.0 + b) / 9007199254740992.0;
+    return (rn >> 11) * (1.0 / 9007199254740992.0);
 }
 
 extern void set_seed(aug_state* state, uint64_t seed);

randomstate/interface/xorshift1024/xorshift1024-shim.h

@@ -47,11 +47,8 @@ inline uint64_t random_raw_values(aug_state* state)
 inline double random_double(aug_state* state)
 {
     uint64_t rn;
-    int32_t a, b;
     rn = random_uint64(state);
-    a = rn >> 37;
-    b = (rn & 0xFFFFFFFFLL) >> 6;
-    return (a * 67108864.0 + b) / 9007199254740992.0;
+    return (rn >> 11) * (1.0 / 9007199254740992.0);
 }
 
 extern void set_seed(aug_state* state, uint64_t seed);

randomstate/interface/xorshift128/xorshift128-shim.h

@@ -47,11 +47,8 @@ inline uint64_t random_raw_values(aug_state* state)
 inline double random_double(aug_state* state)
 {
     uint64_t rn;
-    int32_t a, b;
     rn = random_uint64(state);
-    a = rn >> 37;
-    b = (rn & 0xFFFFFFFFLL) >> 6;
-    return (a * 67108864.0 + b) / 9007199254740992.0;
+    return (rn >> 11) * (1.0 / 9007199254740992.0);
 }
 
 extern void set_seed(aug_state* state, uint64_t seed);

randomstate/tests/test_direct.py

@@ -26,15 +26,9 @@ def uniform_from_uint(x, bits):
         return uniform_from_uint64(x)
     elif bits == 32:
         return uniform_from_uint32(x)
-    elif bits == 63:
-        return uniform_from_uint64(x)
-
 
 def uniform_from_uint64(x):
-    a = x >> 37
-    b = (x & 0xFFFFFFFF) >> 6
-    return (a * 67108864.0 + b) / 9007199254740992.0
-
+    return (x >> np.uint64(11)) * (1.0 / 9007199254740992.0)
 
 def uniform_from_uint32(x):
     out = np.empty(len(x) // 2)
@@ -44,7 +38,6 @@ def uniform_from_uint32(x):
         out[i // 2] = (a * 67108864.0 + b) / 9007199254740992.0
     return out
 
-
 def uint64_from_uint63(x):
     out = np.empty(len(x) // 2, dtype=np.uint64)
     for i in range(0, len(x), 2):