View Javadoc
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       * {@code 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 }