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.gamma.LogBeta;
021
022/**
023 * Natural logarithm 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 LogBinomialCoefficient {
029    /** The maximum n that can be computed without overflow of a long for any m.
030     * {@code C(66, 33) < 2^63}. */
031    private static final int LIMIT_N_LONG = 66;
032    /** The maximum n that can be computed without overflow of a double for an m.
033     * C(1029, 514) ~ 1.43e308. */
034    private static final int LIMIT_N_DOUBLE = 1029;
035    /** The maximum m that can be computed without overflow of a double for any n.
036     * C(2147483647, 37) ~ 1.39e302. */
037    private static final int LIMIT_M_DOUBLE = 37;
038
039    /** Private constructor. */
040    private LogBinomialCoefficient() {
041        // intentionally empty.
042    }
043
044    /**
045     * Computes the logarithm of the binomial coefficient.
046     *
047     * <p>This returns a finite result for any valid {@code n choose k}.
048     *
049     * @param n Size of the set.
050     * @param k Size of the subsets to be counted.
051     * @return {@code log(n choose k)}.
052     * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0} or {@code k > n}.
053     */
054    public static double value(int n, int k) {
055        final int m = BinomialCoefficient.checkBinomial(n, k);
056
057        if (m == 0) {
058            return 0;
059        }
060        if (m == 1) {
061            return Math.log(n);
062        }
063
064        if (n <= LIMIT_N_LONG) {
065            // Delegate to the exact long result
066            return Math.log(BinomialCoefficient.value(n, k));
067        }
068        if (n <= LIMIT_N_DOUBLE || m <= LIMIT_M_DOUBLE) {
069            // Delegate to the double result
070            return Math.log(BinomialCoefficientDouble.value(n, k));
071        }
072
073        //    n!                gamma(n+1)               gamma(k+1) * gamma(n-k+1)
074        // ---------   = ------------------------- = 1 / -------------------------
075        // k! (n-k)!     gamma(k+1) * gamma(n-k+1)              gamma(n+1)
076        //
077        //
078        //             = 1 / (k * beta(k, n-k+1))
079        //
080        // where: beta(a, b) = gamma(a) * gamma(b) / gamma(a+b)
081
082        // Delegate to LogBeta
083        return -Math.log(m) - LogBeta.value(m, n - m + 1);
084    }
085}