c++ - 在 Rcpp 中有效地生成随机比特流

Question

我目前正在构建的 R 包中有一个名为 rbinom01 的辅助函数。请注意，它调用 random(3)。

int rbinom01(int size) {
  if (!size) {
    return 0;
  }

  int64_t result = 0;
  while (size >= 32) {
    result += __builtin_popcount(random());
    size -= 32;
  }

  result += __builtin_popcount(random() & ~(LONG_MAX << size));

  return result;
}

当 R CMD 检查 my_package 时，我收到以下警告:

* checking compiled code ... NOTE
File ‘ my_package/libs/my_package.so’:
  Found ‘_random’, possibly from ‘random’ (C)
    Object: ‘ my_function.o’

Compiled code should not call entry points which might terminate R nor
write to stdout/stderr instead of to the console, nor use Fortran I/O
nor system RNGs.

See ‘Writing portable packages’ in the ‘Writing R Extensions’ manual.

我前往 the Document ，它说我可以使用 *_rand 函数之一，以及一系列 distribution functions .嗯，这很酷，但我的包只需要一个随机位流而不是一个随机的 double。最简单的方法是使用 random(3) 或者从 /dev/urandom 中读取，但这会使我的包“不可移植”。

This post建议使用 sample，但遗憾的是它不适合我的用例。对于我的应用程序，生成随机位显然对性能至关重要，因此我不希望它浪费任何时间调用 unif_rand，将结果乘以 N 并四舍五入。不管怎样，我使用 C++ 的原因是为了利用位级并行性。

当然，我可以手动滚动自己的 PRNG 或复制并粘贴最先进的 PRNG 的代码，如 xoshiro256** ，但在此之前，我想看看是否有更简单的选择。

顺便说一下，有人可以给我一个不错的 Rcpp 简短教程链接吗？ Writing R Extensions 非常全面而且很棒，但我需要数周时间才能完成。我正在寻找一个更简洁的版本，但最好它应该比调用 Rcpp.package.skeleton 提供更多信息。

---

根据 @Ralf Stubner 的建议的答案，我已经重新编写了原始代码如下。但是，我每次都得到相同的结果。我怎样才能正确地播种它并同时保持我的代码“可移植”？

int rbinom01(int size) {
  dqrng::xoshiro256plus rng;

  if (!size) {
    return 0;
  }

  int result = 0;
  while (size >= 64) {
    result += __builtin_popcountll(rng());
    Rcout << sizeof(rng()) << std::endl;
    size -= 64;
  }

  result += __builtin_popcountll(rng() & ((1LLU << size) - 1));

  return result;
}

Answer 1

有不同的 R 包使 PRNG 可用作仅 C++ header 库:

• BH : 来自 boost.random
的一切
• sitmo : 各种 Threefry 版本
• dqrng : PCG 家族，xoshiro256+ 和 xoroshiro128+
• ...

您可以通过将 LinkingTo 添加到您的包的 DECRIPTION 来使用其中的任何一个。通常，这些 PRNG 是在 C++11 random header 之后建模的，这意味着您必须控制它们的生命周期并自行播种。在单线程环境中，我喜欢使用匿名命名空间来进行生命周期控制，例如:

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

namespace {
dqrng::xoshiro256plus rng{};
}

// [[Rcpp::export]]
void set_seed(int seed) {
  rng.seed(seed);
}

// [[Rcpp::export]]
int rbinom01(int size) {
  if (!size) {
    return 0;
  }

  int result = 0;
  while (size >= 64) {
    result += __builtin_popcountll(rng());
    size -= 64;
  }

  result += __builtin_popcountll(rng() & ((1LLU << size) - 1));

  return result;
}

/*** R
set_seed(42)
rbinom01(10)
rbinom01(10)
rbinom01(10)
*/

但是，使用 runif 并不全是坏事，而且肯定比访问 /dev/urandom 更快。在 dqrng 中有一个 convenient wrapper为此。

至于教程:除了 WRE Rcpp package小插图是必读的。 R Packages如果您想使用 devtools-way，Hadley Wickham 也有一章是关于“编译代码”的。