I work on numerical evaluation of special functions. Fyi, typically a better way to evaluate Gamma ratios is to use the Lanczos approximation. See the Boost math documentation [0] for good explanation of the problems with calculating n! / k! as
exp(lgamma(n+1) - lgamma(k+1)).
The form I use for the approximation is:
gamma(x) = factor(x) * lanczos_sum_expg_scaled(x)
where
factor(x) = ((x + lanczos_g - 0.5)/e)*(x - 0.5)
and lanczos_g is a constant.
When dealing with Gamma ratios, you can often cancel the factors nicely, giving a highly accurate approximation which doesn't overflow. It's curious I saw your comment when I did. I was just tinkering with the Lanczos approximation earlier today.
If you use Python, the Lanczos approximation is available in SciPy using
from scipy.special._ufuncs import _lanczos_sum_expg_scaled
and you can use lanczos_g = 6.024680040776729583740234375
exp(lgamma(n+1) - lgamma(k+1)).
The form I use for the approximation is:
gamma(x) = factor(x) * lanczos_sum_expg_scaled(x)
where
factor(x) = ((x + lanczos_g - 0.5)/e)*(x - 0.5)
and lanczos_g is a constant.
When dealing with Gamma ratios, you can often cancel the factors nicely, giving a highly accurate approximation which doesn't overflow. It's curious I saw your comment when I did. I was just tinkering with the Lanczos approximation earlier today.
If you use Python, the Lanczos approximation is available in SciPy using
from scipy.special._ufuncs import _lanczos_sum_expg_scaled
and you can use lanczos_g = 6.024680040776729583740234375
[0] https://www.boost.org/doc/libs/1_82_0/libs/math/doc/html/mat...