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);
+}