randoms.hpp

Go to the documentation of this file.

// vim: set tabstop=2:softtabstop=2:shiftwidth=2:expandtab

#ifndef LIN_GENERATORS_RANDOMS_HPP_
#define LIN_GENERATORS_RANDOMS_HPP_

#include "../core.hpp"

#include <cstdint>
#include <type_traits>

namespace lin {
namespace internal {

class RandomsGenerator {
 private:
  unsigned long long seed;

  double cached_rand;

  bool has_cached;

 public:
  constexpr RandomsGenerator(RandomsGenerator const &) = default;
  constexpr RandomsGenerator(RandomsGenerator &&) = default;
  constexpr RandomsGenerator &operator=(RandomsGenerator const &) = default;
  constexpr RandomsGenerator &operator=(RandomsGenerator &&) = default;

  constexpr RandomsGenerator(unsigned long long seed = 0) : 
    seed(seed ^ 4101842887655102017LL), 
    cached_rand(0.0),
    has_cached(false) { }

  constexpr double rand() {
    seed ^= (seed >> 21);
    seed ^= (seed << 35);
    seed ^= (seed >> 4);

    return 5.42101086242752217E-20 * (seed * 2685821657736338717ULL);
  }

  constexpr double gaussian() {
    if(has_cached) {
      has_cached = false;
      return cached_rand;
    }
    else{
      double R = std::sqrt(-2.0*std::log(rand()));
      double T = 2.0*M_PI*rand();
      cached_rand = R*std::sin(T);
      has_cached = true;
      return R*std::cos(T);
    }
  }
};
}  // namespace internal

template <class C, std::enable_if_t<internal::has_traits<C>::value, size_t> = 0>
constexpr auto rands(internal::RandomsGenerator &rand, size_t r = C::Traits::max_rows, size_t c = C::Traits::max_cols) {
  typename C::Traits::eval_t t(r, c);
  for (lin::size_t i = 0; i < t.size(); i++) t(i) = typename C::Traits::elem_t(rand.rand());
  return t;
}

template <class C, std::enable_if_t<internal::has_traits<C>::value, size_t> = 0>
constexpr auto gaussians(internal::RandomsGenerator &rand, size_t r = C::Traits::max_rows, size_t c = C::Traits::max_cols) {
  typename C::Traits::eval_t t(r, c);
  for (lin::size_t i = 0; i < t.size(); i++) t(i) = typename C::Traits::elem_t(rand.gaussian());
  return t;
}
}  // namespace lin

#endif