A straightforward method for computing the Binomial probability mass function
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.