Fast Pseudo Random Number Generators for R

// [[Rcpp::depends(dqrng)]]
#include <Rcpp.h>
#include <dqrng.h>

using Rcpp::IntegerVector;
using Rcpp::NumericVector;
using Rcpp::sqrt;
using Rcpp::sum;
using dqrng::dqrunif;

// [[Rcpp::export]]
double piCpp(const int n) {
  dqrng::dqset_seed(IntegerVector::create(42));
  NumericVector x = dqrunif(n);
  NumericVector y = dqrunif(n);
  NumericVector d = sqrt(x*x + y*y);
  return 4.0 * sum(d < 1.0) / n;
}
/*** R
system.time(cat("pi ~= ", piCpp(1e7), "\n"))
*/

// [[Rcpp::depends(dqrng)]]
#include <Rcpp.h>
#include <dqrng.h>

// [[Rcpp::export]]
void sampleCpp(const int n) {
  dqrng::dqset_seed(Rcpp::IntegerVector::create(42));
  Rcpp::IntegerVector sample = dqrng::dqsample_int(n, n/100, true);
  Rcpp::Rcout << sample << std::endl;
}
/*** R
sampleCpp(1000)
*/

#include <Rcpp.h>
// [[Rcpp::depends(dqrng, BH)]]
#include <dqrng_distribution.h>

class SplitMix {
public:
  typedef uint64_t result_type;
  SplitMix (result_type seed) : state(seed) {};
  result_type operator() () {
    result_type z = (state += 0x9e3779b97f4a7c15ULL);
    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL;
    z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL;
    return z ^ (z >> 31);
  }
  void seed(result_type seed) {state = seed;}
  static constexpr result_type min() {return 0;};
  static constexpr result_type max() {return UINT64_MAX;};

private:
  result_type state;

public:  
  friend std::ostream& operator<<(std::ostream& ost, const SplitMix& e) {
    return ost << e.state;
  }
  friend std::istream& operator>>(std::istream& ist, SplitMix& e) {
    return ist >> e.state;
  }
};

// [[Rcpp::export]]
Rcpp::NumericVector splitmix_rnorm(const int n, const double mean = 0.0, const double sd = 1.0) {
  auto rng = dqrng::generator<SplitMix>(42);
  Rcpp::NumericVector result(n);
  rng->generate<dqrng::normal_distribution>(result, mean, sd);
  return result;
}
/*** R
splitmix_rnorm(10)
system.time(splitmix_rnorm(1e7))
*/

#include <Rcpp.h>
// [[Rcpp::depends(dqrng, BH, sitmo)]]
#include <dqrng_distribution.h>
#include <threefry.h>

// [[Rcpp::export]]
Rcpp::NumericVector threefry_rnorm(const int n, const double mean = 0.0, const double sd = 1.0) {
  auto rng = dqrng::generator<sitmo::threefry_13_64>(42);
  Rcpp::NumericVector result(n);
  rng->generate<dqrng::normal_distribution>(result, mean, sd);
  return result;
}

/*** R
threefry_rnorm(10)
system.time(threefry_rnorm(1e7))
*/

#include <random>
#include <Rcpp.h>
// [[Rcpp::depends(dqrng)]]
#include <dqrng_generator.h>
#include <xoshiro.h>

// [[Rcpp::export]]
Rcpp::NumericVector std_rnorm(const int n, const double mean = 0.0, const double sd = 1.0) {
  auto rng = dqrng::generator<dqrng::xoroshiro128plusplus>(42);
  Rcpp::NumericVector result(n);
  rng->generate<std::normal_distribution>(result, mean, sd);
  return result;
}

/*** R
std_rnorm(10)
system.time(std_rnorm(1e7))
*/

#include <Rcpp.h>
// [[Rcpp::depends(dqrng, BH)]]
#include <pcg_random.hpp>
#include <boost/random/normal_distribution.hpp>

// [[Rcpp::plugins(cpp11)]]   
// [[Rcpp::export]]
Rcpp::NumericVector boost_pcg_rnorm(const int n, const double mean = 0.0, const double sd = 1.0) {
  pcg32 rng(42);
  boost::random::normal_distribution<double> dist(mean, sd);
  Rcpp::NumericVector result(n);
  std::generate(result.begin(), result.end(), [&dist, &rng](){return dist(rng);});
  return result;
}
/*** R
boost_pcg_rnorm(10)
system.time(boost_pcg_rnorm(1e7))
*/

#include <Rcpp.h>
// [[Rcpp::depends(dqrng, BH)]]
#include <boost/random/binomial_distribution.hpp>
#include <dqrng_distribution.h>

// [[Rcpp::export]]
Rcpp::NumericMatrix multiple_distributions(int n) {
  auto rng = dqrng::generator<dqrng::xoshiro256plusplus>(42);
  Rcpp::NumericMatrix out(n, 3);
  double p = 0.0;
  for (int i = 0; i < n; ++i) {
    p = double(i) / double(n);
    out(i,0) = rng->variate<boost::random::binomial_distribution<int>>(1, p);
    out(i,1) = rng->variate<dqrng::normal_distribution>(p, 1.0);
    out(i,2) = rng->variate<dqrng::normal_distribution>(4.0, 3.0 - p);
  }
  Rcpp::colnames(out) = Rcpp::CharacterVector::create("Bernoulli", "Normal1", "Normal2");
  return out;
}

/*** R
multiple_distributions(5)
*/