I am currently struggling with calculating binomial coefficients for very large numbers, let's say "n choose k" with n < 10,000,000 and n < k. This is necessary in the context of computing hypergeometric probability distributions.

Up to this point, I have tried many approaches to handle the large numbers that result from these calculations. However, the problem is that I don't need to calculate these binomial coefficients once but hundred of thousands of times. That means that normal approaches of computing factorials is by far too expensive and standard data types like `long long int` are too limited and cannot hold these numbers. 

I've already tried multiprecision data types from the `Boost` library but as I mentioned, doing calculations so many times results in extremely slow performance. I've also tried multi-threading using `OpenMP` but the performance benefit was still way too low. As a result, I've switched to calculating the logarithm of binomial coefficients to keep the numbers small. Although this solved the problem of large numbers, this did not speed up the process. That's why I've tried out the Stirling approximation of logarithmic binomial coefficients. My current solution looks like this:

```
#include <math.h>

long double calc_hgeom(unsigned int k, unsigned int n, unsigned int K, unsigned int N)
{
    long double hprob = std::exp((log_C(K, k) + log_C(N-K, n-k)) - log_C(N, n));
    return hprob;
}

long double log_C(unsigned int u, unsigned int m)
{
    long double C = u * std::log(u) - m * std::log(m) - (u-m) * std::log(u-m)) + 0.5 * (std::log(u) - std::log(m) - std::log(u-m) - std::log(2*M_PI));
    return C;
}

```

However, the results differ quite a lot from the actual values, up to 7 %. Hence my question: Is there an efficient way to calculate the logarithm of binomial coefficients or could my approximation be improved to increase accuracy?

Any help would be much appreciated since this calculation is the basis of my whole algorithm.