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
020import org.apache.commons.numbers.core.ArithmeticUtils;
021
022/**
023 * Representation of the <a href="http://mathworld.wolfram.com/BinomialCoefficient.html">
024 * binomial coefficient</a>.
025 * It is "{@code n choose k}", the number of {@code k}-element subsets that
026 * can be selected from an {@code n}-element set.
027 */
028public final class BinomialCoefficient {
029    /** The maximum m that can be computed without overflow of a long.
030     * {@code C(68, 34) > 2^63}. */
031    private static final int MAX_M = 33;
032    /** The maximum n that can be computed without intermediate overflow for any m.
033     * {@code C(61, 30) * 30 < 2^63}. */
034    private static final int SMALL_N = 61;
035    /** The maximum n that can be computed without overflow of a long for any m.
036     * {@code C(66, 33) < 2^63}. */
037    private static final int LIMIT_N = 66;
038
039    /** Private constructor. */
040    private BinomialCoefficient() {
041        // intentionally empty.
042    }
043
044    /**
045     * Computes the binomial coefficient.
046     *
047     * <p>The largest value of {@code n} for which <em>all</em> coefficients can
048     * fit into a {@code long} is 66. Larger {@code n} may result in an
049     * {@link ArithmeticException} depending on the value of {@code k}.
050     *
051     * <p>Any {@code min(k, n - k) >= 34} cannot fit into a {@code long}
052     * and will result in an {@link ArithmeticException}.
053     *
054     * @param n Size of the set.
055     * @param k Size of the subsets to be counted.
056     * @return {@code n choose k}.
057     * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0} or {@code k > n}.
058     * @throws ArithmeticException if the result is too large to be
059     * represented by a {@code long}.
060     */
061    public static long value(int n, int k) {
062        final int m = checkBinomial(n, k);
063
064        if (m == 0) {
065            return 1;
066        }
067        if (m == 1) {
068            return n;
069        }
070
071        // We use the formulae:
072        // (n choose m) = n! / (n-m)! / m!
073        // (n choose m) = ((n-m+1)*...*n) / (1*...*m)
074        // which can be written
075        // (n choose m) = (n-1 choose m-1) * n / m
076        long result = 1;
077        if (n <= SMALL_N) {
078            // For n <= 61, the naive implementation cannot overflow.
079            int i = n - m + 1;
080            for (int j = 1; j <= m; j++) {
081                result = result * i / j;
082                i++;
083            }
084        } else if (n <= LIMIT_N) {
085            // For n > 61 but n <= 66, the result cannot overflow,
086            // but we must take care not to overflow intermediate values.
087            int i = n - m + 1;
088            for (int j = 1; j <= m; j++) {
089                // We know that (result * i) is divisible by j,
090                // but (result * i) may overflow, so we split j:
091                // Filter out the gcd, d, so j/d and i/d are integer.
092                // result is divisible by (j/d) because (j/d)
093                // is relative prime to (i/d) and is a divisor of
094                // result * (i/d).
095                final long d = ArithmeticUtils.gcd(i, j);
096                result = (result / (j / d)) * (i / d);
097                ++i;
098            }
099        } else {
100            if (m > MAX_M) {
101                throw new ArithmeticException(n + " choose " + k);
102            }
103
104            // For n > 66, a result overflow might occur, so we check
105            // the multiplication, taking care to not overflow
106            // unnecessary.
107            int i = n - m + 1;
108            for (int j = 1; j <= m; j++) {
109                final long d = ArithmeticUtils.gcd(i, j);
110                result = Math.multiplyExact(result / (j / d), i / d);
111                ++i;
112            }
113        }
114
115        return result;
116    }
117
118    /**
119     * Check binomial preconditions.
120     *
121     * <p>For convenience in implementations this returns the smaller of
122     * {@code k} or {@code n - k} allowing symmetry to be exploited in
123     * computing the binomial coefficient.
124     *
125     * @param n Size of the set.
126     * @param k Size of the subsets to be counted.
127     * @return min(k, n - k)
128     * @throws IllegalArgumentException if {@code n < 0}.
129     * @throws IllegalArgumentException if {@code k > n} or {@code k < 0}.
130     */
131    static int checkBinomial(int n,
132                             int k) {
133        // Combine all checks with a single branch:
134        // 0 <= n; 0 <= k <= n
135        // Note: If n >= 0 && k >= 0 && n - k < 0 then k > n.
136        final int m = n - k;
137        // Bitwise or will detect a negative sign bit in any of the numbers
138        if ((n | k | m) < 0) {
139            // Raise the correct exception
140            if (n < 0) {
141                throw new CombinatoricsException(CombinatoricsException.NEGATIVE, n);
142            }
143            throw new CombinatoricsException(CombinatoricsException.OUT_OF_RANGE, k, 0, n);
144        }
145        return m < k ? m : k;
146    }
147}