in commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostBeta.java [1044:1096]
static double binomialCoefficient(int n, int k) {
// Assume: n >= 0; 0 <= k <= n
// Use symmetry
final int m = Math.min(k, n - k);
// This function is typically called with m <= 3 so handle these special cases
if (m == 0) {
return 1;
}
if (m == 1) {
return n;
}
if (m == 2) {
return 0.5 * n * (n - 1);
}
if (m == 3) {
// Divide by 3 at the end to avoid using an 1/6 inexact initial multiplier.
return 0.5 * n * (n - 1) * (n - 2) / 3;
}
double result;
if (n <= MAX_FACTORIAL) {
// Use fast table lookup:
result = BoostGamma.uncheckedFactorial(n);
// Smaller m will have a more accurate factorial
result /= BoostGamma.uncheckedFactorial(m);
result /= BoostGamma.uncheckedFactorial(n - m);
} else {
// Updated:
// Do not use the beta function as per <boost/math/special_functions/binomial.hpp>,
// e.g. result = 1 / (m * beta(m, n - m + 1)) == gamma(n+1) / (m * gamma(m) * gamma(n-m+1))
// Use a summation as per BinomialCoefficientDouble which is more accurate
// than using the beta function. This is only used up to m = 39 so
// may overflow *only* on the final term (i.e. m << n when overflow can occur).
result = 1;
for (int i = 1; i < m; i++) {
result *= n - m + i;
result /= i;
}
// Final term may cause overflow.
if (result * n > Double.MAX_VALUE) {
result /= m;
result *= n;
} else {
result *= n;
result /= m;
}
}
// convert to nearest integer:
return Math.ceil(result - 0.5f);
}