The R Project SVN R

#include <math.h>

#include <stdint.h>

//copied from src/main/RNG.c:

//generate a random non-negative integer < 2 ^ bits in 16 bit chunks

static double rbits(int bits)

{

    int_least64_t v = 0;

    for (int n = 0; n <= bits; n += 16) {

	int v1 = (int) floor(unif_rand() * 65536);

	v = 65536 * v + v1;

    }

    // mask out the bits in the result that are not needed

    return (double) (v & ((1L << bits) - 1));

}

double R_unif_index(double dn)

{

    // rejection sampling from integers below the next larger power of two

    if (dn <= 0)

	return 0.0;

    int bits = (int) ceil(log2(dn));

    double dv;

    do { dv = rbits(bits); } while (dn <= dv);

    return dv;

}