Add clone(stream) method for the RNGs

.Rbuildignore

@@ -18,3 +18,5 @@
 ^.vscode$
 ^Meta$
 ^doc$
+^\.cache$
+^compile_commands\.json$

.gitignore

@@ -258,3 +258,5 @@ $RECYCLE.BIN/
 inst/doc
 vignettes/*.R
 /external/tmp
+.cache
+compile_commands.json

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: dqrng
 Type: Package
 Title: Fast Pseudo Random Number Generators
-Version: 0.3.2.2
+Version: 0.3.2.4
 Authors@R: c(
     person("Ralf", "Stubner", email = "ralf.stubner@gmail.com", role = c("aut", "cre")),
     person("daqana GmbH", role = "cph"),

NEWS.md

@@ -4,6 +4,7 @@
 
 * The default RNG has changed from Xoroshiro128+ to Xoroshiro128++. The older generators Xoroshiro128+ and Xoshiro256+ are still available but should only be used for backward compatibility or for generating floating point numbers, i.e. not sampling etc.  ([#57](https://github.com/daqana/dqrng/pull/57) fixing [#56](https://github.com/daqana/dqrng/issues/56))
 * The `dqrng::rng64_t` type has been changed to use `Rcpp::XPtr` instead of `std::shared_ptr` and the functions from `dqrng_sample.h` now expect a reference to `dqrng::random_64bit_generator` instead of `dqrng::rng64_t` ([#70](https://github.com/daqana/dqrng/pull/70) fixing [#63](https://github.com/daqana/dqrng/issues/63))
+* The two argument constructor and `seed`  function from PCG has [surprising properties](https://github.com/imneme/pcg-cpp/issues/91): it is not identical to the one argument version followed by `set_stream(stream)`. For consistency with the new `clone(stream)` method, the two argument versions are no longer used. This influences code the uses multiple stream with PCG together with the tooling from this package, e.g. the example code in the vignette on parallel RNG usage. 
 
 ## Other changes
 
@@ -13,6 +14,9 @@
 * Allow uniform and normal distributions to be registered as user-supplied RNG within R. This happens automatically if the option `dqrng.register_methods` is set to `TRUE`.
 * Add missing inline attributes and limit the included Rcpp headers in `dqrng_types.h` ([#75](https://github.com/daqana/dqrng/pull/75) together with Paul Liétar)
 * Add I/O methods for the RNG's internal state (fixing [#66](https://github.com/daqana/dqrng/issues/66) in [#78](https://github.com/daqana/dqrng/pull/78))
+* Extend `random_64bit_generator` with additional convenience methods (fixing [#64](https://github.com/daqana/dqrng/issues/64) in [#79](https://github.com/daqana/dqrng/pull/79))
+    * A `clone(stream)` method to allow using the global RNG state for parallel computation
+    * New methods `variate<dist>(param)`, `generate<dist>(container, param)` etc. using and inspired by [`randutils`](https://www.pcg-random.org/posts/ease-of-use-without-loss-of-power.html).
 
 
 # dgrng 0.3.2

README.md

@@ -61,8 +61,8 @@ bm[, 1:4]
 #> # A tibble: 2 × 4
 #>   expression      min   median `itr/sec`
 #>   <bch:expr> <bch:tm> <bch:tm>     <dbl>
-#> 1 rnorm(N)      610µs    709µs     1322.
-#> 2 dqrnorm(N)    197µs    222µs     4411.
+#> 1 rnorm(N)    608.9µs  692.8µs     1342.
+#> 2 dqrnorm(N)   82.7µs   86.2µs    10227.
 ```
 
 This is also true for the provided sampling functions with replacement:
@@ -79,10 +79,10 @@ bm[, 1:4]
 #> # A tibble: 4 × 4
 #>   expression                                     min   median `itr/sec`
 #>   <bch:expr>                                <bch:tm> <bch:tm>     <dbl>
-#> 1 sample.int(m, n, replace = TRUE)             7.2ms    7.7ms      120.
-#> 2 sample.int(1000 * m, n, replace = TRUE)      8.8ms   9.84ms      100.
-#> 3 dqsample.int(m, n, replace = TRUE)         435.2µs 497.26µs     1896.
-#> 4 dqsample.int(1000 * m, n, replace = TRUE)    446µs 496.73µs     1780.
+#> 1 sample.int(m, n, replace = TRUE)            7.03ms   7.46ms      133.
+#> 2 sample.int(1000 * m, n, replace = TRUE)     8.87ms    9.4ms      104.
+#> 3 dqsample.int(m, n, replace = TRUE)        292.69µs 316.75µs     2821.
+#> 4 dqsample.int(1000 * m, n, replace = TRUE) 408.94µs 446.19µs     1987.
 ```
 
 And without replacement:
@@ -100,11 +100,11 @@ bm[, 1:4]
 #> # A tibble: 5 × 4
 #>   expression                            min   median `itr/sec`
 #>   <bch:expr>                       <bch:tm> <bch:tm>     <dbl>
-#> 1 sample.int(m, n)                  42.52ms  43.95ms      16.3
-#> 2 sample.int(1000 * m, n)           13.35ms   14.8ms      65.8
-#> 3 sample.int(m, n, useHash = TRUE)   9.72ms  11.34ms      85.8
-#> 4 dqsample.int(m, n)               717.13µs 915.69µs     846. 
-#> 5 dqsample.int(1000 * m, n)          1.47ms   2.13ms     384.
+#> 1 sample.int(m, n)                  25.43ms   27.8ms      31.8
+#> 2 sample.int(1000 * m, n)           12.99ms   17.4ms      57.0
+#> 3 sample.int(m, n, useHash = TRUE)  10.44ms   12.8ms      77.2
+#> 4 dqsample.int(m, n)               617.83µs  762.9µs     882. 
+#> 5 dqsample.int(1000 * m, n)          1.51ms    2.5ms     334.
 ```
 
 Note that sampling from `10^10` elements triggers “long-vector support”

inst/include/dqrng.h

@@ -22,6 +22,6 @@
 #include "dqrng_RcppExports.h"
 
 namespace dqrng {
-random_64bit_accessor::random_64bit_accessor() : gen(dqrng::get_rng()) {}
+inline random_64bit_accessor::random_64bit_accessor() : gen(dqrng::get_rng()) {}
 } // namespace dqrng
 #endif // dqrng_H

inst/include/dqrng_distribution.h

@@ -1,4 +1,5 @@
 // Copyright 2018-2019 Ralf Stubner (daqana GmbH)
+// Copyright 2023-2024 Ralf Stubner
 //
 // This file is part of dqrng.
 //
@@ -45,17 +46,13 @@ namespace detail {
 template<>
 inline std::pair<double, int> generate_int_float_pair<double, 8, dqrng::random_64bit_generator>(dqrng::random_64bit_generator& eng)
 {
-  dqrng::random_64bit_generator::result_type x = eng();
-  double r = dqrng::uniform01(x);
-  // shift x due to weakness of lowest bits for xoshiro/xoroshiro with used "+" scrambler
-  int bucket = (x >> 3) & 0xFF;
-  return std::make_pair(r, bucket);
+  return eng.generate_double_8bit_pair();
 }
 
 template<>
 inline double generate_uniform_real<dqrng::random_64bit_generator, double>(dqrng::random_64bit_generator& eng, double min, double max)
 {
-  return dqrng::uniform01(eng()) * (max - min) + min;
+  return eng.uniform01() * (max - min) + min;
 }
 
 template<>

inst/include/dqrng_generator.h

@@ -1,5 +1,5 @@
 // Copyright 2018-2019 Ralf Stubner (daqana GmbH)
-// Copyright 2023 Ralf Stubner
+// Copyright 2023-2024 Ralf Stubner
 //
 // This file is part of dqrng.
 //
@@ -26,6 +26,12 @@
 #include <pcg_random.hpp>
 #include <Rcpp.h>
 
+#if defined(__cpp_lib_make_unique) && (__cpp_lib_make_unique >= 201304)
+using std::make_unique;
+#else
+#include <boost/smart_ptr/make_unique.hpp>
+using boost::make_unique;
+#endif
 
 namespace dqrng {
 using rng64_t = Rcpp::XPtr<random_64bit_generator>;
@@ -39,155 +45,62 @@ class random_64bit_wrapper : public random_64bit_generator {
   static_assert(RNG::min() == 0, "Provided RNG has wrong minimum.");
 private:
   RNG gen;
-  bool has_cache{false};
-  uint32_t cache;
-
-  uint64_t bit64() {return gen();}
-  uint32_t bit32() {
-    if (has_cache) {
-      has_cache = false;
-      return cache;
-    }
-    uint64_t random = this->bit64();
-    cache = uint32_t(random);
-    has_cache = true;
-    return random >> 32;
-  }
 
 protected:
   virtual void output(std::ostream& ost) const override {ost << gen;}
   virtual void input(std::istream& ist) override {ist >> gen;}
+  void set_stream(result_type stream) {throw std::runtime_error("Stream handling not supported for this RNG!");}
 
 public:
   random_64bit_wrapper() : gen() {};
+  random_64bit_wrapper(RNG _gen) : gen(_gen) {};
   random_64bit_wrapper(result_type seed) : gen(seed) {};
-  random_64bit_wrapper(result_type seed, result_type stream) : gen(seed, stream) {};
-  virtual result_type operator() () override {return this->bit64();}
+  random_64bit_wrapper(result_type seed, result_type stream) : gen(seed) {this->set_stream(stream);};
+  virtual result_type operator() () override {return gen();}
   virtual void seed(result_type seed) override {cache = false; gen.seed(seed);}
-  virtual void seed(result_type seed, result_type stream) override {throw std::runtime_error("Stream handling not supported for this RNG!");}
-
-  /*
-   * https://raw.githubusercontent.com/imneme/bounded-rands/3d71f53c975b1e5b29f2f3b05a74e26dab9c3d84/bounded32.cpp
-   * https://raw.githubusercontent.com/imneme/bounded-rands/3d71f53c975b1e5b29f2f3b05a74e26dab9c3d84/bounded64.cpp
-   * A C++ implementation methods and benchmarks for random numbers in a range
-   * (64 and 32-bit version)
-   *
-   * The MIT License (MIT)
-   *
-   * Copyright (c) 2018 Melissa E. O'Neill
-   *
-   * Permission is hereby granted, free of charge, to any person obtaining a
-   * copy of this software and associated documentation files (the "Software"),
-   * to deal in the Software without restriction, including without limitation
-   * the rights to use, copy, modify, merge, publish, distribute, sublicense,
-   * and/or sell copies of the Software, and to permit persons to whom the
-   * Software is furnished to do so, subject to the following conditions:
-   *
-   * The above copyright notice and this permission notice shall be included in
-   * all copies or substantial portions of the Software.
-   *
-   * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-   * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-   * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-   * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-   * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
-   * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
-   * DEALINGS IN THE SOFTWARE.
-   */
-
-  virtual uint32_t operator() (uint32_t range) override {
-    uint32_t x = this->bit32();
-    uint64_t m = uint64_t(x) * uint64_t(range);
-    uint32_t l = uint32_t(m);
-    if (l < range) {
-      uint32_t t = -range;
-      if (t >= range) {
-        t -= range;
-        if (t >= range)
-          t %= range;
-      }
-      while (l < t) {
-        x = this->bit32();
-        m = uint64_t(x) * uint64_t(range);
-        l = uint32_t(m);
-      }
-    }
-    return m >> 32;
-  }
-
-#ifdef LONG_VECTOR_SUPPORT
-  virtual uint64_t operator() (uint64_t range) override {
-    using pcg_extras::pcg128_t;
-    uint64_t x = this->bit64();
-    pcg128_t m = pcg128_t(x) * pcg128_t(range);
-    uint64_t l = uint64_t(m);
-    if (l < range) { // # nocov start
-      uint64_t t = -range;
-      if (t >= range) {
-        t -= range;
-        if (t >= range)
-          t %= range;
-      }
-      while (l < t) {
-        x = this->bit64();
-        m = pcg128_t(x) * pcg128_t(range);
-        l = uint64_t(m);
-      }
-    } // # nocov end
-    return m >> 64;
+  virtual void seed(result_type seed, result_type stream) override {cache = false; gen.seed(seed); this->set_stream(stream);}
+  virtual std::unique_ptr<random_64bit_generator> clone(result_type stream) override {
+    auto rng = make_unique<random_64bit_wrapper<RNG>>(gen);
+    rng->set_stream(stream);
+    return rng;
   }
-#endif
 };
 
 template<>
-inline void random_64bit_wrapper<::dqrng::xoroshiro128plus>::seed(result_type seed, result_type stream) {
-    gen.seed(seed);
-    gen.jump(stream);
-    cache = false;
+inline void random_64bit_wrapper<::dqrng::xoroshiro128plus>::set_stream(result_type stream) {
+  gen.jump(stream);
 }
 
 template<>
-inline void random_64bit_wrapper<::dqrng::xoroshiro128plusplus>::seed(result_type seed, result_type stream) {
-  gen.seed(seed);
+inline void random_64bit_wrapper<::dqrng::xoroshiro128plusplus>::set_stream(result_type stream) {
   gen.jump(stream);
-  cache = false;
 }
 
 template<>
-inline void random_64bit_wrapper<::dqrng::xoroshiro128starstar>::seed(result_type seed, result_type stream) {
-  gen.seed(seed);
+inline void random_64bit_wrapper<::dqrng::xoroshiro128starstar>::set_stream(result_type stream) {
   gen.jump(stream);
-  cache = false;
 }
 
 template<>
-inline void random_64bit_wrapper<::dqrng::xoshiro256plus>::seed(result_type seed, result_type stream) {
-  gen.seed(seed);
+inline void random_64bit_wrapper<::dqrng::xoshiro256plus>::set_stream(result_type stream) {
   gen.long_jump(stream);
-  cache = false;
 }
 
 template<>
-inline void random_64bit_wrapper<::dqrng::xoshiro256plusplus>::seed(result_type seed, result_type stream) {
-  gen.seed(seed);
+inline void random_64bit_wrapper<::dqrng::xoshiro256plusplus>::set_stream(result_type stream) {
   gen.long_jump(stream);
-  cache = false;
 }
 
 template<>
-inline void random_64bit_wrapper<::dqrng::xoshiro256starstar>::seed(result_type seed, result_type stream) {
-    gen.seed(seed);
-    gen.long_jump(stream);
-    cache = false;
+inline void random_64bit_wrapper<::dqrng::xoshiro256starstar>::set_stream(result_type stream) {
+  gen.long_jump(stream);
 }
 
 template<>
-inline void random_64bit_wrapper<pcg64>::seed(result_type seed, result_type stream) {
-    gen.seed(seed, stream);
-    cache = false;
+inline void random_64bit_wrapper<pcg64>::set_stream(result_type stream) {
+  gen.set_stream(stream);
 }
 
-
 template<typename RNG = default_64bit_generator>
 typename std::enable_if<!std::is_base_of<random_64bit_generator, RNG>::value, rng64_t>::type
 generator () {

inst/include/dqrng_threefry.h

@@ -1,5 +1,5 @@
 // Copyright 2018-2019 Ralf Stubner (daqana GmbH)
-// Copyright 2023 Ralf Stubner
+// Copyright 2023-2024 Ralf Stubner
 //
 // This file is part of dqrng.
 //
@@ -19,16 +19,23 @@
 #ifndef DQRNG_THREEFRY_H
 #define DQRNG_THREEFRY_H 1
 
+#include <sstream>
+#include <vector>
 #include <threefry.h>
 #include <dqrng_generator.h>
 
 namespace dqrng {
 
 template<>
-inline void random_64bit_wrapper<sitmo::threefry_20_64>::seed(result_type seed, result_type stream) {
-  gen.seed(seed);
-  gen.set_counter(0, 0, 0, stream);
-  cache = false;
+inline void random_64bit_wrapper<sitmo::threefry_20_64>::set_stream(result_type stream) {
+  uint64_t number;
+  std::vector<uint64_t> state;
+  std::stringstream iss;
+  iss << gen;
+  while (iss >> number)
+    state.push_back(number);
+  // state[4:7] is the current counter, the highest part is incremented by stream
+  gen.set_counter(state[4], state[5], state[6], state[7] + stream);
 }
 } // namespace dqrng