A more naive but less complex way to avoid overflows would be to use the original factorial expression but work with logs - in a rough experiment with python for k=10^8 and N=10^17 (numbers plucked out of thin air), the calculation takes ~20s on an M1 MacBook Pro.
This didn't feel too egregious until I looked at the space of inputs and outputs for SHA-256, which is 2^256 out and an absolutely stonking 2^(2^64) - 1 possible inputs.
It feels obvious but I'm still surprised by the efficiency of approximations.
https://en.wikipedia.org/wiki/Stirling's_approximation#Speed...
You can build on that to build good, composable approximations of any the standard combinatorial functions, in log-space (and recover the approximations you want by simply exponentiating back). For example, if you've implemented an approximate ln-fac(n) ~ ln(fac(n)), you immediately get ln-choose, the logarithm of (n!)/(k!(n-k)!), as simply ln-fac(n) - ln-fac(k) - ln-fac(n-k). Fully composable: if ln-fac() is a numerically good approximation, then is so any reasonable sum or difference.
Or: the log of the binomial distribution PDF is simply ln-choose(n,k) + k*ln(p) + (n-k)*ln(1-p).