static double binomialCoefficient()

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);
    }