random.h

Go to the documentation of this file.

 
#if !defined(ULTRA_RANDOM_H)
#define      ULTRA_RANDOM_H
 
#include "kernel/interval.h"
 
#include "utility/assert.h"
#include "utility/misc.h"
#include "utility/xoshiro256ss.h"
 
#include <random>
#include <type_traits>
 
namespace ultra::random
{
 
enum class distribution {uniform, normal};

using engine_t = vigna::xoshiro256ss;
 
[[nodiscard]] engine_t &engine();
[[nodiscard]] std::size_t ring(std::size_t, std::size_t, std::size_t);
void randomize();

template<std::floating_point T>
[[nodiscard]] T between(T min, T sup)
{
  Expects(min < sup);
 
  std::uniform_real_distribution<T> d(min, sup);
  return d(engine());
}

template<std::integral T>
[[nodiscard]] T between(T min, T sup)
{
  Expects(min < sup);
 
  std::uniform_int_distribution<T> d(min, sup - 1);
  return d(engine());
}
 
template<class T>
requires std::is_enum_v<T>
[[nodiscard]] T between(T min, T sup)
{
  Expects(min < sup);
 
  return static_cast<T>(between<std::underlying_type_t<T>>(min, sup));
}

template<class T>
requires std::is_arithmetic_v<T> || std::is_enum_v<T>
[[nodiscard]] T sup(T sup)
{
  return between(static_cast<T>(0), sup);
}

template<std::ranges::sized_range C>
[[nodiscard]] std::ranges::range_reference_t<const C> element(const C &c)
{
  Expects(c.size());
 
  return *std::next(
    c.begin(),
    static_cast<typename C::difference_type>(sup(c.size())));
}

template<ArithmeticScalar A>
[[nodiscard]] A element(const interval<A> &i)
{
  return between(i.min, i.sup);
}

template<std::ranges::sized_range C>
[[nodiscard]] typename C::value_type &element(C &c)
{
  Expects(c.size());
 
  return *std::next(
    c.begin(),
    static_cast<typename C::difference_type>(sup(c.size())));
}

[[nodiscard]] inline bool boolean(double p)
{
  Expects(in_0_1(p));
 
  std::bernoulli_distribution d(p);
  return d(engine());
}

[[nodiscard]] inline bool boolean()
{
  // Fast path for fair coin flips: cache one PRNG result and consume it one
  // bit at a time, amortising one engine call across all bits in
  // `engine_t::result_type`.
  class thread_local_coin_flip
  {
  public:
    [[nodiscard]] bool operator()()
    {
      if (!bits_left_)  [[unlikely]]
      {
        cache_ = random::engine()();
        bits_left_ = std::numeric_limits<word_t>::digits;
      }
 
      const bool result((cache_ >> (bits_left_ - 1)) & 1u);
      --bits_left_;
      return result;
    }
 
  private:
    using word_t = random::engine_t::result_type;
 
    word_t cache_ {random::engine()()};
    unsigned bits_left_ {std::numeric_limits<word_t>::digits};
  };
 
  thread_local thread_local_coin_flip coin_flip;
  return coin_flip();
}

template<class T>
requires std::is_arithmetic_v<T>
[[nodiscard]] T ephemeral(distribution d, T p1, T p2)
{
  Expects(p1 < p2);
 
  switch (d)
  {
  case distribution::normal:
    return std::normal_distribution<T>(std::midpoint(p1, p2),
                                       p2 - p1)(engine());
  default:
    return between(p1, p2);
  }
}
 
}  // namespace ultra::random
 
#endif  // include guard