001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018package org.apache.commons.numbers.combinatorics;
019
020/**
021 * Representation of the <a href="http://mathworld.wolfram.com/BinomialCoefficient.html">
022 * binomial coefficient</a>, as a {@code double}.
023 * It is "{@code n choose k}", the number of {@code k}-element subsets that
024 * can be selected from an {@code n}-element set.
025 */
026public final class BinomialCoefficientDouble {
027    /** The maximum factorial that can be represented as a double. */
028    private static final int MAX_FACTORIAL = 170;
029    /** The maximum n that can be computed without overflow of a long for any m.
030     * {@code C(66, 33) < 2^63}. */
031    private static final int LIMIT_N_LONG = 66;
032    /** The maximum m that can be computed without overflow of a double.
033     * C(1030, 515) ~ 2.85e308. */
034    private static final int MAX_M = 514;
035    /** The maximum n that can be computed without intermediate overflow for any m.
036     * C(1020, 510) * 510 ~ 1.43e308. */
037    private static final int SMALL_N = 1020;
038    /** The maximum m that can be computed without intermediate overflow for any n.
039     * C(2147483647, 37) * 37 ~ 5.13e303. */
040    private static final int SMALL_M = 37;
041
042    /** Private constructor. */
043    private BinomialCoefficientDouble() {
044        // intentionally empty.
045    }
046
047    /**
048     * Computes the binomial coefficient.
049     *
050     * <p>The largest value of {@code n} for which <em>all</em> coefficients can
051     * fit into a {@code double} is 1029. Larger {@code n} may result in
052     * infinity depending on the value of {@code k}.
053     *
054     * <p>Any {@code min(k, n - k) >= 515} cannot fit into a {@code double}
055     * and will result in infinity.
056     *
057     * @param n Size of the set.
058     * @param k Size of the subsets to be counted.
059     * @return {@code n choose k}.
060     * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0} or {@code k > n}.
061     */
062    public static double value(int n, int k) {
063        if (n <= LIMIT_N_LONG) {
064            // Delegate to the exact long result
065            return BinomialCoefficient.value(n, k);
066        }
067        final int m = BinomialCoefficient.checkBinomial(n, k);
068
069        if (m == 0) {
070            return 1;
071        }
072        if (m == 1) {
073            return n;
074        }
075
076        double result;
077        if (n <= MAX_FACTORIAL) {
078            // Small factorials are tabulated exactly
079            // n! / m! / (n-m)!
080            result = Factorial.uncheckedFactorial(n) /
081                     Factorial.uncheckedFactorial(m) /
082                     Factorial.uncheckedFactorial(n - m);
083        } else {
084            // Compute recursively using:
085            // (n choose m) = (n-1 choose m-1) * n / m
086
087            if (n <= SMALL_N || m <= SMALL_M) {
088                // No overflow possible
089                result = 1;
090                for (int i = 1; i <= m; i++) {
091                    result *= n - m + i;
092                    result /= i;
093                }
094            } else {
095                if (m > MAX_M) {
096                    return Double.POSITIVE_INFINITY;
097                }
098
099                // Compute the initial part without overflow checks
100                result = 1;
101                for (int i = 1; i <= SMALL_M; i++) {
102                    result *= n - m + i;
103                    result /= i;
104                }
105                // Careful of overflow
106                for (int i = SMALL_M + 1; i <= m; i++) {
107                    final double next = result * (n - m + i);
108                    if (next > Double.MAX_VALUE) {
109                        // Reverse order of terms
110                        result /= i;
111                        result *= n - m + i;
112                        if (result > Double.MAX_VALUE) {
113                            // Definite overflow
114                            return Double.POSITIVE_INFINITY;
115                        }
116                    } else {
117                        result = next / i;
118                    }
119                }
120            }
121        }
122
123        return Math.floor(result + 0.5);
124    }
125}