A straightforward method for computing the Binomial probability mass function

I've been writing up code to test an alternate way of calculating the Binomial probability mass function. The Binomial random variable is the number of heads out of N coin tosses with coins that have probability p of being heads.

The algorithm is pretty simple, we loop through the N coin tosses and update the probability that there are k successes as:

pmf[k] = pmf[k-1] * p +  pmf[k] * (1 - p)


Since we know pmf[k] = 0 for k > n (there can't be more successes than coins tossed so far), we evaluate this line N(N+1)/2 times.

The algorithm works well, is exact, and avoids a number of common pitfalls when computing the Binomial distribution. This is the best implementation I've managed so far, which takes ~190 ms to run on my Dell XPS-15 for N = 15000 and long double p = 0.5.

template <typename T>
std::vector<T> ComputePmf(T p, unsigned int N) {
std::vector<T> pmf(N + 1, 0.0);
pmf[0] = 1.0;

auto k = 0;
for (auto n = 0; n < N; n++) {
for (k = n + 1; k > 0; --k) {
pmf[k] += p * (pmf[k - 1] - pmf[k]);
}
pmf[0] *= (1 - p);
}
return pmf;
}



Ordinarily I'd hang up my hat and move on, but when I implemented the same algorithm in Julia it ran in 150 ms! Since then I've been stuck trying to figure out how to make my c++ version faster - I've tried:

1. using different data types (long double works best, maybe due to SIMD alignment?),
2. std::array versus std::vector (no difference),
3. other ways of factoring the inner loop (this is the best I've come up with),
4. clang and gcc compilers (gcc does slightly better),

It seems like the only real information that I'm not exploiting is the very sequential nature of the data access in this algorithm. Any suggestions?

• It's funny that long double was best, it prevents SIMD from being used (Clang doesn't seem too keen on auto-vectorizing this anyway though) – harold Jul 16 at 0:18
• The difference is remarkable - using double takes ~7 times longer (~1.4 seconds). Maybe this hints at the problem? – combo Jul 16 at 0:21
• Maybe it's running into the denormal penalty? Long double benefits from its significantly larger exponent range then – harold Jul 16 at 0:26
• That would make sense - using float is even slower. Most of these values are very close to zero, some eventually round to zero even with long double values. – combo Jul 16 at 0:32
• @harold you are right about denormals. Setting appropriate compiler flags/macros and switching to doubles made the code run in 50ms. Of course the long double code is still 200ms, but I think julia must be using doubles. If you write up an answer I'll accept it – combo Jul 16 at 3:32

High powers of small-magnitude numbers tend toward zero. When you get subnormal floating-point numbers, then there's a heavy speed penalty; this point is reached sooner with double than with long double (and reached even sooner still using float).
A simple way to improve the speed may be to disable subnormals for this function - e.g. compile with g++ -ffast-math to push subnormal numbers to zero.