在 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 check 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.

我前往文档，它说我可以使用*_rand函数之一，以及一系列分发函数。嗯，这很酷，但我的包只需要一个随机位流而不是一个随机double。我可以拥有它的最简单的方法是使用random(3)或从/dev/urandom读取，但这会使我的包"不可移植"。

这篇文章建议使用sample，但不幸的是它不适合我的用例。对于我的应用程序，生成随机位显然对性能至关重要，所以我不希望它浪费任何时间调用unif_rand，将结果乘以N并对其进行四舍五入。无论如何，我使用 C++ 的原因是为了利用位级并行性。

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

顺便说一句，有人可以将 Rcpp 的一个很好的简短教程链接到我吗？编写 R 扩展是全面且很棒的，但我需要几周才能完成。我正在寻找一个更简洁的版本，但最好它应该比打电话给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++标头库：

• BH：一切都来自boost.random
• sitmo：各种三炸版本
• dqrng：PCG家族，xoshiro256+和xoroshiro128+
• 。

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

#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有一个方便的包装器。

至于教程：除了WRE，Rcpp包小插曲是必读的。Hadley Wickham 的 R Packages 也有一章是关于"编译代码"的，如果你想走devtools路的话。