random.h

// -*-c++-*-
#ifndef STORMM_RANDOM_H
#define STORMM_RANDOM_H
 
#include <vector>
#include "copyright.h"
#include "Accelerator/hybrid.h"
#include "Constants/scaling.h"
#include "Constants/hpc_bounds.h"
#include "DataTypes/common_types.h"
#include "DataTypes/stormm_vector_types.h"
#include "Math/rounding.h"
#include "Numerics/split_fixed_precision.h"
#include "random_enumerators.h"
 
namespace stormm {
namespace random {
 
using card::Hybrid;
using card::HybridKind;
using card::HybridTargetLevel;
using constants::giga_zu;
using data_types::isFloatingPointScalarType;
using stmath::roundUp;

constexpr int im1 = 2147483563;
constexpr int im2 = 2147483399;
constexpr double am = 1.0 / im1;
constexpr int im1_minus_one = im1 - 1;
constexpr int ia1 = 40014;
constexpr int ia2 = 40692;
constexpr int iq1 = 53668;
constexpr int iq2 = 52774;
constexpr int ir1 = 12211;
constexpr int ir2 = 3791;
constexpr int ntab = 32;
constexpr int ndiv = 1 + im1_minus_one / ntab;
constexpr double double_increment = 2.2205e-16;
constexpr double ran2_max = 1.0 - double_increment;

constexpr int default_random_seed = 827493;

constexpr int default_xoroshiro128p_scrub = 25;
constexpr int default_xoshiro256pp_scrub  = 50;
  
constexpr int max_xo_long_jumps = 1024;

constexpr double rng_unit_bin_offset  = 1.1102230246251565404236316680908203125e-16;
constexpr float rng_unit_bin_offset_f = 5.960464477539062500e-08f;
  
class Ran2Generator {
public:

  Ran2Generator(int igseed = default_random_seed);  

  double uniformRandomNumber();

  float spUniformRandomNumber();

  double gaussianRandomNumber();

  float spGaussianRandomNumber();
 
private:
  int iunstable_a;
  int iunstable_b;
  int state_sample;
  int state_vector[ntab + 8];
};

constexpr ullint xrs128p_jump_i      = 0xdf900294d8f554a5LLU;
constexpr ullint xrs128p_jump_ii     = 0x170865df4b3201fcLLU;
constexpr ullint xrs128p_longjump_i  = 0xd2a98b26625eee7bLLU;
constexpr ullint xrs128p_longjump_ii = 0xdddf9b1090aa7ac1LLU;

constexpr ullint xrs256pp_jump_i       = 0x180ec6d33cfd0abaLLU;
constexpr ullint xrs256pp_jump_ii      = 0xd5a61266f0c9392cLLU;
constexpr ullint xrs256pp_jump_iii     = 0xa9582618e03fc9aaLLU;
constexpr ullint xrs256pp_jump_iv      = 0x39abdc4529b1661cLLU;
constexpr ullint xrs256pp_longjump_i   = 0x76e15d3efefdcbbfLLU;
constexpr ullint xrs256pp_longjump_ii  = 0xc5004e441c522fb3LLU;
constexpr ullint xrs256pp_longjump_iii = 0x77710069854ee241LLU;
constexpr ullint xrs256pp_longjump_iv  = 0x39109bb02acbe635LLU;

class Xoroshiro128pGenerator {
public:

  Xoroshiro128pGenerator(int igseed = default_random_seed, int niter = 25);
  Xoroshiro128pGenerator(const ullint2 state_in);
  
  double uniformRandomNumber();

  float spUniformRandomNumber();

  double gaussianRandomNumber();

  float spGaussianRandomNumber();

  void jump();
  
  void longJump();

  ullint2 revealState() const;

  ullint revealBitString() const;

  void setState(const ullint2 state_in);
  
private:
  ullint2 state;  
  
  ullint2 seed128(int igseed);
  
  ullint next();

  void fastForward(const ullint2 stride);
};

class Xoshiro256ppGenerator {
public:

  Xoshiro256ppGenerator(int igseed = default_random_seed, int niter = 50);
  Xoshiro256ppGenerator(const ullint4 state_in);

  double uniformRandomNumber();

  float spUniformRandomNumber();

  double gaussianRandomNumber();

  float spGaussianRandomNumber();

  void jump();
  
  void longJump();

  ullint4 revealState() const;

  ullint revealBitString() const;

  void setState(const ullint4 state_in);
  
private:
  ullint4 state;  
  
  ullint4 seed256(int igseed);
  
  ullint next();

  void fastForward(const ullint4 stride);
};

template <typename T> class RandomNumberMill {
public:

  RandomNumberMill(const ullint2 state_in, size_t generators_in = 1LLU, size_t depth_in = 2LLU,
                   RandomNumberKind init_kind = RandomNumberKind::GAUSSIAN,
                   size_t bank_limit = constants::giga_zu);
 
  RandomNumberMill(const ullint4 state_in, size_t generators_in = 1LLU, size_t depth_in = 2LLU,
                   RandomNumberKind init_kind = RandomNumberKind::GAUSSIAN,
                   size_t bank_limit = constants::giga_zu);
 
  RandomNumberMill(size_t generators_in = 1LLU, size_t depth_in = 2LLU,
                   RandomAlgorithm style_in = RandomAlgorithm::XOSHIRO_256PP,
                   RandomNumberKind init_kind = RandomNumberKind::GAUSSIAN,
                   int igseed_in = default_random_seed, int niter = default_xoshiro256pp_scrub,
                   size_t bank_limit = constants::giga_zu);

  size_t getGeneratorCount() const;

  size_t getDepth() const;
  
  size_t getRefreshStride() const;

  T getBankValue(size_t generator_index, size_t layer_index) const;
  
  void uniformRandomNumbers(size_t first_gen, size_t last_gen);

  void gaussianRandomNumbers(size_t first_gen, size_t last_gen);
 
private:
  RandomAlgorithm style;     
  size_t generators;         
  size_t depth;              
  size_t refresh_stride;     
  Hybrid<ullint2> state_xy;  
  Hybrid<ullint2> state_zw;  
  Hybrid<T> bank;            

  void checkDimensions(size_t bank_limit);

  void initializeBank(RandomNumberKind init_kind);
};

template <typename Tprng> double uniformRand(Tprng *rng, double scale = 1.0);
 
template <typename Tprng> std::vector<double> uniformRand(Tprng *rng, size_t count,
                                                          double scale);
 
template <typename Tprng>
std::vector<double> uniformRand(Tprng *rng, size_t rows, size_t columns, double scale,
                                RngFillMode mode = RngFillMode::COLUMNS);
 
template <typename Tprng, typename Tprod>
void uniformRand(Tprng *rng, std::vector<Tprod> *xv, double scale = 1.0, double fp_scale = 1.0);

template <typename Tprng> float spUniformRand(Tprng *rng, float scale = 1.0);
 
template <typename Tprng> std::vector<float> spUniformRand(Tprng *rng, size_t count, float scale);
 
template <typename Tprng>
std::vector<float> spUniformRand(Tprng *rng, size_t rows, size_t columns, float scale,
                                 RngFillMode mode = RngFillMode::COLUMNS);
 
template <typename Tprng, typename Tprod>
void spUniformRand(Tprng *rng, std::vector<Tprod> *xv, float scale = 1.0, float fp_scale = 1.0);

template <typename Tprng> double guassianRand(Tprng *rng, double scale = 1.0);
 
template <typename Tprng> std::vector<double> guassianRand(Tprng *rng, size_t count, double scale);
 
template <typename Tprng>
std::vector<double> guassianRand(Tprng *rng, size_t rows, size_t columns, double scale,
                                 RngFillMode mode = RngFillMode::COLUMNS);
 
template <typename Tprng, typename Tprod>
void guassianRand(Tprng *rng, std::vector<Tprod> *xv, double fp_scale = 1.0, double scale = 1.0);

template <typename Tprng> float spGuassianRand(Tprng *rng, float scale = 1.0);
 
template <typename Tprng> std::vector<float> spGuassianRand(Tprng *rng, size_t count, float scale);
 
template <typename Tprng>
std::vector<float> spGuassianRand(Tprng *rng, size_t rows, size_t columns, float scale,
                                  RngFillMode mode = RngFillMode::COLUMNS);
 
template <typename Tprng, typename Tprod>
void spGuassianRand(Tprng *rng, std::vector<Tprod> *xv, float fp_scale = 1.0, float scale = 1.0);

void initXoroshiro128pArray(ullint2* state_vector, int n_generators, int igseed,
                            int scrub_cycles = default_xoroshiro128p_scrub);
void initXoroshiro128pArray(std::vector<ullint2> *state_vector, int igseed,
                            int scrub_cycles = default_xoroshiro128p_scrub);
void initXoroshiro128pArray(Hybrid<ullint2> *state_vector, int igseed,
                            int scrub_cycles = default_xoroshiro128p_scrub);

void initXoshiro256ppArray(ullint2* state_xy, ullint2* state_zw, int n_generators, int igseed,
                           int scrub_cycles = default_xoshiro256pp_scrub);
void initXoshiro256ppArray(std::vector<ullint2> *state_xy, std::vector<ullint2> *state_zw,
                           int igseed, int scrub_cycles = default_xoshiro256pp_scrub);
void initXoshiro256ppArray(Hybrid<ullint2> *state_xy, Hybrid<ullint2> *state_zw, int igseed,
                           int scrub_cycles = default_xoshiro256pp_scrub);

template <typename T>
void fillRandomCache(ullint2* state_xy, ullint2* state_zw, T* cache, size_t length, size_t depth,
                     RandomAlgorithm method, RandomNumberKind product, size_t index_start,
                     size_t index_end);
 
template <typename T>
void fillRandomCache(std::vector<ullint2> *state_xy, std::vector<ullint2> *state_zw,
                     std::vector<T> *cache, size_t length, size_t depth, RandomAlgorithm method,
                     RandomNumberKind product, size_t index_start, size_t index_end);
 
template <typename T>
void fillRandomCache(Hybrid<ullint2> *state_xy, Hybrid<ullint2> *state_zw, std::vector<T> *cache,
                     size_t length, size_t depth, RandomAlgorithm method, RandomNumberKind product,
                     size_t index_start, size_t index_end);

template <typename Trng, typename Tvar>
void addRandomNoise(Trng *prng, Tvar* x, size_t length, double mult, double scale = 1.0,
                    RandomNumberKind kind = RandomNumberKind::GAUSSIAN);
 
template <typename Trng, typename Tvar>
void addRandomNoise(Trng *prng, std::vector<Tvar> *x, double mult, double scale = 1.0,
                    RandomNumberKind kind = RandomNumberKind::GAUSSIAN);
 
template <typename Trng, typename Tvar>
void addRandomNoise(Trng *prng, Hybrid<Tvar> *x, double mult, double scale = 1.0,
                    RandomNumberKind kind = RandomNumberKind::GAUSSIAN);
 
template <typename Trng, typename Tvar>
void addRandomNoise(Trng *prng, Tvar* x, Tvar* y, Tvar *z, size_t length, double mult,
                    double scale = 1.0, RandomNumberKind kind = RandomNumberKind::GAUSSIAN);
 
template <typename Trng, typename Tvar>
void addRandomNoise(Trng *prng, std::vector<Tvar> *x, std::vector<Tvar> *y, std::vector<Tvar> *z,
                    double mult, double scale = 1.0,
                    RandomNumberKind kind = RandomNumberKind::GAUSSIAN);
 
template <typename Trng, typename Tvar>
void addRandomNoise(Trng *prng, Hybrid<Tvar> *x, Hybrid<Tvar> *y, Hybrid<Tvar> *z, double mult,
                    double scale = 1.0, RandomNumberKind kind = RandomNumberKind::GAUSSIAN);
 
template <typename Trng>
void addRandomNoise(Trng *prng, llint* x, int* x_ovrf, size_t length, double mult, double scale,
                    RandomNumberKind kind = RandomNumberKind::GAUSSIAN);
 
template <typename Trng>
void addRandomNoise(Trng *prng, std::vector<llint> *x, std::vector<int> *x_ovrf, double mult,
                    double scale, RandomNumberKind kind = RandomNumberKind::GAUSSIAN);
 
template <typename Trng>
void addRandomNoise(Trng *prng, Hybrid<llint> *x, Hybrid<int> *x_ovrf, double mult, double scale,
                    RandomNumberKind kind = RandomNumberKind::GAUSSIAN);
 
template <typename Trng>
void addRandomNoise(Trng *prng, llint* x, int* x_ovrf, llint* y, int* y_ovrf, llint* z,
                    int* z_ovrf, size_t length, double mult, double scale,
                    RandomNumberKind kind = RandomNumberKind::GAUSSIAN);
 
template <typename Trng>
void addRandomNoise(Trng *prng, std::vector<llint> *x, std::vector<int> *x_ovrf,
                    std::vector<llint> *y, std::vector<int> *y_ovrf, std::vector<llint> *z,
                    std::vector<int> *z_ovrf, double mult, double scale,
                    RandomNumberKind kind = RandomNumberKind::GAUSSIAN);
 
template <typename Trng>
void addRandomNoise(Trng *prng, Hybrid<llint> *x, Hybrid<int> *x_ovrf, Hybrid<llint> *y,
                    Hybrid<int> *y_ovrf, Hybrid<llint> *z, Hybrid<int> *z_ovrf, double mult,
                    double scale, RandomNumberKind kind = RandomNumberKind::GAUSSIAN);
  
} // namespace random
} // namespace stormm
 
#include "random.tpp"
 
#endif