extmod/modrandom: Add a 64-bits native PRNG for 64-bits systems.

Inizia con extmod/modrandom.c e con le macro di configurazione e l’implementazione del PRNG proposte nell’issue. Traccia come getrandbits, l’inizializzazione del seed, randrange e la generazione in virgola mobile usano il generatore corrente su target a 32 bit e a 64 bit. Il lavoro è completato quando un generatore nativo a 64 bit attivabile opzionalmente funziona sulle architetture previste, mentre Yasmarang rimane quello predefinito per compatibilità.

Description

Right now MicroPython is using a PRNG that is meant to emit numbers in the 0..2³²-1 range (Yasmarang). Whilst this is fine on most MicroPython targets, for 64-bits architecture this may be a problem. See also #6718 for an existing description of other potential issues.

For example, if I got my maths right, a naïve implementation of a 64-bits PRNG calling Yasmarang twice to fill both halves of a 64-bits word would have one chance in 2³² to emit values that have the top 32 bits set to 0 (and similarly increasing probabilities for smaller upper contiguous sequences of zeroes). This can be fine for certain situations, but there are better algorithms out there that strive to be uniform across the whole 0..2⁶⁴-1 range.

https://prng.di.unimi.it/ contains a selection of fast, lean, and portable generators that pass most statistical randomness tests, and it is my understanding that those generators have similar capabilities as Yasmarang - albeit with larger periods for 64-bits generators. It could be useful to add one of those 64-bits native generators for 64-bits architectures on a opt-in basis (ie. 64-bits Unix targets, Win64, and QEMU/RV64), and keep Yasmarang around as a default to maintain compatibility.

Code Size

As a test I've added a xoroshiro128** implementation to extmod/modrandom.c and on my Linux/x64 machine the size increase is 128 bytes [1]. Most of that extra space is taken by the state seeder method suggested by the PRNG's paper, as it uses a few 64-bits constants and the compiler unrolls the small count iterations loop as well (#pragma GCC unroll 1/_Pragma("GCC unroll 1") doesn't seem to work). Still, ~128 bytes is relatively little overhead for an already larger binary when compared to the rest of the supported targets.

[1]

diff --git a/extmod/modrandom.c b/extmod/modrandom.c
index 79a1b18ba..cacce8b8c 100644
--- a/extmod/modrandom.c
+++ b/extmod/modrandom.c
@@ -38,6 +38,17 @@
 #define SEED_ON_IMPORT (0)
 #endif
 
+#if MICROPY_PY_RANDOM_YASMARANG && MICROPY_PY_RANDOM_XOROSHIRO128SS
+#error "Must pick only one PRNG algorithm."
+#endif
+
+// If Xoshiro128** is not explicitly defined use Yasmarang as a default.
+#if !MICROPY_PY_RANDOM_YASMARANG && !MICROPY_PY_RANDOM_XOROSHIRO128SS
+#define MICROPY_PY_RANDOM_YASMARANG (1)
+#endif
+
+#if MICROPY_PY_RANDOM_YASMARANG
+
 // Yasmarang random number generator
 // by Ilya Levin
 // http://www.literatecode.com/yasmarang
@@ -67,36 +78,103 @@ static uint32_t yasmarang(void) {
 
 // End of Yasmarang
 
+#define RNDINITSEED(seed) \
+    do { \
+        yasmarang_pad = (uint32_t)seed; \
+        yasmarang_n = 69; \
+        yasmarang_d = 233; \
+        yasmarang_dat = 0; \
+    } while (0)
+#define RND() yasmarang()
+#define RNDBITS 32
 #if MICROPY_PY_RANDOM_EXTRA_FUNCS
+#define RNDTYPE uint32_t
+#if MICROPY_PY_BUILTINS_FLOAT
+#define RND32() yasmarang()
+#define RND64() (((uint64_t)yasmarang() << 32) | (uint64_t)yasmarang())
+#endif
+#endif
 
-// returns an unsigned integer below the given argument
-// n must not be zero
-static uint32_t yasmarang_randbelow(uint32_t n) {
-    uint32_t mask = 1;
-    while ((n & mask) < n) {
-        mask = (mask << 1) | 1;
-    }
-    uint32_t r;
-    do {
-        r = yasmarang() & mask;
-    } while (r >= n);
-    return r;
+#elif MICROPY_PY_RANDOM_XOROSHIRO128SS
+
+#if !MICROPY_ENABLE_DYNRUNTIME
+#if SEED_ON_IMPORT
+// If the state is seeded on import then keep these variables in the BSS.
+static uint64_t xoroshiro128_state[2];
+#else
+// Without seed-on-import these variables must be initialised via the data section.
+static uint64_t xoroshiro128_state[2] = {
+    // Generated with splitmix64 using 0xbba490ef784fbcbf as the start value.
+    0xE8244547623F34E2, 0x2EAE10A6E969FC43
+};
+#endif
+#endif
+
+#if __has_builtin(__builtin_stdc_rotate_left)
+#define XOROSHIRO_ROL(x, k) __builtin_stdc_rotate_left(x, k)
+#else
+static uint64_t xoroshiro_rol(uint64_t x, uint64_t k) {
+    return (x << k) | (x >> (64 - k));
 }
+#define XOROSHIRO_ROL(x, k) xoroshiro_rol(x, k)
+#endif
+
+// Xoroshiro128** by David Blackman and Sebastiano Vigna
+// https://prng.di.unimi.it/xoroshiro128starstar.c
+// Original implementation released as public domain
+static uint64_t xoroshiro128starstar() {
+    const uint64_t s0 = xoroshiro128_state[0];
+    uint64_t s1 = xoroshiro128_state[1];
+    const uint64_t result = XOROSHIRO_ROL(s0 * 5, 7) * 9;
+
+    s1 ^= s0;
+    xoroshiro128_state[0] = XOROSHIRO_ROL(s0, 24) ^ s1 ^ (s1 << 16);
+    xoroshiro128_state[1] = XOROSHIRO_ROL(s1, 37);
+
+    return result;
+}
+
+// Splitmix64 by Sebastiano Vigna
+// https://prng.di.unimi.it/splitmix64.c
+// Original implementation released as public domain
+#define RNDINITSEED(seed) \
+    do { \
+        uint64_t x = (uint64_t)seed; \
+        for (size_t i = 0; i < 2; ++i) { \
+            uint64_t z = (x += 0x9e3779b97f4a7c15); \
+            z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9; \
+            z = (z ^ (z >> 27)) * 0x94d049bb133111eb; \
+            xoroshiro128_state[i] = z ^ (z >> 31); \
+        } \
+    } while (0)
+#define RND() xoroshiro128starstar()
+#define RNDBITS 64
+#if MICROPY_PY_RANDOM_EXTRA_FUNCS
+#define RNDTYPE uint64_t
+#if MICROPY_PY_BUILTINS_FLOAT
+#define RND32() ((uint32_t)(xoroshiro128starstar() >> 32))
+#define RND64() xoroshiro128starstar()
+#endif
+#endif
+
+#else
+
+#error "Please define either MICROPY_PY_RANDOM_YASMARANG or MICROPY_PY_RANDOM_XOROSHIRO128SS."
 
 #endif
 
 static mp_obj_t mod_random_getrandbits(mp_obj_t num_in) {
     mp_int_t n = mp_obj_get_int(num_in);
-    if (n > 32 || n < 0) {
-        mp_raise_ValueError(MP_ERROR_TEXT("bits must be 32 or less"));
+    if (n > RNDBITS || n < 0) {
+        mp_raise_ValueError(MP_ERROR_TEXT("bits must be " MP_STRINGIFY(RNDBITS) " or less"));
     }
     if (n == 0) {
         return MP_OBJ_NEW_SMALL_INT(0);
     }
-    uint32_t mask = ~0;
+    RNDTYPE mask = (RNDTYPE)-1;
     // Beware of C undefined behavior when shifting by >= than bit size
-    mask >>= (32 - n);
-    return mp_obj_new_int_from_uint(yasmarang() & mask);
+    mask >>= (RNDBITS - n);
+    return mp_obj_new_int_from_uint(RND() & mask);
 }
 static MP_DEFINE_CONST_FUN_OBJ_1(mod_random_getrandbits_obj, mod_random_getrandbits);
 
@@ -111,22 +189,33 @@ static mp_obj_t mod_random_seed(size_t n_args, const mp_obj_t *args) {
     } else {
         seed = mp_obj_get_int_truncated(args[0]);
     }
-    yasmarang_pad = (uint32_t)seed;
-    yasmarang_n = 69;
-    yasmarang_d = 233;
-    yasmarang_dat = 0;
+    RNDINITSEED(seed);
     return mp_const_none;
 }
 static MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(mod_random_seed_obj, 0, 1, mod_random_seed);
 
 #if MICROPY_PY_RANDOM_EXTRA_FUNCS
 
+// returns an unsigned integer below the given argument
+// n must not be zero
+static RNDTYPE randbelow(RNDTYPE n) {
+    RNDTYPE mask = 1;
+    while ((n & mask) < n) {
+        mask = (mask << 1) | 1;
+    }
+    RNDTYPE r;
+    do {
+        r = RND() & mask;
+    } while (r >= n);
+    return r;
+}
+
 static mp_obj_t mod_random_randrange(size_t n_args, const mp_obj_t *args) {
     mp_int_t start = mp_obj_get_int(args[0]);
     if (n_args == 1) {
         // range(stop)
         if (start > 0) {
-            return mp_obj_new_int(yasmarang_randbelow((uint32_t)start));
+            return mp_obj_new_int(randbelow(start));
         } else {
             goto error;
         }
@@ -135,7 +224,7 @@ static mp_obj_t mod_random_randrange(size_t n_args, const mp_obj_t *args) {
         if (n_args == 2) {
             // range(start, stop)
             if (start < stop) {
-                return mp_obj_new_int(start + yasmarang_randbelow((uint32_t)(stop - start)));
+                return mp_obj_new_int(start + randbelow(stop - start));
             } else {
                 goto error;
             }
@@ -151,7 +240,7 @@ static mp_obj_t mod_random_randrange(size_t n_args, const mp_obj_t *args) {
                 goto error;
             }
             if (n > 0) {
-                return mp_obj_new_int(start + step * yasmarang_randbelow((uint32_t)n));
+                return mp_obj_new_int(start + step * randbelow(n));
             } else {
                 goto error;
             }
@@ -167,7 +256,7 @@ static mp_obj_t mod_random_randint(mp_obj_t a_in, mp_obj_t b_in) {
     mp_int_t a = mp_obj_get_int(a_in);
     mp_int_t b = mp_obj_get_int(b_in);
     if (a <= b) {
-        return mp_obj_new_int(a + yasmarang_randbelow((uint32_t)(b - a + 1)));
+        return mp_obj_new_int(a + randbelow(b - a + 1));
     } else {
         mp_raise_ValueError(NULL);
     }
@@ -177,7 +266,7 @@ static MP_DEFINE_CONST_FUN_OBJ_2(mod_random_randint_obj, mod_random_randint);
 static mp_obj_t mod_random_choice(mp_obj_t seq) {
     mp_int_t len = mp_obj_get_int(mp_obj_len(seq));
     if (len > 0) {
-        return mp_obj_subscr(seq, mp_obj_new_int(yasmarang_randbelow((uint32_t)len)), MP_OBJ_SENTINEL);
+        return mp_obj_subscr(seq, mp_obj_new_int(randbelow(len)), MP_OBJ_SENTINEL);
     } else {
         mp_raise_type(&mp_type_IndexError);
     }
@@ -187,27 +276,27 @@ static MP_DEFINE_CONST_FUN_OBJ_1(mod_random_choice_obj, mod_random_choice);
 #if MICROPY_PY_BUILTINS_FLOAT
 
 // returns a number in the range [0..1) using Yasmarang to fill in the fraction bits
-static mp_float_t yasmarang_float(void) {
+static mp_float_t random_float(void) {
     mp_float_union_t u;
     u.p.sgn = 0;
     u.p.exp = (1 << (MP_FLOAT_EXP_BITS - 1)) - 1;
     if (MP_FLOAT_FRAC_BITS <= 32) {
-        u.p.frc = yasmarang();
+        u.p.frc = RND32();
     } else {
-        u.p.frc = ((uint64_t)yasmarang() << 32) | (uint64_t)yasmarang();
+        u.p.frc = RND64();
     }
     return u.f - 1;
 }
 
 static mp_obj_t mod_random_random(void) {
-    return mp_obj_new_float(yasmarang_float());
+    return mp_obj_new_float(random_float());
 }
 static MP_DEFINE_CONST_FUN_OBJ_0(mod_random_random_obj, mod_random_random);
 
 static mp_obj_t mod_random_uniform(mp_obj_t a_in, mp_obj_t b_in) {
     mp_float_t a = mp_obj_get_float(a_in);
     mp_float_t b = mp_obj_get_float(b_in);
-    return mp_obj_new_float(a + (b - a) * yasmarang_float());
+    return mp_obj_new_float(a + (b - a) * random_float());
 }
 static MP_DEFINE_CONST_FUN_OBJ_2(mod_random_uniform_obj, mod_random_uniform);

Implementation

I intend to implement this feature and would submit a Pull Request if desirable

Code of Conduct

Yes, I agree