the PRNG is from code here: http://prng.di.unimi.it ; information on how to code URNG template class for STL's random header is here https://stackoverflow.com/questions/25129500/using-stdshuffle-with-custom-rng and https://www.boost.org/doc/libs/1_72_0/doc/html/boost_random/reference.html ; useful discussion is found on the dev's boards https://…

xorshift.cpp

#include <random>

#include <array>
#include <iostream>
#include <limits>
#include <random>

#include <Rcpp.h>

// [[Rcpp::plugins(cpp14)]]

// for details on what a URNG compatible with STL must have
// see: https://www.boost.org/doc/libs/1_72_0/doc/html/boost_random/reference.html 
constexpr std::uint64_t rotl(std::uint64_t x, int k) noexcept
{
  return (x << k) | (x >> (64 - k));
}

class xoshiro256ss
{
public:
  using result_type = std::uint64_t;
  
  // min value that the sequence will produce
  constexpr static result_type min() noexcept { return 0; }
  
  // max value that the sequence will produce
  constexpr static result_type max() noexcept
  { return std::numeric_limits<result_type>::max(); }
  
  // gives value in [min,max]
  result_type operator()() noexcept
  {
    const auto result_starstar(rotl(state[1] * 5, 7) * 9);
    const auto t(state[1] << 17);
    
    state[2] ^= state[0];
    state[3] ^= state[1];
    state[1] ^= state[2];
    state[0] ^= state[3];
    
    state[2] ^= t;
    
    state[3] = rotl(state[3], 45);
    
    return result_starstar;
  }
  
  // set seed
  void seed(const std::array<std::uint64_t, 4> &s) noexcept {
    state = s;
  };
  
  bool operator==(const xoshiro256ss &rhs) const noexcept
  { return state == rhs.state; }
  bool operator!=(const xoshiro256ss &rhs) const noexcept
  { return !(*this == rhs); }
  
  friend std::ostream &operator<<(std::ostream &, const xoshiro256ss &);
  friend std::istream &operator>>(std::istream &, xoshiro256ss &);
  
private:  
  std::array<std::uint64_t, 4> state;
};

// [[Rcpp::export]]
Rcpp::NumericVector test_prng(
  const std::uint64_t s1,
  const std::uint64_t s2,
  const std::uint64_t s3,
  const std::uint64_t s4,
  const int n
){
  
  std::array<std::uint64_t, 4> seeds = {s1,s2,s3,s4};
  xoshiro256ss gen;
  gen.seed(seeds);
  
  Rcpp::NumericVector out(n);
  
  std::uniform_real_distribution<double> distribution(0.0,1.0);
  
  for(int i=0; i<n; i++){
    out[i] = distribution(gen);
  }
  
  return out;
}

/*** R
x <- test_prng(s1 = 4534534L,s2 = 312L,s3 = 90956L,s4 = 342L,n = 1e4)
*/