在 Rcpp 中有效地生成随机比特流
Generating a random bit stream in Rcpp efficiently
我正在构建的 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.
我前往 the Document, and it says I can use one of the *_rand
function, along with a family of distribution functions。嗯,这很酷,但我的包只需要一个随机位流而不是一个随机 double
。我可以获得它的最简单方法是使用 random(3)
或从 /dev/urandom
读取,但这使我的包 "unportable".
This post 建议使用 sample
,但不幸的是它不适合我的用例。对于我的应用程序,生成随机位显然对性能至关重要,因此我不希望它浪费任何时间调用 unif_rand
,将结果乘以 N
并对其进行舍入。无论如何,我使用 C++ 的原因是为了利用位级并行性。
当然我可以手动滚动我自己的 PRNG 或复制并粘贴最先进的 PRNG 的代码,如 xoshiro256**,但在这样做之前我想看看是否有任何更简单的选择。
顺便说一句,有人可以 link 给我一个不错的 Rcpp 简短教程吗? 编写 R 扩展 非常全面而且很棒,但我需要数周才能完成。我正在寻找更简洁的版本,但最好它应该比调用 Rcpp.package.skeleton
.
提供更多信息
根据@Ralf Stubner的回答,我将原代码重写如下。但是,我每次都得到相同的结果。我怎样才能正确播种并同时保留我的代码 "portable"?
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;
}
有不同的 R 包使 PRNG 可用作 C++ header 仅库:
您可以通过将 LinkingTo
添加到您的包的 DECRIPTION
来使用其中的任何一个。通常,这些 PRNG 是在 C++11 random
header 之后建模的,这意味着您必须控制它们的 life-cycle 并自己播种。在 single-threaded 环境中,我喜欢使用匿名名称空间进行 life-cycle 控制,例如:
#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,Hadley Wickham 的 Rcpp package vignette is a must read. R Packages 也有一章是关于 "compiled code" 的,如果你想走 devtools
-way。
我正在构建的 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.
我前往 the Document, and it says I can use one of the *_rand
function, along with a family of distribution functions。嗯,这很酷,但我的包只需要一个随机位流而不是一个随机 double
。我可以获得它的最简单方法是使用 random(3)
或从 /dev/urandom
读取,但这使我的包 "unportable".
This post 建议使用 sample
,但不幸的是它不适合我的用例。对于我的应用程序,生成随机位显然对性能至关重要,因此我不希望它浪费任何时间调用 unif_rand
,将结果乘以 N
并对其进行舍入。无论如何,我使用 C++ 的原因是为了利用位级并行性。
当然我可以手动滚动我自己的 PRNG 或复制并粘贴最先进的 PRNG 的代码,如 xoshiro256**,但在这样做之前我想看看是否有任何更简单的选择。
顺便说一句,有人可以 link 给我一个不错的 Rcpp 简短教程吗? 编写 R 扩展 非常全面而且很棒,但我需要数周才能完成。我正在寻找更简洁的版本,但最好它应该比调用 Rcpp.package.skeleton
.
根据@Ralf Stubner的回答,我将原代码重写如下。但是,我每次都得到相同的结果。我怎样才能正确播种并同时保留我的代码 "portable"?
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;
}
有不同的 R 包使 PRNG 可用作 C++ header 仅库:
您可以通过将 LinkingTo
添加到您的包的 DECRIPTION
来使用其中的任何一个。通常,这些 PRNG 是在 C++11 random
header 之后建模的,这意味着您必须控制它们的 life-cycle 并自己播种。在 single-threaded 环境中,我喜欢使用匿名名称空间进行 life-cycle 控制,例如:
#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,Hadley Wickham 的 Rcpp package vignette is a must read. R Packages 也有一章是关于 "compiled code" 的,如果你想走 devtools
-way。