1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 package org.apache.commons.numbers.combinatorics;
19
20 import org.apache.commons.numbers.core.ArithmeticUtils;
21
22 /**
23 * Representation of the <a href="http://mathworld.wolfram.com/BinomialCoefficient.html">
24 * binomial coefficient</a>.
25 * It is "{@code n choose k}", the number of {@code k}-element subsets that
26 * can be selected from an {@code n}-element set.
27 */
28 public final class BinomialCoefficient {
29 /** The maximum m that can be computed without overflow of a long.
30 * C(68, 34) > 2^63. */
31 private static final int MAX_M = 33;
32 /** The maximum n that can be computed without intermediate overflow for any m.
33 * C(61, 30) * 30 < 2^63. */
34 private static final int SMALL_N = 61;
35 /** The maximum n that can be computed without overflow of a long for any m.
36 * C(66, 33) < 2^63. */
37 private static final int LIMIT_N = 66;
38
39 /** Private constructor. */
40 private BinomialCoefficient() {
41 // intentionally empty.
42 }
43
44 /**
45 * Computes the binomial coefficient.
46 *
47 * <p>The largest value of {@code n} for which <em>all</em> coefficients can
48 * fit into a {@code long} is 66. Larger {@code n} may result in an
49 * {@link ArithmeticException} depending on the value of {@code k}.
50 *
51 * <p>Any {@code min(k, n - k) >= 34} cannot fit into a {@code long}
52 * and will result in an {@link ArithmeticException}.
53 *
54 * @param n Size of the set.
55 * @param k Size of the subsets to be counted.
56 * @return {@code n choose k}.
57 * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0} or {@code k > n}.
58 * @throws ArithmeticException if the result is too large to be
59 * represented by a {@code long}.
60 */
61 public static long value(int n, int k) {
62 final int m = checkBinomial(n, k);
63
64 if (m == 0) {
65 return 1;
66 }
67 if (m == 1) {
68 return n;
69 }
70
71 // We use the formulae:
72 // (n choose m) = n! / (n-m)! / m!
73 // (n choose m) = ((n-m+1)*...*n) / (1*...*m)
74 // which can be written
75 // (n choose m) = (n-1 choose m-1) * n / m
76 long result = 1;
77 if (n <= SMALL_N) {
78 // For n <= 61, the naive implementation cannot overflow.
79 int i = n - m + 1;
80 for (int j = 1; j <= m; j++) {
81 result = result * i / j;
82 i++;
83 }
84 } else if (n <= LIMIT_N) {
85 // For n > 61 but n <= 66, the result cannot overflow,
86 // but we must take care not to overflow intermediate values.
87 int i = n - m + 1;
88 for (int j = 1; j <= m; j++) {
89 // We know that (result * i) is divisible by j,
90 // but (result * i) may overflow, so we split j:
91 // Filter out the gcd, d, so j/d and i/d are integer.
92 // result is divisible by (j/d) because (j/d)
93 // is relative prime to (i/d) and is a divisor of
94 // result * (i/d).
95 final long d = ArithmeticUtils.gcd(i, j);
96 result = (result / (j / d)) * (i / d);
97 ++i;
98 }
99 } 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 }