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 /**
21 * Representation of the <a href="http://mathworld.wolfram.com/BinomialCoefficient.html">
22 * binomial coefficient</a>, as a {@code double}.
23 * It is "{@code n choose k}", the number of {@code k}-element subsets that
24 * can be selected from an {@code n}-element set.
25 */
26 public final class BinomialCoefficientDouble {
27 /** The maximum factorial that can be represented as a double. */
28 private static final int MAX_FACTORIAL = 170;
29 /** The maximum n that can be computed without overflow of a long for any m.
30 * C(66, 33) < 2^63. */
31 private static final int LIMIT_N_LONG = 66;
32 /** The maximum m that can be computed without overflow of a double.
33 * C(1030, 515) ~ 2.85e308. */
34 private static final int MAX_M = 514;
35 /** The maximum n that can be computed without intermediate overflow for any m.
36 * C(1020, 510) * 510 ~ 1.43e308. */
37 private static final int SMALL_N = 1020;
38 /** The maximum m that can be computed without intermediate overflow for any n.
39 * C(2147483647, 37) * 37 ~ 5.13e303. */
40 private static final int SMALL_M = 37;
41
42 /** Private constructor. */
43 private BinomialCoefficientDouble() {
44 // intentionally empty.
45 }
46
47 /**
48 * Computes the binomial coefficient.
49 *
50 * <p>The largest value of {@code n} for which <em>all</em> coefficients can
51 * fit into a {@code double} is 1029. Larger {@code n} may result in
52 * infinity depending on the value of {@code k}.
53 *
54 * <p>Any {@code min(k, n - k) >= 515} cannot fit into a {@code double}
55 * and will result in infinity.
56 *
57 * @param n Size of the set.
58 * @param k Size of the subsets to be counted.
59 * @return {@code n choose k}.
60 * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0} or {@code k > n}.
61 */
62 public static double value(int n, int k) {
63 if (n <= LIMIT_N_LONG) {
64 // Delegate to the exact long result
65 return BinomialCoefficient.value(n, k);
66 }
67 final int m = BinomialCoefficient.checkBinomial(n, k);
68
69 if (m == 0) {
70 return 1;
71 }
72 if (m == 1) {
73 return n;
74 }
75
76 double result;
77 if (n <= MAX_FACTORIAL) {
78 // Small factorials are tabulated exactly
79 // n! / m! / (n-m)!
80 result = Factorial.uncheckedFactorial(n) /
81 Factorial.uncheckedFactorial(m) /
82 Factorial.uncheckedFactorial(n - m);
83 } else {
84 // Compute recursively using:
85 // (n choose m) = (n-1 choose m-1) * n / m
86
87 if (n <= SMALL_N || m <= SMALL_M) {
88 // No overflow possible
89 result = 1;
90 for (int i = 1; i <= m; i++) {
91 result *= n - m + i;
92 result /= i;
93 }
94 } else {
95 if (m > MAX_M) {
96 return Double.POSITIVE_INFINITY;
97 }
98
99 // 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 }