rng.hpp

// ==========================================================================
// Dedmonwakeen's Raid DPS/TPS Simulator.
// Send questions to [email protected]
// ==========================================================================

#pragma once

#include "config.hpp"

#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>

#include "util/timespan.hpp"

namespace rng {

template <typename Engine>
struct basic_rng_t
{
public:
  explicit basic_rng_t() = default;

  basic_rng_t(const basic_rng_t&) = delete;
  basic_rng_t& operator=(const basic_rng_t&) = delete;

  const char* name() const {
    return engine.name();
  }

  void seed( uint64_t s ) {
    engine.seed( s );
  }

  uint64_t reseed();

  void reset();

  double real();

  bool roll( double chance );

  double range( double min, double max );

  template <typename T, typename = std::enable_if_t<std::is_integral<T>::value>>
  T range( T min, T max ) {
    return static_cast<T>(range(static_cast<double>(min), static_cast<double>(max)));
  }

  template <typename T, typename = std::enable_if_t<std::is_integral<T>::value>>
  T range( T max ) {
    return range<T>( T{}, max );
  }

  double gauss( double mean, double stddev, bool truncate_low_end = false );

  double exponential( double nu );

  double exgauss( double gauss_mean, double gauss_stddev, double exp_nu );

  timespan_t range( timespan_t min, timespan_t max );

  timespan_t gauss( timespan_t mean, timespan_t stddev );

  timespan_t exgauss( timespan_t mean, timespan_t stddev, timespan_t nu );

private:
  Engine engine;
  // Allow re-use of unused ( but necessary ) random number of a previous call to gauss()
  double gauss_pair_value = 0.0;
  bool   gauss_pair_use = false;
};

template <typename Engine>
uint64_t basic_rng_t<Engine>::reseed()
{
  const uint64_t s = engine.next();
  seed( s );
  reset();
  return s;
}

template <typename Engine>
void basic_rng_t<Engine>::reset()
{
  gauss_pair_value = 0;
  gauss_pair_use = false;
}

// ==========================================================================
// Probability Distributions
// ==========================================================================

template <typename Engine>
inline double basic_rng_t<Engine>::real()
{
  uint64_t ui64 = engine.next();
  ui64 &= 0x000fffffffffffff;
  ui64 |= 0x3ff0000000000000;
  union { uint64_t ui64; double d; } u;
  u.ui64 = ui64;
  return u.d - 1.0;
}

template <typename Engine>
bool basic_rng_t<Engine>::roll( double chance )
{
  if ( chance <= 0 ) return false;
  if ( chance >= 1 ) return true;
  return real() < chance;
}

template <typename Engine>
double basic_rng_t<Engine>::range( double min, double max )
{
  assert( min <= max );
  return min + real() * ( max - min );
}

template <typename Engine>
double basic_rng_t<Engine>::gauss( double mean, double stddev, bool truncate_low_end )
{
  assert(stddev >= 0 && "Calling gauss with negative stddev");

  double z;

  if ( stddev != 0 )
  {
    if ( gauss_pair_use )
    {
      z = gauss_pair_value;
      gauss_pair_use = false;
    }
    else
    {
      double x1, x2, w, y1, y2;
      do
      {
        x1 = 2.0 * real() - 1.0;
        x2 = 2.0 * real() - 1.0;
        w = x1 * x1 + x2 * x2;
      }
      while ( w >= 1.0 || w == 0.0 );

      w = std::sqrt( ( -2.0 * std::log( w ) ) / w );
      y1 = x1 * w;
      y2 = x2 * w;

      z = y1;
      gauss_pair_value = y2;
      gauss_pair_use = true;
    }
  }
  else
    z = 0.0;

  double result = mean + z * stddev;

  // True gaussian distribution can of course yield any number at some probability.  So truncate on the low end.
  if ( truncate_low_end && result < 0 )
    result = 0;

  return result;
}

template <typename Engine>
double basic_rng_t<Engine>::exponential( double nu )
{
  double x;
  do { x = real(); } while ( x >= 1.0 ); // avoid ln(0)
  return - std::log( 1 - x ) * nu;
}

template <typename Engine>
double basic_rng_t<Engine>::exgauss( double gauss_mean, double gauss_stddev, double exp_nu )
{
  double v = gauss( gauss_mean, gauss_stddev ) + exponential( exp_nu );
  return v > 0.0 ? v : 0.0;
}

template <typename Engine>
timespan_t basic_rng_t<Engine>::range( timespan_t min, timespan_t max )
{
  return timespan_t::from_native( range( timespan_t::to_native( min ), timespan_t::to_native( max ) ) );
}

template <typename Engine>
timespan_t basic_rng_t<Engine>::gauss( timespan_t mean, timespan_t stddev )
{
  return timespan_t::from_native( gauss( static_cast<double>( timespan_t::to_native( mean ) ),
                                         static_cast<double>( timespan_t::to_native( stddev ) ) ) );
}

template <typename Engine>
timespan_t basic_rng_t<Engine>::exgauss( timespan_t mean, timespan_t stddev, timespan_t nu )
{
  return timespan_t::from_native( exgauss( static_cast<double>( timespan_t::to_native( mean   ) ),
                                           static_cast<double>( timespan_t::to_native( stddev ) ),
                                           static_cast<double>( timespan_t::to_native( nu ) ) ) );
}


struct xorshift128_t
{
  uint64_t next() noexcept;
  void seed( uint64_t start ) noexcept;
  const char* name() const noexcept;
private:
  std::array<uint64_t, 2> s;
};

struct xoshiro256plus_t
{
  uint64_t next() noexcept;
  void seed( uint64_t start ) noexcept;
  const char* name() const noexcept;
private:
  std::array<uint64_t, 4> s;
};

struct xorshift1024_t
{
  uint64_t next() noexcept;
  void seed( uint64_t start ) noexcept;
  const char* name() const noexcept;
private:
  std::array<uint64_t, 16> s;
  int p;
};

// "Default" rng
// Explicitly *NOT* a type alias to allow forward declaraions
struct rng_t : public basic_rng_t<xoshiro256plus_t> {};

double stdnormal_cdf( double );
double stdnormal_inv( double );

} // rng