jeudi 28 mars 2019

Generating a random bit stream in Rcpp efficiently

I have an auxiliary function in the R package I'm currently building named rbinom01. Note that it calls 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;
}

When R CMD check my_package, I got the following warning:

* 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.

I headed to the Document, and it says I can use one of the *_rand function, along with a family of distribution functions. Well that's cool, but my package simply needs a stream of random bits rather than a random double. The easiest way I can have it is by using random(3) or maybe reading from /dev/urandom, but that makes my package "unportable".

This post suggests using sample, but unfortunately it doesn't fit into my use case. For my application, generating random bits is apparently critical to the performance, so I don't want it waste any time calling unif_rand, multiply the result by N and round it. Anyway, the reason I'm using C++ is to exploit bit-level parallelism.

Surely I can hand-roll my own PRNG or copy and paste the code of a state-of-the-art PRNG like xoshiro256**, but before doing that I would like to see if there are any easier alternatives.

Incidentally, could someone please link a nice short tutorial of Rcpp to me? Writing R Extensions is comprehensive and awesome but it would take my weeks to finish. I'm looking for a more concise version, but preferably it should be more informative than a call to Rcpp.package.skeleton.




Aucun commentaire:

Enregistrer un commentaire