# Generating random number from a binomial distribution

In my application of Monte Carlo simulation, the crucial part is to generate a random number from a binomial distribution with parameters n = size and p = 0.5. Here is my current implementation

#include <stdio.h>
#include <stdint.h>
#include <x86intrin.h>

int64_t rbinom(int64_t size) {
if (!size) {
return 0;
}

int64_t result = 0;
while (size >= 64) {
uint64_t random64;
while (!_rdrand64_step(&random64)) {
}
result += _popcnt64(random64);
size -= 64;
}

uint64_t random64;
while (!_rdrand64_step(&random64)) {
}
result += _popcnt64(random64 & ~(UINT64_MAX << size));

return result;
}


However, the result of benchmarking terrifies me:

I'm spending 99.68% of the runtime on this function! How can I optimize it?

The result doesn't need to be cryptographically secure, as long as it's good enough for Monte Carlo simulations.

• What is the value of size typically? Your function doesn't seem efficient for large values of size because it requires generating size/64 random numbers. – JS1 Mar 23 '19 at 1:28
• @JS1 It could be very large! Like millions. Generating 64 random bits in a batch doesn’t help a lot in that case. – nalzok Mar 23 '19 at 1:30

The RDRAND instruction generated from the _rdrand64_step intrinsic is actually very slow, though high quality (modulo some concerns about back doors). Depending on the processor it may take hundreds (Ivy Bridge through Skylake) or even thousands (Intel Atom, AMD) of cycles per RDRAND. So just replacing the random number generator will help a lot.

For example, xoroshiro128+ is a relatively fast PRNG, it has some weaknesses but they don't seem too bad for this purpose. An interesting aspect is that it contains no operation that must go to execution port 1 on Intel processors, so its operations do not "fight" the popcnt much, in contrast to PRNGs that contain multiplication.

So overall, something like this:

static inline uint64_t rotl(const uint64_t x, int k) {
return (x << k) | (x >> (64 - k));
}

static uint64_t s[2];

uint64_t next(void) {
const uint64_t s0 = s[0];
uint64_t s1 = s[1];
const uint64_t result = s0 + s1;

s1 ^= s0;
s[0] = rotl(s0, 24) ^ s1 ^ (s1 << 16); // a, b
s[1] = rotl(s1, 37); // c

return result;
}

int64_t rbinom(int64_t size) {
if (!size) {
return 0;
}

int64_t result = 0;
while (size >= 64) {
result += _popcnt64(next());
size -= 64;
}

result += _popcnt64(next() & ~(UINT64_MAX << size));

return result;
}


Elsewhere in the application, the state s must be seeded with a non-zero random-enough number. For example, you could use _rdrand64_step to seed it once, at the start of the application.

But different strategies are possible. With a size in the thousands or even millions (as indicated in the comments), SIMD could be used both to generate pseudo-random bits and to accumulate the pop-counts. Using some techniques from Faster Population Counts Using AVX2 Instructions (mainly, reducing the amount of actual pop-counting by using carry-save addition) and Xorshift+ as PRNG (I avoid rotate because AVX2 does not have them built in, and multiplication because AVX2 also has no 64bit integer multiply built in), it could look like this:

__m256i bigstate0, bigstate1;

__m256i xorshift128plus_avx2(__m256i *state0, __m256i *state1)
{
__m256i s1 = *state0;
const __m256i s0 = *state1;
*state0 = s0;
s1 = _mm256_xor_si256(s1, _mm256_slli_epi64(s1, 23));
*state1 = _mm256_xor_si256(_mm256_xor_si256(_mm256_xor_si256(s1, s0),
_mm256_srli_epi64(s1, 18)),
_mm256_srli_epi64(s0, 5));
}

__m256i popcnt_AVX2(__m256i x) {
const __m256i popcntLUT = _mm256_setr_epi8(
0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
);
const __m256i zero = _mm256_setzero_si256();

__m256i L = _mm256_shuffle_epi8(popcntLUT, _mm256_and_si256(x, nibmask));
x = _mm256_srli_epi16(x, 4);
__m256i H = _mm256_shuffle_epi8(popcntLUT, _mm256_and_si256(x, nibmask));
}

__m256i CSA(__m256i a, __m256i b, __m256i c, __m256i *carry) {
__m256i t0 = _mm256_xor_si256(a, b);
__m256i t1 = _mm256_xor_si256(t0, c);
*carry = _mm256_or_si256(_mm256_and_si256(a, b), _mm256_and_si256(t0, c));
return t1;
}

int64_t rbinom_AVX2(int64_t size) {
if (!size) {
return 0;
}

int64_t result = 0;

__m256i sum1 = _mm256_setzero_si256();
__m256i sum2 = sum1;
__m256i sum4 = sum1;
__m256i sum = sum1;
while (size >= 2048) {
__m256i sample0 = xorshift128plus_avx2(&bigstate0, &bigstate1);
__m256i sample1 = xorshift128plus_avx2(&bigstate0, &bigstate1);
__m256i sample2 = xorshift128plus_avx2(&bigstate0, &bigstate1);
__m256i sample3 = xorshift128plus_avx2(&bigstate0, &bigstate1);
__m256i sample4 = xorshift128plus_avx2(&bigstate0, &bigstate1);
__m256i sample5 = xorshift128plus_avx2(&bigstate0, &bigstate1);
__m256i sample6 = xorshift128plus_avx2(&bigstate0, &bigstate1);
__m256i sample7 = xorshift128plus_avx2(&bigstate0, &bigstate1);
// reduce weight 1
__m256i c0, c1, c2, c3;
__m256i t0 = CSA(sample0, sample1, sample2, &c0);
__m256i t1 = CSA(sample3, sample4, sample5, &c1);
__m256i t2 = CSA(sample6, sample7, sum1, &c2);
sum1 = CSA(t0, t1, t2, &c3);
// reduce weight 2
__m256i c4, c5;
__m256i t3 = CSA(c0, c1, c2, &c4);
sum2 = CSA(c3, t3, sum2, &c5);
// reduce weight 4
__m256i c6;
sum4 = CSA(sum4, c4, c5, &c6);
size -= 2048;
}
sum1 = popcnt_AVX2(sum1);
sum2 = popcnt_AVX2(sum2);
sum4 = popcnt_AVX2(sum4);
result += _mm256_extract_epi64(sum, 0);
result += _mm256_extract_epi64(sum, 1);
result += _mm256_extract_epi64(sum, 2);
result += _mm256_extract_epi64(sum, 3);

while (size >= 64) {
result += _mm_popcnt_u64(next());
size -= 64;
}

result += _mm_popcnt_u64(next() & ~(UINT64_MAX << size));

return result;
}


Algorithmic tricks such as the alias method may be appropriate. I have no experience with this so I cannot explain it or even really recommend it, but it's something to look into.

• Thanks! Is it OK to seed bigstate0 and bigstate1 with arc4random, or do I must seed all static buffers with the same RNG, i.e. _rdrand64_step? I find using _rdrand64_step a little troublesome, due to those error-checking parts. – nalzok Mar 23 '19 at 13:53
• @nalzok arc4random should be OK too – harold Mar 23 '19 at 14:33
• Cool! Another quirky point is I thought I could replace rtol with _rtol64, but surprisingly it doesn't compile on my machine (Core i5), which supports AVX2. I guess I would stay with rtol until I figure out how to write inline assembly. – nalzok Mar 23 '19 at 14:37
• @nalzok it should be fine like this too, many compilers (eg GCC, Clang, MSVC) recognize that that function does a rotate and translate it into one instruction. – harold Mar 23 '19 at 14:40
• @nalzok it would be better, for more reasons than just wider vectors. For example, the CSA function can use _mm512_ternarylogic_epi32 to save some operations. Also future processors will have _mm512_popcnt_epi64 which is obviously great for this. I didn't use it because support for it is still uncommon, but if your processor has it then go for it. – harold Mar 23 '19 at 14:53