A straightforward method for computing the Binomial probability mass function

Question

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?

Answer 1

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.