random.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/util/assertion.hpp>
#include <kernel/util/type_traits.hpp>
 
// includes, system
#include <algorithm>
#include <limits>
#include <cstdint>
#include <chrono>
 
namespace FEAT
{
  namespace Intern
  {
    // Random number class; this class's gen() function takes a Random object reference and
    // uses its next() function to generate a random number of the specified type T_.
    template<typename T_, typename TypeClass_ = typename Type::Traits<T_>::TypeClass>
    class RandomNumber;
  } // namespace Intern
 
  class Random
  {
  public:
    typedef std::uint64_t SeedType;
 
    static constexpr SeedType xor_mult = SeedType(2685821657736338717ull);
 
  private:
    std::uint64_t _seed;
    std::uint64_t _x;
 
  public:
    explicit Random() :
      _seed(0),
      _x(0)
    {
      // use the high-resolution clock from the standard library go get the epoch time in nanoseconds as the seed
      auto now = std::chrono::high_resolution_clock::now();
      _seed = _x = std::chrono::duration_cast<std::chrono::duration<SeedType, std::nano>>(now.time_since_epoch()).count() + SeedType(17u);
    }
 
    explicit Random(SeedType seed) :
      _seed(seed),
      _x(seed)
    {
      XASSERTM(seed != SeedType(0), "Random seed must not be zero!");
    }
 
    SeedType get_seed() const
    {
      return _seed;
    }
 
    uint64_t next()
    {
      _x ^= _x >> 12;
      _x ^= _x << 25;
      _x ^= _x >> 27;
      return _x * xor_mult;
    }
 
    template<typename T_>
    Random& operator>>(T_& x)
    {
      x = FEAT::Intern::RandomNumber<T_>::gen(*this);
      return *this;
    }
 
    template<typename T_>
    T_ operator()(T_ a, T_ b)
    {
      return FEAT::Intern::RandomNumber<T_>::gen_ranged(*this, std::min(a, b), std::max(a, b));
    }
  }; // class Random
 
  namespace Intern
  {
    // Random integer class; this class is the base class for the integer specialization
    // of the RandomNumber class defined below.
    template<typename T_, size_t num_bytes_ = sizeof(T_)>
    class RandomInteger;
 
    // Helper class: return non-negative integer
    template<typename T_, bool signed_ = (Type::Traits<T_>::is_signed)>
    class NonNegInt
    {
    public:
      static T_ make(T_ x)
      {
        return x;
      }
    };
 
    // specialization for signed integer types
    template<typename T_>
    class NonNegInt<T_, true>
    {
    public:
      static T_ make(T_ x)
      {
        return (x < T_(0) ? T_(-(x + T_(1))) : x);
      }
    };
 
    // specialization for 8-bit integers
    template<typename T_>
    class RandomInteger<T_, 1>
    {
    public:
      static T_ gen(Random& rng)
      {
        return (T_)(rng.next() & 0xFFu);
      }
    };
 
    // specialization for 16-bit integers
    template<typename T_>
    class RandomInteger<T_, 2>
    {
    public:
      static T_ gen(Random& rng)
      {
        return (T_)(rng.next() & 0xFFFFu);
      }
    };
 
    // specialization for 32-bit integers
    template<typename T_>
    class RandomInteger<T_, 4>
    {
    public:
      static T_ gen(Random& rng)
      {
        return (T_)(rng.next() & 0xFFFFFFFFu);
      }
    };
 
    // specialization for 64-bit integers
    template<typename T_>
    class RandomInteger<T_, 8>
    {
    public:
      static T_ gen(Random& rng)
      {
        return (T_)(rng.next());
      }
    };
 
    // specialization for integers
    template<typename T_>
    class RandomNumber<T_, Type::IntegralClass> :
      public RandomInteger<T_>
    {
    public:
      static T_ gen_ranged(Random& rng, T_ a, T_ b)
      {
        // call gen() prior to the if-clause to advance the rng in any case
        // also, ensure that the value is non-negative
        T_ x = NonNegInt<T_>::make(RandomInteger<T_>::gen(rng));
        if(a < b)
          // Note: The following casting orgy is necessary on some compilers if the type 'T_' is
          // actually shorter than 'int', as the arithmetic operations return values of type 'int'
          // in this case.
          return T_(a + T_(x % T_(T_(b - a) + T_(1))));
        else
          return a;
      }
    };
 
    // specialization for floating point numbers
    template<typename T_>
    class RandomNumber<T_, Type::FloatingClass>
    {
    public:
      static T_ gen(Random& rng)
      {
        static constexpr std::uint64_t mask = std::uint64_t(0xFFFFFFFFull);
        const double d = 1.0 / double(mask);
        return T_(double(rng.next() & mask) * d);
      }
 
      static T_ gen_ranged(Random& rng, T_ a, T_ b)
      {
        // machine exactness; we do not want to divide by anything smaller than that...
        static const T_ eps = std::numeric_limits<T_>::epsilon();
 
        // call gen() prior to the if-clause to advance the rng in any case
        T_ x(gen(rng));
        if(a + eps < b)
          return a + x * (b - a);
        else
          return a;
      }
    };
 
    // specialization for bool
    template<typename T_>
    class RandomNumber<T_, Type::BooleanClass>
    {
    public:
      static T_ gen(Random& rng)
      {
        return T_(int(rng.next() & 0x8ull) != 0);
      }
 
      static T_ gen_ranged(Random& rng, T_ a, T_ b)
      {
        // call gen() prior to the if-clause to advance the rng in any case
        T_ x(gen(rng));
        // The following case may look pointless for bool, but if a and b are equal
        // we want to return that value here - otherwise its either true or false.
        if(a != b)
          return x;
        else
          return a;
      }
    };
  } // namespace Intern
} // namespace FEAT