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.gamma.LogBeta;
21
22 /**
23 * Natural logarithm 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 LogBinomialCoefficient {
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 n that can be computed without overflow of a double for an m.
33 * C(1029, 514) ~ 1.43e308. */
34 private static final int LIMIT_N_DOUBLE = 1029;
35 /** The maximum m that can be computed without overflow of a double for any n.
36 * C(2147483647, 37) ~ 1.39e302. */
37 private static final int LIMIT_M_DOUBLE = 37;
38
39 /** Private constructor. */
40 private LogBinomialCoefficient() {
41 // intentionally empty.
42 }
43
44 /**
45 * Computes the logarithm of the binomial coefficient.
46 *
47 * <p>This returns a finite result for any valid {@code n choose k}.
48 *
49 * @param n Size of the set.
50 * @param k Size of the subsets to be counted.
51 * @return {@code log(n choose k)}.
52 * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0} or {@code k > n}.
53 */
54 public static double value(int n, int k) {
55 final int m = BinomialCoefficient.checkBinomial(n, k);
56
57 if (m == 0) {
58 return 0;
59 }
60 if (m == 1) {
61 return Math.log(n);
62 }
63
64 if (n <= LIMIT_N_LONG) {
65 // Delegate to the exact long result
66 return Math.log(BinomialCoefficient.value(n, k));
67 }
68 if (n <= LIMIT_N_DOUBLE || m <= LIMIT_M_DOUBLE) {
69 // Delegate to the double result
70 return Math.log(BinomialCoefficientDouble.value(n, k));
71 }
72
73 // n! gamma(n+1) gamma(k+1) * gamma(n-k+1)
74 // --------- = ------------------------- = 1 / -------------------------
75 // k! (n-k)! gamma(k+1) * gamma(n-k+1) gamma(n+1)
76 //
77 //
78 // = 1 / (k * beta(k, n-k+1))
79 //
80 // where: beta(a, b) = gamma(a) * gamma(b) / gamma(a+b)
81
82 // Delegate to LogBeta
83 return -Math.log(m) - LogBeta.value(m, n - m + 1);
84 }
85 }