Merge pull request #130 from cactusdynamics/xorshift

Real-time-safe random number generator.

CMakeLists.txt

@@ -190,6 +190,7 @@ if(${CMAKE_PROJECT_NAME} STREQUAL ${PROJECT_NAME})
     add_subdirectory(examples/signal_handling_example)
     add_subdirectory(examples/simple_deadline_example)
     add_subdirectory(examples/simple_example)
+    add_subdirectory(examples/random_example)
 
     if (ENABLE_TRACING)
       add_subdirectory(examples/tracing_protos_example)

examples/random_example/CMakeLists.txt

@@ -0,0 +1,11 @@
+add_executable(random_example
+  main.cc
+)
+
+target_link_libraries(random_example
+  PRIVATE
+  cactus_rt
+)
+
+setup_cactus_rt_target_options(random_example)
+

examples/random_example/main.cc

@@ -0,0 +1,48 @@
+#include <cactus_rt/experimental/random.h>
+
+#include <array>
+#include <iomanip>
+#include <iostream>
+#include <random>
+
+template <size_t N>
+class Histogram {
+  std::array<size_t, N> hist_;
+
+ public:
+  Histogram() : hist_({0}) {}
+
+  void Record(float value) {
+    const auto i = static_cast<size_t>(value * N);
+    hist_.at(i)++;
+  }
+
+  void Display() {
+    constexpr float width = 1.0F / static_cast<float>(N);
+    float           current_bucket = 0.0F;
+    for (size_t i = 0; i < N; i++) {
+      std::cout << std::setprecision(4) << current_bucket << ": " << hist_[i] << "\n";
+      current_bucket += width;
+    }
+  }
+};
+
+int main() {
+  const uint64_t seed = std::random_device{}();
+  std::cout << "Seed: " << seed << "\n";
+
+  Histogram<20> hist;
+
+  cactus_rt::experimental::Xorshift64Rand rng(seed);
+
+  for (int i = 0; i < 1'000'000; i++) {
+    const float num = cactus_rt::experimental::RandomRealNumber(rng);
+    if (num >= 1.0F || num < 0.0F) {
+      std::cerr << "ERROR: seed = " << seed << " i = " << i << " num = " << num << " is out of range \n";
+      return 1;
+    }
+    hist.Record(num);
+  }
+  hist.Display();
+  return 0;
+}

include/cactus_rt/experimental/random.h

@@ -0,0 +1,77 @@
+#ifndef CACTUS_RT_EXPERIMENTAL_RAND_H_
+#define CACTUS_RT_EXPERIMENTAL_RAND_H_
+
+#include <cstdint>
+#include <limits>
+
+/**
+ * C++'s random number engines and random number distributions are amortized
+ * O(1) which as we know is greater than O(1) (theoretically O(inf)? But
+ * statistically very unlikely?). See Real-time Programming with the C++
+ * Standard Library - Timur Doumler - CppCon 2021[1].
+ *
+ * Thus we implement the Xorshift algorithm[2][3] to generate random numbers.
+ * This is not a cryptographically safe random number generator. Notably the
+ * uniform distribution implemented here do not guarantee perfect uniformity as
+ * it never discard numbers to generate another number (to ensure we can
+ * generate in O(1) and not amortized O(1)).
+ *
+ * [1]: https://youtu.be/Tof5pRedskI?t=2514
+ * [2]: https://en.wikipedia.org/wiki/Xorshift
+ * [3]: https://doi.org/10.18637/jss.v008.i14
+ */
+
+namespace cactus_rt::experimental {
+class Xorshift64Rand {
+  uint64_t x_;
+
+ public:
+  using result_type = uint64_t;
+
+  // Xorshift cannot have an initial state of 0. So we set it to 4 as it was chosen by a random die.
+  // (https://xkcd.com/221/)
+  explicit Xorshift64Rand(result_type initial_state) : x_(initial_state == 0 ? 4 : initial_state) {
+  }
+
+  result_type operator()() {
+    x_ ^= (x_ << 13);
+    x_ ^= (x_ >> 7);
+    x_ ^= (x_ << 17);
+    return x_;
+  };
+
+  static constexpr result_type min() {
+    return 1;
+  }
+
+  static constexpr result_type max() {
+    return std::numeric_limits<uint64_t>::max();
+  }
+};
+
+/**
+ * @brief Return a random number between [0, 1). Similar to
+ * std::uniform_real_distribution but not an object as it has no state. This is
+ * not a perfect uniform distribution and has some minor amount of bias, which
+ * is OK for real-time usage. It will also repeat as it doesn't allow you
+ * reseed.
+ *
+ * @tparam T The data type of the return result, default to float.
+ * @tparam Generator The random engine, default to Xorshift64Rand which is real-time safe.
+ * @param rng The RNG generator instance.
+ * @return T A random number between [0, 1)
+ */
+template <typename T = float, typename Generator = Xorshift64Rand>
+T RandomRealNumber(Generator& rng) {
+  T v = static_cast<T>(rng() - Generator::min()) / static_cast<T>(Generator::max() - Generator::min());
+  if (v == static_cast<T>(1.0)) {
+    // Random numbers are supposed to be [0, 1). This is a hack to make sure we never land on 1.
+    return static_cast<T>(0.0);
+  }
+
+  return v;
+}
+
+}  // namespace cactus_rt::experimental
+
+#endif

tests/CMakeLists.txt

@@ -4,6 +4,7 @@ enable_testing()
 
 add_executable(cactus_rt_tests
   utils_test.cc
+  experimental/random_test.cc
   experimental/lockless/atomic_bitset_test.cc
   experimental/lockless/spsc/realtime_readable_value_test.cc
   experimental/lockless/spsc/realtime_writable_value_test.cc

tests/experimental/random_test.cc

@@ -0,0 +1,66 @@
+#include "cactus_rt/experimental/random.h"
+
+#include <gtest/gtest.h>
+
+#include <random>
+
+using cactus_rt::experimental::RandomRealNumber;
+using cactus_rt::experimental::Xorshift64Rand;
+
+TEST(RandomRealNumber, Generate) {
+  const uint64_t seed = std::random_device{}();
+  Xorshift64Rand rng(seed);
+
+  for (int i = 0; i < 1'000'000; i++) {
+    const float current = RandomRealNumber(rng);
+    if (current < 0.0F || current >= 1.0F) {
+      ADD_FAILURE() << "number generated out of range: " << current << " (seed = " << seed << ", i = " << i << ")";
+    }
+  }
+
+  for (int i = 0; i < 1'000'000; i++) {
+    const auto current = RandomRealNumber<double>(rng);
+    if (current < 0.0 || current >= 1.0) {
+      ADD_FAILURE() << "number generated out of range: " << current << " (seed = " << seed << ", i = " << i << ")";
+    }
+  }
+}
+
+TEST(RandomRealNumber, GenerateZeroSeed) {
+  Xorshift64Rand rng(0);
+
+  for (int i = 0; i < 1'000'000; i++) {
+    const float current = RandomRealNumber(rng);
+    if (current < 0.0F || current >= 1.0F) {
+      ADD_FAILURE() << "number generated out of range: " << current << " (seed = " << 0 << ", i = " << i << ")";
+    }
+  }
+
+  for (int i = 0; i < 1'000'000; i++) {
+    const auto current = RandomRealNumber<double>(rng);
+    if (current < 0.0 || current >= 1.0) {
+      ADD_FAILURE() << "number generated out of range: " << current << " (seed = " << 0 << ", i = " << i << ")";
+    }
+  }
+}
+
+TEST(RandomRealNumber, DoesNotGenerate1) {
+  struct MaxGenerator {
+    using result_type = uint64_t;
+
+    static constexpr result_type max() {
+      return std::numeric_limits<uint64_t>::max();
+    }
+
+    static constexpr result_type min() {
+      return 1;
+    }
+
+    result_type operator()() {
+      return std::numeric_limits<uint64_t>::max();
+    }
+  };
+
+  MaxGenerator rng;
+  EXPECT_EQ(RandomRealNumber(rng), 0.0F);
+}