runtime: fix gr::random API to be fixed-width 64 bit, use XOROSHIRO128+

Seeding being inconsistent between gr::random and things like fastnoise source is an outstanding issue (#1559). Fix that by fixing the API, and pivoting to our built-in XOROSHIRO128+, which is also substantially faster than MT19937. For that purpose, introduce, and python-bind, a wrapper class that can be used with STL distribution shapers. Add a unit test for the (P)RNG. Signed-off-by: Marcus Müller <mmueller@gnuradio.org>
author: Marcus Müller <mmueller@gnuradio.org> 2021-03-22 00:37:47 +0100
committer: mormj <34754695+mormj@users.noreply.github.com> 2021-04-06 07:37:17 -0400
commit: eb91fb04b3d0ca8124b9bea375e3bf72963fa172 (patch)
tree: dcf96c9e232b7019a7d81cd16901a95bda8b5471
parent: 42349ff2855c1bb0148c8191339189b0a4cee675 (diff)
6 files changed, 139 insertions, 18 deletions
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 270e5d14d5..a41496cc7c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -32,6 +32,10 @@ Older Logs can be found in `docs/RELEASE-NOTES-*`.
 - `fastnoise_source`: Use a simple bitmask if the random pool length is a power
   of 2 to determine indices, instead of `%`, which consumed considerable CPU
 
+#### gnuradio-runtime
+
+- `gr::random` uses xoroshiro128+ internally, takes `uint64_t` seed
+
 ### Added
 
 - New in-tree module gr-pdu
diff --git a/gnuradio-runtime/include/gnuradio/random.h b/gnuradio-runtime/include/gnuradio/random.h
index 504b36a933..ba216f3117 100644
--- a/gnuradio-runtime/include/gnuradio/random.h
+++ b/gnuradio-runtime/include/gnuradio/random.h
@@ -13,14 +13,53 @@
 
 #include <gnuradio/api.h>
 #include <gnuradio/gr_complex.h>
+#include <gnuradio/xoroshiro128p.h>
 
-#include <cstdlib>
-#include <ctime>
+#include <limits>
 #include <random>
 
 namespace gr {
 
 /*!
+ * \brief wrapper for XOROSHIRO128+ PRNG for use in std::distributions
+ * Fulfills C++ named requirements for UniformRandomBitGenerator
+ * \ingroup math_blk
+ */
+class GR_RUNTIME_API xoroshiro128p_prng
+{
+public:
+    using result_type = uint64_t; //! \brief value type is uint64
+
+private:
+    result_type state[2];
+
+public:
+    /*!
+     * \brief minimum value
+     */
+    static constexpr result_type min() { return std::numeric_limits<result_type>::min(); }
+    /*!
+     * \brief maximum value
+     */
+    static constexpr result_type max() { return std::numeric_limits<result_type>::max(); }
+
+    /*!
+     * \brief constructor. Expects a seed.
+     */
+    xoroshiro128p_prng(uint64_t init) { seed(init); }
+
+
+    /*!
+     * \brief yield a random value and advance state
+     */
+    result_type operator()() { return xoroshiro128p_next(state); }
+
+    /*!
+     * \brief set new seed
+     */
+    void seed(uint64_t seed) { xoroshiro128p_seed(state, seed); }
+};
+/*!
  * \brief pseudo random number generator
  * \ingroup math_blk
  */
@@ -31,32 +70,32 @@ protected:
     bool d_gauss_stored;
     float d_gauss_value;
 
-    std::mt19937 d_rng; // mersenne twister as random number generator
+    xoroshiro128p_prng d_rng; // mersenne twister as random number generator
     std::uniform_real_distribution<float>
         d_uniform; // choose uniform distribution, default is [0,1)
-    std::uniform_int_distribution<> d_integer_dis;
+    std::uniform_int_distribution<int64_t> d_integer_dis;
 
 public:
-    random(unsigned int seed = 0, int min_integer = 0, int max_integer = 2);
+    random(uint64_t seed = 0, int64_t min_integer = 0, int64_t max_integer = 2);
     ~random();
 
     /*!
      * \brief Change the seed for the initialized number generator. seed = 0 initializes
      * the random number generator with the system time.
      */
-    void reseed(unsigned int seed);
+    void reseed(uint64_t seed);
 
     /*!
      * set minimum and maximum for integer random number generator.
      * Limits are [minimum, maximum)
      * Default: [0, std::numeric_limits< IntType >::max)]
      */
-    void set_integer_limits(const int minimum, const int maximum);
+    void set_integer_limits(int64_t minimum, int64_t maximum);
 
     /*!
      * Uniform random integers in the range set by 'set_integer_limits' [min, max).
      */
-    int ran_int();
+    int64_t ran_int();
 
     /*!
      * \brief Uniform random numbers in the range [0.0, 1.0)
diff --git a/gnuradio-runtime/lib/math/random.cc b/gnuradio-runtime/lib/math/random.cc
index 8b67532207..6423a39f6c 100644
--- a/gnuradio-runtime/lib/math/random.cc
+++ b/gnuradio-runtime/lib/math/random.cc
@@ -35,13 +35,12 @@
 
 namespace gr {
 
-random::random(unsigned int seed, int min_integer, int max_integer)
-    : d_rng(), d_integer_dis(0, 1)
+random::random(uint64_t seed, int64_t min_integer, int64_t max_integer)
+    : d_rng(seed), d_integer_dis(0, 1)
 {
     d_gauss_stored = false; // set gasdev (gauss distributed numbers) on calculation state
 
     // Setup random number generators
-    reseed(seed); // set seed for random number generator
     set_integer_limits(min_integer, max_integer);
 }
 
@@ -51,7 +50,7 @@ random::~random() {}
  * Seed is initialized with time if the given seed is 0. Otherwise the seed is taken
  * directly. Sets the seed for the random number generator.
  */
-void random::reseed(unsigned int seed)
+void random::reseed(uint64_t seed)
 {
     d_seed = seed;
     if (d_seed == 0) {
@@ -63,17 +62,17 @@ void random::reseed(unsigned int seed)
     }
 }
 
-void random::set_integer_limits(const int minimum, const int maximum)
+void random::set_integer_limits(int64_t minimum, int64_t maximum)
 {
     // boost expects integer limits defined as [minimum, maximum] which is unintuitive.
     // use the expected half open interval behavior! [minimum, maximum)!
-    d_integer_dis = std::uniform_int_distribution<>(minimum, maximum - 1);
+    d_integer_dis = std::uniform_int_distribution<int64_t>(minimum, maximum - 1);
 }
 
 /*!
  * Uniform random integers in the range set by 'set_integer_limits' [min, max).
  */
-int random::ran_int() { return d_integer_dis(d_rng); }
+int64_t random::ran_int() { return d_integer_dis(d_rng); }
 
 /*
  * Returns uniformly distributed numbers in [0,1) taken from boost.random using a Mersenne
diff --git a/gnuradio-runtime/python/gnuradio/gr/bindings/docstrings/random_pydoc_template.h b/gnuradio-runtime/python/gnuradio/gr/bindings/docstrings/random_pydoc_template.h
index 71fdb84d39..d02cadf872 100644
--- a/gnuradio-runtime/python/gnuradio/gr/bindings/docstrings/random_pydoc_template.h
+++ b/gnuradio-runtime/python/gnuradio/gr/bindings/docstrings/random_pydoc_template.h
@@ -1,5 +1,5 @@
 /*
- * Copyright 2020 Free Software Foundation, Inc.
+ * Copyright 2021 Free Software Foundation, Inc.
  *
  * This file is part of GNU Radio
  *
@@ -15,6 +15,24 @@
  */
 
 
+static const char* __doc_gr_xoroshiro128p_prng = R"doc()doc";
+
+
+static const char* __doc_gr_xoroshiro128p_prng_xoroshiro128p_prng_0 = R"doc()doc";
+
+
+static const char* __doc_gr_xoroshiro128p_prng_xoroshiro128p_prng_1 = R"doc()doc";
+
+
+static const char* __doc_gr_xoroshiro128p_prng_min = R"doc()doc";
+
+
+static const char* __doc_gr_xoroshiro128p_prng_max = R"doc()doc";
+
+
+static const char* __doc_gr_xoroshiro128p_prng_seed = R"doc()doc";
+
+
 static const char* __doc_gr_random = R"doc()doc";
 
 
diff --git a/gnuradio-runtime/python/gnuradio/gr/bindings/random_python.cc b/gnuradio-runtime/python/gnuradio/gr/bindings/random_python.cc
index 8c823e5977..57fbbb2495 100644
--- a/gnuradio-runtime/python/gnuradio/gr/bindings/random_python.cc
+++ b/gnuradio-runtime/python/gnuradio/gr/bindings/random_python.cc
@@ -14,7 +14,7 @@
 /* BINDTOOL_GEN_AUTOMATIC(0)                                                       */
 /* BINDTOOL_USE_PYGCCXML(0)                                                        */
 /* BINDTOOL_HEADER_FILE(random.h)                                        */
-/* BINDTOOL_HEADER_FILE_HASH(95215055c6ad6f78b88d532d3994d29a)                     */
+/* BINDTOOL_HEADER_FILE_HASH(489c2917f344c75e2001d0c670de7784)                     */
 /***********************************************************************************/
 
 #include <pybind11/complex.h>
@@ -29,10 +29,35 @@ namespace py = pybind11;
 
 void bind_random(py::module& m)
 {
-
+    using xoroshiro128p_prng = ::gr::xoroshiro128p_prng;
     using random = ::gr::random;
 
 
+    py::class_<xoroshiro128p_prng, std::shared_ptr<xoroshiro128p_prng>>(
+        m, "xoroshiro128p_prng", D(xoroshiro128p_prng))
+
+        .def(py::init<uint64_t>(),
+             py::arg("init"),
+             D(xoroshiro128p_prng, xoroshiro128p_prng, 0))
+        .def(py::init<gr::xoroshiro128p_prng const&>(),
+             py::arg("arg0"),
+             D(xoroshiro128p_prng, xoroshiro128p_prng, 1))
+
+
+        .def_static("min", &xoroshiro128p_prng::min, D(xoroshiro128p_prng, min))
+
+
+        .def_static("max", &xoroshiro128p_prng::max, D(xoroshiro128p_prng, max))
+
+
+        .def("seed",
+             &xoroshiro128p_prng::seed,
+             py::arg("seed"),
+             D(xoroshiro128p_prng, seed))
+
+        .def("__call__", &xoroshiro128p_prng::operator());
+
+
     py::class_<random, std::shared_ptr<random>>(m, "random", D(random))
 
         .def(py::init<uint64_t, int64_t, int64_t>(),
diff --git a/gnuradio-runtime/python/gnuradio/gr/qa_random.py b/gnuradio-runtime/python/gnuradio/gr/qa_random.py
index fa8c44917b..d7f4bb65a5 100644
--- a/gnuradio-runtime/python/gnuradio/gr/qa_random.py
+++ b/gnuradio-runtime/python/gnuradio/gr/qa_random.py
@@ -63,6 +63,42 @@ class test_random(gr_unittest.TestCase):
         self.assertGreaterEqual(minimum, np.min(rnd_vals))
         self.assertLess(np.max(rnd_vals), maximum)
 
+    def test_005_xoroshiro128p_seed_stability(self):
+        """
+        Test that seeding is stable.
+        It's basically an API break if it isn't.
+
+        We simply check for the first value of a sequence
+        being the same as it was when the module was integrated.
+        """
+        rng = gr.xoroshiro128p_prng(42)
+        self.assertEqual(3520422898491873512, rng())
+
+    def test_006_xoroshiro128p_reproducibility(self):
+        """
+        Make sure two RNGs with the same seed yield the same
+        sequence
+        """
+        seed = 123456
+        N = 10000
+        rng1 = gr.xoroshiro128p_prng(123456)
+        rng2 = gr.xoroshiro128p_prng(123456)
+        self.assertSequenceEqual(
+            tuple(rng1() for _ in range(N)),
+            tuple(rng2() for _ in range(N)))
+
+    def test_007_xoroshiro128p_range(self):
+        """
+        Check bounds.
+        Check whether a long sequence of values are within that bounds.
+        """
+        N = 10**6
+
+        self.assertEqual(gr.xoroshiro128p_prng.min(), 0)
+        self.assertEqual(gr.xoroshiro128p_prng.max(), 2**64-1)
+        rng = gr.xoroshiro128p_prng(42)
+        arr = all((0 <= rng() <= 2**64 - 1 for _ in range(N)))
+        self.assertTrue(arr)
 
 if __name__ == '__main__':
     gr_unittest.run(test_random)
author	Marcus Müller <mmueller@gnuradio.org>	2021-03-22 00:37:47 +0100
committer	mormj <34754695+mormj@users.noreply.github.com>	2021-04-06 07:37:17 -0400
commit	eb91fb04b3d0ca8124b9bea375e3bf72963fa172 (patch)
tree	dcf96c9e232b7019a7d81cd16901a95bda8b5471
parent	42349ff2855c1bb0148c8191339189b0a4cee675 (diff)