Here you can find the source of binomialCoefficientLog(final int n, final int k)
Returns the natural log of the <a href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial Coefficient</a>, " n choose k ", the number of k -element subsets that can be selected from an n -element set.
| Parameter | Description |
|---|---|
| n | the size of the set |
| k | the size of the subsets to be counted |
| Parameter | Description |
|---|---|
| NotPositiveException | if n < 0. |
| NumberIsTooLargeException | if k > n. |
| MathArithmeticException | if the result is too large to berepresented by a long integer. |
public static double binomialCoefficientLog(final int n, final int k) throws NotPositiveException, NumberIsTooLargeException, MathArithmeticException
/*/*from w ww . j av a 2 s . c om*/ * Licensed to the Apache Software Foundation (ASF) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. * The ASF licenses this file to You under the Apache License, Version 2.0 * (the "License"); you may not use this file except in compliance with * the License. You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ import java.math.BigInteger; import java.util.concurrent.atomic.AtomicReference; import org.apache.commons.math3.exception.MathArithmeticException; import org.apache.commons.math3.exception.NotPositiveException; import org.apache.commons.math3.exception.NumberIsTooLargeException; import org.apache.commons.math3.exception.util.Localizable; import org.apache.commons.math3.exception.util.LocalizedFormats; public class Main{ /** * Returns the natural {@code log} of the <a * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial * Coefficient</a>, "{@code n choose k}", the number of * {@code k}-element subsets that can be selected from an * {@code n}-element set. * <p> * <Strong>Preconditions</strong>: * <ul> * <li> {@code 0 <= k <= n } (otherwise * {@code IllegalArgumentException} is thrown)</li> * </ul></p> * * @param n the size of the set * @param k the size of the subsets to be counted * @return {@code n choose k} * @throws NotPositiveException if {@code n < 0}. * @throws NumberIsTooLargeException if {@code k > n}. * @throws MathArithmeticException if the result is too large to be * represented by a long integer. */ public static double binomialCoefficientLog(final int n, final int k) throws NotPositiveException, NumberIsTooLargeException, MathArithmeticException { ArithmeticUtils.checkBinomial(n, k); if ((n == k) || (k == 0)) { return 0; } if ((k == 1) || (k == n - 1)) { return FastMath.log(n); } /* * For values small enough to do exact integer computation, * return the log of the exact value */ if (n < 67) { return FastMath.log(binomialCoefficient(n, k)); } /* * Return the log of binomialCoefficientDouble for values that will not * overflow binomialCoefficientDouble */ if (n < 1030) { return FastMath.log(binomialCoefficientDouble(n, k)); } if (k > n / 2) { return binomialCoefficientLog(n, n - k); } /* * Sum logs for values that could overflow */ double logSum = 0; // n!/(n-k)! for (int i = n - k + 1; i <= n; i++) { logSum += FastMath.log(i); } // divide by k! for (int i = 2; i <= k; i++) { logSum -= FastMath.log(i); } return logSum; } /** * Check binomial preconditions. * * @param n Size of the set. * @param k Size of the subsets to be counted. * @throws NotPositiveException if {@code n < 0}. * @throws NumberIsTooLargeException if {@code k > n}. */ private static void checkBinomial(final int n, final int k) throws NumberIsTooLargeException, NotPositiveException { if (n < k) { throw new NumberIsTooLargeException( LocalizedFormats.BINOMIAL_INVALID_PARAMETERS_ORDER, k, n, true); } if (n < 0) { throw new NotPositiveException( LocalizedFormats.BINOMIAL_NEGATIVE_PARAMETER, n); } } /** * Returns an exact representation of the <a * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial * Coefficient</a>, "{@code n choose k}", the number of * {@code k}-element subsets that can be selected from an * {@code n}-element set. * <p> * <Strong>Preconditions</strong>: * <ul> * <li> {@code 0 <= k <= n } (otherwise * {@code IllegalArgumentException} is thrown)</li> * <li> The result is small enough to fit into a {@code long}. The * largest value of {@code n} for which all coefficients are * {@code < Long.MAX_VALUE} is 66. If the computed value exceeds * {@code Long.MAX_VALUE} an {@code ArithMeticException} is * thrown.</li> * </ul></p> * * @param n the size of the set * @param k the size of the subsets to be counted * @return {@code n choose k} * @throws NotPositiveException if {@code n < 0}. * @throws NumberIsTooLargeException if {@code k > n}. * @throws MathArithmeticException if the result is too large to be * represented by a long integer. */ public static long binomialCoefficient(final int n, final int k) throws NotPositiveException, NumberIsTooLargeException, MathArithmeticException { ArithmeticUtils.checkBinomial(n, k); if ((n == k) || (k == 0)) { return 1; } if ((k == 1) || (k == n - 1)) { return n; } // Use symmetry for large k if (k > n / 2) { return binomialCoefficient(n, n - k); } // We use the formula // (n choose k) = n! / (n-k)! / k! // (n choose k) == ((n-k+1)*...*n) / (1*...*k) // which could be written // (n choose k) == (n-1 choose k-1) * n / k long result = 1; if (n <= 61) { // For n <= 61, the naive implementation cannot overflow. int i = n - k + 1; for (int j = 1; j <= k; j++) { result = result * i / j; i++; } } else if (n <= 66) { // For n > 61 but n <= 66, the result cannot overflow, // but we must take care not to overflow intermediate values. int i = n - k + 1; for (int j = 1; j <= k; j++) { // We know that (result * i) is divisible by j, // but (result * i) may overflow, so we split j: // Filter out the gcd, d, so j/d and i/d are integer. // result is divisible by (j/d) because (j/d) // is relative prime to (i/d) and is a divisor of // result * (i/d). final long d = gcd(i, j); result = (result / (j / d)) * (i / d); i++; } } else { // For n > 66, a result overflow might occur, so we check // the multiplication, taking care to not overflow // unnecessary. int i = n - k + 1; for (int j = 1; j <= k; j++) { final long d = gcd(i, j); result = mulAndCheck(result / (j / d), i / d); i++; } } return result; } /** * Returns a {@code double} representation of the <a * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial * Coefficient</a>, "{@code n choose k}", the number of * {@code k}-element subsets that can be selected from an * {@code n}-element set. * <p> * <Strong>Preconditions</strong>: * <ul> * <li> {@code 0 <= k <= n } (otherwise * {@code IllegalArgumentException} is thrown)</li> * <li> The result is small enough to fit into a {@code double}. The * largest value of {@code n} for which all coefficients are < * Double.MAX_VALUE is 1029. If the computed value exceeds Double.MAX_VALUE, * Double.POSITIVE_INFINITY is returned</li> * </ul></p> * * @param n the size of the set * @param k the size of the subsets to be counted * @return {@code n choose k} * @throws NotPositiveException if {@code n < 0}. * @throws NumberIsTooLargeException if {@code k > n}. * @throws MathArithmeticException if the result is too large to be * represented by a long integer. */ public static double binomialCoefficientDouble(final int n, final int k) throws NotPositiveException, NumberIsTooLargeException, MathArithmeticException { ArithmeticUtils.checkBinomial(n, k); if ((n == k) || (k == 0)) { return 1d; } if ((k == 1) || (k == n - 1)) { return n; } if (k > n / 2) { return binomialCoefficientDouble(n, n - k); } if (n < 67) { return binomialCoefficient(n, k); } double result = 1d; for (int i = 1; i <= k; i++) { result *= (double) (n - k + i) / (double) i; } return FastMath.floor(result + 0.5); } /** * Computes the greatest common divisor of the absolute value of two * numbers, using a modified version of the "binary gcd" method. * See Knuth 4.5.2 algorithm B. * The algorithm is due to Josef Stein (1961). * <br/> * Special cases: * <ul> * <li>The invocations * {@code gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)}, * {@code gcd(Integer.MIN_VALUE, 0)} and * {@code gcd(0, Integer.MIN_VALUE)} throw an * {@code ArithmeticException}, because the result would be 2^31, which * is too large for an int value.</li> * <li>The result of {@code gcd(x, x)}, {@code gcd(0, x)} and * {@code gcd(x, 0)} is the absolute value of {@code x}, except * for the special cases above.</li> * <li>The invocation {@code gcd(0, 0)} is the only one which returns * {@code 0}.</li> * </ul> * * @param p Number. * @param q Number. * @return the greatest common divisor (never negative). * @throws MathArithmeticException if the result cannot be represented as * a non-negative {@code int} value. * @since 1.1 */ public static int gcd(int p, int q) throws MathArithmeticException { int a = p; int b = q; if (a == 0 || b == 0) { if (a == Integer.MIN_VALUE || b == Integer.MIN_VALUE) { throw new MathArithmeticException( LocalizedFormats.GCD_OVERFLOW_32_BITS, p, q); } return FastMath.abs(a + b); } long al = a; long bl = b; boolean useLong = false; if (a < 0) { if (Integer.MIN_VALUE == a) { useLong = true; } else { a = -a; } al = -al; } if (b < 0) { if (Integer.MIN_VALUE == b) { useLong = true; } else { b = -b; } bl = -bl; } if (useLong) { if (al == bl) { throw new MathArithmeticException( LocalizedFormats.GCD_OVERFLOW_32_BITS, p, q); } long blbu = bl; bl = al; al = blbu % al; if (al == 0) { if (bl > Integer.MAX_VALUE) { throw new MathArithmeticException( LocalizedFormats.GCD_OVERFLOW_32_BITS, p, q); } return (int) bl; } blbu = bl; // Now "al" and "bl" fit in an "int". b = (int) al; a = (int) (blbu % al); } return gcdPositive(a, b); } /** * <p> * Gets the greatest common divisor of the absolute value of two numbers, * using the "binary gcd" method which avoids division and modulo * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef * Stein (1961). * </p> * Special cases: * <ul> * <li>The invocations * {@code gcd(Long.MIN_VALUE, Long.MIN_VALUE)}, * {@code gcd(Long.MIN_VALUE, 0L)} and * {@code gcd(0L, Long.MIN_VALUE)} throw an * {@code ArithmeticException}, because the result would be 2^63, which * is too large for a long value.</li> * <li>The result of {@code gcd(x, x)}, {@code gcd(0L, x)} and * {@code gcd(x, 0L)} is the absolute value of {@code x}, except * for the special cases above. * <li>The invocation {@code gcd(0L, 0L)} is the only one which returns * {@code 0L}.</li> * </ul> * * @param p Number. * @param q Number. * @return the greatest common divisor, never negative. * @throws MathArithmeticException if the result cannot be represented as * a non-negative {@code long} value. * @since 2.1 */ public static long gcd(final long p, final long q) throws MathArithmeticException { long u = p; long v = q; if ((u == 0) || (v == 0)) { if ((u == Long.MIN_VALUE) || (v == Long.MIN_VALUE)) { throw new MathArithmeticException( LocalizedFormats.GCD_OVERFLOW_64_BITS, p, q); } return FastMath.abs(u) + FastMath.abs(v); } // keep u and v negative, as negative integers range down to // -2^63, while positive numbers can only be as large as 2^63-1 // (i.e. we can't necessarily negate a negative number without // overflow) /* assert u!=0 && v!=0; */ if (u > 0) { u = -u; } // make u negative if (v > 0) { v = -v; } // make v negative // B1. [Find power of 2] int k = 0; while ((u & 1) == 0 && (v & 1) == 0 && k < 63) { // while u and v are // both even... u /= 2; v /= 2; k++; // cast out twos. } if (k == 63) { throw new MathArithmeticException( LocalizedFormats.GCD_OVERFLOW_64_BITS, p, q); } // B2. Initialize: u and v have been divided by 2^k and at least // one is odd. long t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */; // t negative: u was odd, v may be even (t replaces v) // t positive: u was even, v is odd (t replaces u) do { /* assert u<0 && v<0; */ // B4/B3: cast out twos from t. while ((t & 1) == 0) { // while t is even.. t /= 2; // cast out twos } // B5 [reset max(u,v)] if (t > 0) { u = -t; } else { v = t; } // B6/B3. at this point both u and v should be odd. t = (v - u) / 2; // |u| larger: t positive (replace u) // |v| larger: t negative (replace v) } while (t != 0); return -u * (1L << k); // gcd is u*2^k } /** * Multiply two integers, checking for overflow. * * @param x Factor. * @param y Factor. * @return the product {@code x * y}. * @throws MathArithmeticException if the result can not be * represented as an {@code int}. * @since 1.1 */ public static int mulAndCheck(int x, int y) throws MathArithmeticException { long m = ((long) x) * ((long) y); if (m < Integer.MIN_VALUE || m > Integer.MAX_VALUE) { throw new MathArithmeticException(); } return (int) m; } /** * Multiply two long integers, checking for overflow. * * @param a Factor. * @param b Factor. * @return the product {@code a * b}. * @throws MathArithmeticException if the result can not be represented * as a {@code long}. * @since 1.2 */ public static long mulAndCheck(long a, long b) throws MathArithmeticException { long ret; if (a > b) { // use symmetry to reduce boundary cases ret = mulAndCheck(b, a); } else { if (a < 0) { if (b < 0) { // check for positive overflow with negative a, negative b if (a >= Long.MAX_VALUE / b) { ret = a * b; } else { throw new MathArithmeticException(); } } else if (b > 0) { // check for negative overflow with negative a, positive b if (Long.MIN_VALUE / b <= a) { ret = a * b; } else { throw new MathArithmeticException(); } } else { // assert b == 0 ret = 0; } } else if (a > 0) { // assert a > 0 // assert b > 0 // check for positive overflow with positive a, positive b if (a <= Long.MAX_VALUE / b) { ret = a * b; } else { throw new MathArithmeticException(); } } else { // assert a == 0 ret = 0; } } return ret; } /** * Computes the greatest common divisor of two <em>positive</em> numbers * (this precondition is <em>not</em> checked and the result is undefined * if not fulfilled) using the "binary gcd" method which avoids division * and modulo operations. * See Knuth 4.5.2 algorithm B. * The algorithm is due to Josef Stein (1961). * <br/> * Special cases: * <ul> * <li>The result of {@code gcd(x, x)}, {@code gcd(0, x)} and * {@code gcd(x, 0)} is the value of {@code x}.</li> * <li>The invocation {@code gcd(0, 0)} is the only one which returns * {@code 0}.</li> * </ul> * * @param a Positive number. * @param b Positive number. * @return the greatest common divisor. */ private static int gcdPositive(int a, int b) { if (a == 0) { return b; } else if (b == 0) { return a; } // Make "a" and "b" odd, keeping track of common power of 2. final int aTwos = Integer.numberOfTrailingZeros(a); a >>= aTwos; final int bTwos = Integer.numberOfTrailingZeros(b); b >>= bTwos; final int shift = Math.min(aTwos, bTwos); // "a" and "b" are positive. // If a > b then "gdc(a, b)" is equal to "gcd(a - b, b)". // If a < b then "gcd(a, b)" is equal to "gcd(b - a, a)". // Hence, in the successive iterations: // "a" becomes the absolute difference of the current values, // "b" becomes the minimum of the current values. while (a != b) { final int delta = a - b; b = Math.min(a, b); a = Math.abs(delta); // Remove any power of 2 in "a" ("b" is guaranteed to be odd). a >>= Integer.numberOfTrailingZeros(a); } // Recover the common power of 2. return a << shift; } }