Author: psteitz
Date: Mon Jan 19 15:52:02 2009
New Revision: 735879
URL: http://svn.apache.org/viewvc?rev=735879&view=rev
Log:
Fixed error in binomial coefficient computation
JIRA: MATH-241
Reported and patched by Christian Semrau
Modified:
commons/proper/math/trunk/src/java/org/apache/commons/math/util/MathUtils.java
commons/proper/math/trunk/src/site/xdoc/changes.xml
commons/proper/math/trunk/src/test/org/apache/commons/math/util/MathUtilsTest.java
Modified:
commons/proper/math/trunk/src/java/org/apache/commons/math/util/MathUtils.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/util/MathUtils.java?rev=735879&r1=735878&r2=735879&view=diff
==============================================================================
---
commons/proper/math/trunk/src/java/org/apache/commons/math/util/MathUtils.java
(original)
+++
commons/proper/math/trunk/src/java/org/apache/commons/math/util/MathUtils.java
Mon Jan 19 15:52:02 2009
@@ -181,11 +181,42 @@
if ((k == 1) || (k == n - 1)) {
return n;
}
-
- long result = Math.round(binomialCoefficientDouble(n, k));
- if (result == Long.MAX_VALUE) {
- throw new ArithmeticException(
- "result too large to represent in a long integer");
+ // 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.
+ for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
+ result = result * i / j;
+ }
+ } else if (n <= 66) {
+ // For n > 61 but n <= 66, the result cannot overflow,
+ // but we must take care not to overflow intermediate values.
+ for (int j = 1, i = n - k + 1; j <= k; i++, 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).
+ long d = gcd(i, j);
+ result = (result / (j / d)) * (i / d);
+ }
+ } else {
+ // For n > 66, a result overflow might occur, so we check
+ // the multiplication, taking care to not overflow
+ // unnecessary.
+ for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
+ long d = gcd(i, j);
+ result = mulAndCheck((result / (j / d)), (i / d));
+ }
}
return result;
}
@@ -213,7 +244,33 @@
* @throws IllegalArgumentException if preconditions are not met.
*/
public static double binomialCoefficientDouble(final int n, final int k) {
- return Math.floor(Math.exp(binomialCoefficientLog(n, k)) + 0.5);
+ if (n < k) {
+ throw new IllegalArgumentException(
+ "must have n >= k for binomial coefficient (n,k)");
+ }
+ if (n < 0) {
+ throw new IllegalArgumentException(
+ "must have n >= 0 for binomial coefficient (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 Math.floor(result + 0.5);
}
/**
@@ -247,8 +304,28 @@
return 0;
}
if ((k == 1) || (k == n - 1)) {
- return Math.log((double)n);
+ return Math.log((double) n);
}
+
+ /*
+ * For values small enough to do exact integer computation,
+ * return the log of the exact value
+ */
+ if (n < 67) {
+ return Math.log(binomialCoefficient(n,k));
+ }
+
+ /*
+ * Return the log of binomialCoefficientDouble for values that will not
+ * overflow binomialCoefficientDouble
+ */
+ if (n < 1030) {
+ return Math.log(binomialCoefficientDouble(n, k));
+ }
+
+ /*
+ * Sum logs for values that could overflow
+ */
double logSum = 0;
// n!/k!
@@ -261,7 +338,7 @@
logSum -= Math.log((double)i);
}
- return logSum;
+ return logSum;
}
/**
Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=735879&r1=735878&r2=735879&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Mon Jan 19 15:52:02 2009
@@ -39,6 +39,9 @@
</properties>
<body>
<release version="2.0" date="TBD" description="TBD">
+ <action dev="psteitz" type="fix" issue="MATH-241" due-to="Christian
Semrau">
+ Fixed error in binomial coefficient computation.
+ </action>
<action dev="luc" type="add" >
Added a Legendre-Gauss integrator.
</action>
Modified:
commons/proper/math/trunk/src/test/org/apache/commons/math/util/MathUtilsTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/util/MathUtilsTest.java?rev=735879&r1=735878&r2=735879&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/org/apache/commons/math/util/MathUtilsTest.java
(original)
+++
commons/proper/math/trunk/src/test/org/apache/commons/math/util/MathUtilsTest.java
Mon Jan 19 15:52:02 2009
@@ -14,6 +14,10 @@
package org.apache.commons.math.util;
import java.math.BigDecimal;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
import junit.framework.Test;
import junit.framework.TestCase;
@@ -39,17 +43,37 @@
return suite;
}
+ /** cached binomial coefficients */
+ private static List<Map<Integer, Long>> binomialCache = new
ArrayList<Map<Integer, Long>>();
+
/**
- * Exact recursive implementation to test against
+ * Exact (caching) recursive implementation to test against
*/
- private long binomialCoefficient(int n, int k) {
+ private long binomialCoefficient(int n, int k) throws ArithmeticException {
+ if (binomialCache.size() > n) {
+ Long cachedResult = binomialCache.get(n).get(new Integer(k));
+ if (cachedResult != null) {
+ return cachedResult.longValue();
+ }
+ }
+ long result = -1;
if ((n == k) || (k == 0)) {
- return 1;
+ result = 1;
+ } else if ((k == 1) || (k == n - 1)) {
+ result = n;
+ } else {
+ result = MathUtils.addAndCheck(binomialCoefficient(n - 1, k - 1),
+ binomialCoefficient(n - 1, k));
+ }
+ if (result == -1) {
+ throw new ArithmeticException(
+ "error computing binomial coefficient");
}
- if ((k == 1) || (k == n - 1)) {
- return n;
+ for (int i = binomialCache.size(); i < n + 1; i++) {
+ binomialCache.add(new HashMap<Integer, Long>());
}
- return binomialCoefficient(n - 1, k - 1) + binomialCoefficient(n - 1,
k);
+ binomialCache.get(n).put(new Integer(k), new Long(result));
+ return result;
}
/**
@@ -141,12 +165,63 @@
}
}
- /*
- * Takes a long time for recursion to unwind, but succeeds and yields
- * exact value = 2,333,606,220
- * assertEquals(MathUtils.binomialCoefficient(34,17),
- * binomialCoefficient(34,17));
- */
+ assertEquals(binomialCoefficient(34, 17), MathUtils
+ .binomialCoefficient(34, 17));
+ }
+
+ /**
+ * Tests correctness for large n and sharpness of upper bound in API doc
+ * JIRA: MATH-241
+ */
+ public void testBinomialCoefficientLarge() throws Exception {
+ // This tests all legal and illegal values for n <= 200.
+ for (int n = 0; n <= 200; n++) {
+ for (int k = 0; k <= n; k++) {
+ long ourResult = -1;
+ long exactResult = -1;
+ boolean shouldThrow = false;
+ boolean didThrow = false;
+ try {
+ ourResult = MathUtils.binomialCoefficient(n, k);
+ } catch (ArithmeticException ex) {
+ didThrow = true;
+ }
+ try {
+ exactResult = binomialCoefficient(n, k);
+ } catch (ArithmeticException ex) {
+ shouldThrow = true;
+ }
+ assertEquals(n+","+k, shouldThrow, didThrow);
+ assertEquals(n+","+k, exactResult, ourResult);
+ assertTrue(n+","+k, (n > 66 || !didThrow));
+ }
+ }
+
+ long ourResult = MathUtils.binomialCoefficient(300, 3);
+ long exactResult = binomialCoefficient(300, 3);
+ assertEquals(exactResult, ourResult);
+
+ ourResult = MathUtils.binomialCoefficient(700, 697);
+ exactResult = binomialCoefficient(700, 697);
+ assertEquals(exactResult, ourResult);
+
+ // This one should throw
+ try {
+ MathUtils.binomialCoefficient(700, 300);
+ fail("Expecting ArithmeticException");
+ } catch (ArithmeticException ex) {
+ // Expected
+ }
+
+ // Larger values cannot be computed directly by our
+ // test implementation because of stack limitations,
+ // so we make little jumps to fill the cache.
+ for (int i = 2000; i <= 10000; i += 2000) {
+ ourResult = MathUtils.binomialCoefficient(i, 3);
+ exactResult = binomialCoefficient(i, 3);
+ assertEquals(exactResult, ourResult);
+ }
+
}
public void testBinomialCoefficientFail() {
@@ -171,13 +246,20 @@
;
}
try {
+ MathUtils.binomialCoefficient(67, 30);
+ fail("expecting ArithmeticException");
+ } catch (ArithmeticException ex) {
+ ;
+ }
+ try {
MathUtils.binomialCoefficient(67, 34);
fail("expecting ArithmeticException");
} catch (ArithmeticException ex) {
;
}
double x = MathUtils.binomialCoefficientDouble(1030, 515);
- assertTrue("expecting infinite binomial coefficient",
Double.isInfinite(x));
+ assertTrue("expecting infinite binomial coefficient", Double
+ .isInfinite(x));
}
public void testCosh() {
