From 6a5fe463eacb5087fade3e6a1e58371a8789bafe Mon Sep 17 00:00:00 2001
From: Sebastien Brisard
Date: Thu, 13 Oct 2011 05:29:28 +0000
Subject: [PATCH] Finished moving methods from MathUtils to ArithmeticsUtils
(MATH-689)
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1182658 13f79535-47bb-0310-9956-ffa450edef68
---
.../polynomials/PolynomialsUtils.java | 4 +-
.../HypergeometricDistributionImpl.java | 8 +-
.../distribution/PascalDistributionImpl.java | 4 +-
.../commons/math/fraction/Fraction.java | 13 +-
.../commons/math/util/ArithmeticsUtils.java | 487 ++++++++++++++---
.../apache/commons/math/util/MathUtils.java | 349 ------------
.../polynomials/PolynomialsUtilsTest.java | 4 +-
.../math/linear/InverseHilbertMatrix.java | 16 +-
.../math/util/ArithmeticsUtilsTest.java | 514 +++++++++++++++---
.../commons/math/util/MathUtilsTest.java | 364 -------------
10 files changed, 880 insertions(+), 883 deletions(-)
diff --git a/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.java b/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.java
index 65ad670d3..b26ad45dd 100644
--- a/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.java
+++ b/src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.java
@@ -22,8 +22,8 @@ import java.util.List;
import java.util.Map;
import org.apache.commons.math.fraction.BigFraction;
+import org.apache.commons.math.util.ArithmeticsUtils;
import org.apache.commons.math.util.FastMath;
-import org.apache.commons.math.util.MathUtils;
/**
* A collection of static methods that operate on or return polynomials.
@@ -326,7 +326,7 @@ public class PolynomialsUtils {
final int[][] coeff = new int[dp1][dp1];
for (int i = 0; i < dp1; i++){
for(int j = 0; j <= i; j++){
- coeff[i][j] = (int) MathUtils.binomialCoefficient(i, j);
+ coeff[i][j] = (int) ArithmeticsUtils.binomialCoefficient(i, j);
}
}
diff --git a/src/main/java/org/apache/commons/math/distribution/HypergeometricDistributionImpl.java b/src/main/java/org/apache/commons/math/distribution/HypergeometricDistributionImpl.java
index 5c726fa36..d32db04f5 100644
--- a/src/main/java/org/apache/commons/math/distribution/HypergeometricDistributionImpl.java
+++ b/src/main/java/org/apache/commons/math/distribution/HypergeometricDistributionImpl.java
@@ -23,7 +23,7 @@ import org.apache.commons.math.exception.NotPositiveException;
import org.apache.commons.math.exception.NotStrictlyPositiveException;
import org.apache.commons.math.exception.NumberIsTooLargeException;
import org.apache.commons.math.exception.util.LocalizedFormats;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.ArithmeticsUtils;
import org.apache.commons.math.util.FastMath;
/**
@@ -231,9 +231,9 @@ public class HypergeometricDistributionImpl extends AbstractIntegerDistribution
* @return PMF for the distribution.
*/
private double probability(int n, int m, int k, int x) {
- return FastMath.exp(MathUtils.binomialCoefficientLog(m, x) +
- MathUtils.binomialCoefficientLog(n - m, k - x) -
- MathUtils.binomialCoefficientLog(n, k));
+ return FastMath.exp(ArithmeticsUtils.binomialCoefficientLog(m, x) +
+ ArithmeticsUtils.binomialCoefficientLog(n - m, k - x) -
+ ArithmeticsUtils.binomialCoefficientLog(n, k));
}
/**
diff --git a/src/main/java/org/apache/commons/math/distribution/PascalDistributionImpl.java b/src/main/java/org/apache/commons/math/distribution/PascalDistributionImpl.java
index d9f614411..a81e0524c 100644
--- a/src/main/java/org/apache/commons/math/distribution/PascalDistributionImpl.java
+++ b/src/main/java/org/apache/commons/math/distribution/PascalDistributionImpl.java
@@ -22,7 +22,7 @@ import org.apache.commons.math.exception.OutOfRangeException;
import org.apache.commons.math.exception.NotPositiveException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.special.Beta;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.ArithmeticsUtils;
import org.apache.commons.math.util.FastMath;
/**
@@ -128,7 +128,7 @@ public class PascalDistributionImpl extends AbstractIntegerDistribution
if (x < 0) {
ret = 0.0;
} else {
- ret = MathUtils.binomialCoefficientDouble(x +
+ ret = ArithmeticsUtils.binomialCoefficientDouble(x +
numberOfSuccesses - 1, numberOfSuccesses - 1) *
FastMath.pow(probabilityOfSuccess, numberOfSuccesses) *
FastMath.pow(1.0 - probabilityOfSuccess, x);
diff --git a/src/main/java/org/apache/commons/math/fraction/Fraction.java b/src/main/java/org/apache/commons/math/fraction/Fraction.java
index 393cc4792..9711a5ca3 100644
--- a/src/main/java/org/apache/commons/math/fraction/Fraction.java
+++ b/src/main/java/org/apache/commons/math/fraction/Fraction.java
@@ -24,7 +24,6 @@ import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.exception.MathArithmeticException;
import org.apache.commons.math.exception.NullArgumentException;
import org.apache.commons.math.util.ArithmeticsUtils;
-import org.apache.commons.math.util.MathUtils;
import org.apache.commons.math.util.FastMath;
/**
@@ -490,12 +489,12 @@ public class Fraction
int d1 = ArithmeticsUtils.gcd(denominator, fraction.denominator);
if (d1==1) {
// result is ( (u*v' +/- u'v) / u'v')
- int uvp = MathUtils.mulAndCheck(numerator, fraction.denominator);
- int upv = MathUtils.mulAndCheck(fraction.numerator, denominator);
+ int uvp = ArithmeticsUtils.mulAndCheck(numerator, fraction.denominator);
+ int upv = ArithmeticsUtils.mulAndCheck(fraction.numerator, denominator);
return new Fraction
(isAdd ? ArithmeticsUtils.addAndCheck(uvp, upv) :
ArithmeticsUtils.subAndCheck(uvp, upv),
- MathUtils.mulAndCheck(denominator, fraction.denominator));
+ ArithmeticsUtils.mulAndCheck(denominator, fraction.denominator));
}
// the quantity 't' requires 65 bits of precision; see knuth 4.5.1
// exercise 7. we're going to use a BigInteger.
@@ -517,7 +516,7 @@ public class Fraction
w);
}
return new Fraction (w.intValue(),
- MathUtils.mulAndCheck(denominator/d1,
+ ArithmeticsUtils.mulAndCheck(denominator/d1,
fraction.denominator/d2));
}
@@ -543,8 +542,8 @@ public class Fraction
int d1 = ArithmeticsUtils.gcd(numerator, fraction.denominator);
int d2 = ArithmeticsUtils.gcd(fraction.numerator, denominator);
return getReducedFraction
- (MathUtils.mulAndCheck(numerator/d1, fraction.numerator/d2),
- MathUtils.mulAndCheck(denominator/d2, fraction.denominator/d1));
+ (ArithmeticsUtils.mulAndCheck(numerator/d1, fraction.numerator/d2),
+ ArithmeticsUtils.mulAndCheck(denominator/d2, fraction.denominator/d1));
}
/**
diff --git a/src/main/java/org/apache/commons/math/util/ArithmeticsUtils.java b/src/main/java/org/apache/commons/math/util/ArithmeticsUtils.java
index d757d9cf7..e2248a03c 100644
--- a/src/main/java/org/apache/commons/math/util/ArithmeticsUtils.java
+++ b/src/main/java/org/apache/commons/math/util/ArithmeticsUtils.java
@@ -18,6 +18,7 @@ package org.apache.commons.math.util;
import org.apache.commons.math.exception.MathArithmeticException;
import org.apache.commons.math.exception.NotPositiveException;
+import org.apache.commons.math.exception.NumberIsTooLargeException;
import org.apache.commons.math.exception.util.Localizable;
import org.apache.commons.math.exception.util.LocalizedFormats;
@@ -77,92 +78,193 @@ public final class ArithmeticsUtils {
}
/**
- * Add two long integers, checking for overflow.
+ * Returns an exact representation of the Binomial
+ * Coefficient, "{@code n choose k}", the number of
+ * {@code k}-element subsets that can be selected from an
+ * {@code n}-element set.
+ *
+ * Preconditions:
+ *
+ * - {@code 0 <= k <= n } (otherwise
+ * {@code IllegalArgumentException} is thrown)
+ * - 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.
+ *
*
- * @param a Addend.
- * @param b Addend.
- * @param pattern Pattern to use for any thrown exception.
- * @return the sum {@code a + b}.
- * @throws MathArithmeticException if the result cannot be represented
- * as a {@code long}.
- * @since 1.2
+ * @param n the size of the set
+ * @param k the size of the subsets to be counted
+ * @return {@code n choose k}
+ * @throws MathIllegalArgumentException if preconditions are not met.
+ * @throws MathArithmeticException if the result is too large to be
+ * represented by a long integer.
*/
- private static long addAndCheck(long a, long b, Localizable pattern) {
- long ret;
- if (a > b) {
- // use symmetry to reduce boundary cases
- ret = addAndCheck(b, a, pattern);
+ public static long binomialCoefficient(final int n, final int k) {
+ ArithmeticsUtils.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 {
- // assert a <= b
-
- if (a < 0) {
- if (b < 0) {
- // check for negative overflow
- if (Long.MIN_VALUE - b <= a) {
- ret = a + b;
- } else {
- throw new MathArithmeticException(pattern, a, b);
- }
- } else {
- // opposite sign addition is always safe
- ret = a + b;
- }
- } else {
- // assert a >= 0
- // assert b >= 0
-
- // check for positive overflow
- if (a <= Long.MAX_VALUE - b) {
- ret = a + b;
- } else {
- throw new MathArithmeticException(pattern, a, b);
- }
+ // 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 ret;
+ return result;
}
/**
- * Subtract two integers, checking for overflow.
+ * Returns a {@code double} representation of the Binomial
+ * Coefficient, "{@code n choose k}", the number of
+ * {@code k}-element subsets that can be selected from an
+ * {@code n}-element set.
+ *
+ * Preconditions:
+ *
+ * - {@code 0 <= k <= n } (otherwise
+ * {@code IllegalArgumentException} is thrown)
+ * - 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
+ *
*
- * @param x Minuend.
- * @param y Subtrahend.
- * @return the difference {@code x - y}.
- * @throws MathArithmeticException if the result can not be represented
- * as an {@code int}.
- * @since 1.1
+ * @param n the size of the set
+ * @param k the size of the subsets to be counted
+ * @return {@code n choose k}
+ * @throws IllegalArgumentException if preconditions are not met.
*/
- public static int subAndCheck(int x, int y) {
- long s = (long)x - (long)y;
- if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
- throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_SUBTRACTION, x, y);
+ public static double binomialCoefficientDouble(final int n, final int k) {
+ ArithmeticsUtils.checkBinomial(n, k);
+ if ((n == k) || (k == 0)) {
+ return 1d;
}
- return (int)s;
+ 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);
}
/**
- * Subtract two long integers, checking for overflow.
+ * Returns the natural {@code log} of the Binomial
+ * Coefficient, "{@code n choose k}", the number of
+ * {@code k}-element subsets that can be selected from an
+ * {@code n}-element set.
+ *
+ * Preconditions:
+ *
+ * - {@code 0 <= k <= n } (otherwise
+ * {@code IllegalArgumentException} is thrown)
+ *
*
- * @param a Value.
- * @param b Value.
- * @return the difference {@code a - b}.
- * @throws MathArithmeticException if the result can not be represented as a
- * {@code long}.
- * @since 1.2
+ * @param n the size of the set
+ * @param k the size of the subsets to be counted
+ * @return {@code n choose k}
+ * @throws IllegalArgumentException if preconditions are not met.
*/
- public static long subAndCheck(long a, long b) {
- long ret;
- if (b == Long.MIN_VALUE) {
- if (a < 0) {
- ret = a - b;
- } else {
- throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, a, -b);
- }
- } else {
- // use additive inverse
- ret = addAndCheck(a, -b, LocalizedFormats.OVERFLOW_IN_ADDITION);
+ public static double binomialCoefficientLog(final int n, final int k) {
+ ArithmeticsUtils.checkBinomial(n, k);
+ if ((n == k) || (k == 0)) {
+ return 0;
}
- return ret;
+ 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;
}
/**
@@ -419,4 +521,251 @@ public final class ArithmeticsUtils {
} while (t != 0);
return -u * (1L << k); // gcd is u*2^k
}
+
+ /**
+ *
+ * Returns the least common multiple of the absolute value of two numbers,
+ * using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}.
+ *
+ * Special cases:
+ *
+ * - The invocations {@code lcm(Integer.MIN_VALUE, n)} and
+ * {@code lcm(n, Integer.MIN_VALUE)}, where {@code abs(n)} is a
+ * power of 2, throw an {@code ArithmeticException}, because the result
+ * would be 2^31, which is too large for an int value.
+ * - The result of {@code lcm(0, x)} and {@code lcm(x, 0)} is
+ * {@code 0} for any {@code x}.
+ *
+ *
+ * @param a Number.
+ * @param b Number.
+ * @return the least common multiple, never negative.
+ * @throws MathArithmeticException if the result cannot be represented as
+ * a non-negative {@code int} value.
+ * @since 1.1
+ */
+ public static int lcm(int a, int b) {
+ if (a == 0 || b == 0){
+ return 0;
+ }
+ int lcm = FastMath.abs(ArithmeticsUtils.mulAndCheck(a / gcd(a, b), b));
+ if (lcm == Integer.MIN_VALUE) {
+ throw new MathArithmeticException(LocalizedFormats.LCM_OVERFLOW_32_BITS,
+ a, b);
+ }
+ return lcm;
+ }
+
+ /**
+ *
+ * Returns the least common multiple of the absolute value of two numbers,
+ * using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}.
+ *
+ * Special cases:
+ *
+ * - The invocations {@code lcm(Long.MIN_VALUE, n)} and
+ * {@code lcm(n, Long.MIN_VALUE)}, where {@code abs(n)} is a
+ * power of 2, throw an {@code ArithmeticException}, because the result
+ * would be 2^63, which is too large for an int value.
+ * - The result of {@code lcm(0L, x)} and {@code lcm(x, 0L)} is
+ * {@code 0L} for any {@code x}.
+ *
+ *
+ * @param a Number.
+ * @param b Number.
+ * @return the least common multiple, never negative.
+ * @throws MathArithmeticException if the result cannot be represented
+ * as a non-negative {@code long} value.
+ * @since 2.1
+ */
+ public static long lcm(long a, long b) {
+ if (a == 0 || b == 0){
+ return 0;
+ }
+ long lcm = FastMath.abs(ArithmeticsUtils.mulAndCheck(a / gcd(a, b), b));
+ if (lcm == Long.MIN_VALUE){
+ throw new MathArithmeticException(LocalizedFormats.LCM_OVERFLOW_64_BITS,
+ a, b);
+ }
+ return lcm;
+ }
+
+ /**
+ * 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) {
+ 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) {
+ 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;
+ }
+
+ /**
+ * Subtract two integers, checking for overflow.
+ *
+ * @param x Minuend.
+ * @param y Subtrahend.
+ * @return the difference {@code x - y}.
+ * @throws MathArithmeticException if the result can not be represented
+ * as an {@code int}.
+ * @since 1.1
+ */
+ public static int subAndCheck(int x, int y) {
+ long s = (long)x - (long)y;
+ if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
+ throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_SUBTRACTION, x, y);
+ }
+ return (int)s;
+ }
+
+ /**
+ * Subtract two long integers, checking for overflow.
+ *
+ * @param a Value.
+ * @param b Value.
+ * @return the difference {@code a - b}.
+ * @throws MathArithmeticException if the result can not be represented as a
+ * {@code long}.
+ * @since 1.2
+ */
+ public static long subAndCheck(long a, long b) {
+ long ret;
+ if (b == Long.MIN_VALUE) {
+ if (a < 0) {
+ ret = a - b;
+ } else {
+ throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, a, -b);
+ }
+ } else {
+ // use additive inverse
+ ret = addAndCheck(a, -b, LocalizedFormats.OVERFLOW_IN_ADDITION);
+ }
+ return ret;
+ }
+
+ /**
+ * Add two long integers, checking for overflow.
+ *
+ * @param a Addend.
+ * @param b Addend.
+ * @param pattern Pattern to use for any thrown exception.
+ * @return the sum {@code a + b}.
+ * @throws MathArithmeticException if the result cannot be represented
+ * as a {@code long}.
+ * @since 1.2
+ */
+ private static long addAndCheck(long a, long b, Localizable pattern) {
+ long ret;
+ if (a > b) {
+ // use symmetry to reduce boundary cases
+ ret = addAndCheck(b, a, pattern);
+ } else {
+ // assert a <= b
+
+ if (a < 0) {
+ if (b < 0) {
+ // check for negative overflow
+ if (Long.MIN_VALUE - b <= a) {
+ ret = a + b;
+ } else {
+ throw new MathArithmeticException(pattern, a, b);
+ }
+ } else {
+ // opposite sign addition is always safe
+ ret = a + b;
+ }
+ } else {
+ // assert a >= 0
+ // assert b >= 0
+
+ // check for positive overflow
+ if (a <= Long.MAX_VALUE - b) {
+ ret = a + b;
+ } else {
+ throw new MathArithmeticException(pattern, a, b);
+ }
+ }
+ }
+ return ret;
+ }
+
+ /**
+ * 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) {
+ 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);
+ }
+ }
}
diff --git a/src/main/java/org/apache/commons/math/util/MathUtils.java b/src/main/java/org/apache/commons/math/util/MathUtils.java
index e33b0bf30..0145c25fd 100644
--- a/src/main/java/org/apache/commons/math/util/MathUtils.java
+++ b/src/main/java/org/apache/commons/math/util/MathUtils.java
@@ -26,7 +26,6 @@ import org.apache.commons.math.exception.MathIllegalArgumentException;
import org.apache.commons.math.exception.NotFiniteNumberException;
import org.apache.commons.math.exception.NotPositiveException;
import org.apache.commons.math.exception.NullArgumentException;
-import org.apache.commons.math.exception.NumberIsTooLargeException;
import org.apache.commons.math.exception.util.Localizable;
import org.apache.commons.math.exception.util.LocalizedFormats;
@@ -76,214 +75,6 @@ public final class MathUtils {
super();
}
- /**
- * Returns an exact representation of the Binomial
- * Coefficient, "{@code n choose k}", the number of
- * {@code k}-element subsets that can be selected from an
- * {@code n}-element set.
- *
- * Preconditions:
- *
- * - {@code 0 <= k <= n } (otherwise
- * {@code IllegalArgumentException} is thrown)
- * - 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.
- *
- *
- * @param n the size of the set
- * @param k the size of the subsets to be counted
- * @return {@code n choose k}
- * @throws MathIllegalArgumentException if preconditions are not met.
- * @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) {
- 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 = ArithmeticsUtils.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 = ArithmeticsUtils.gcd(i, j);
- result = mulAndCheck(result / (j / d), i / d);
- i++;
- }
- }
- return result;
- }
-
- /**
- * Returns a {@code double} representation of the Binomial
- * Coefficient, "{@code n choose k}", the number of
- * {@code k}-element subsets that can be selected from an
- * {@code n}-element set.
- *
- * Preconditions:
- *
- * - {@code 0 <= k <= n } (otherwise
- * {@code IllegalArgumentException} is thrown)
- * - 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
- *
- *
- * @param n the size of the set
- * @param k the size of the subsets to be counted
- * @return {@code n choose k}
- * @throws IllegalArgumentException if preconditions are not met.
- */
- public static double binomialCoefficientDouble(final int n, final int k) {
- 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);
- }
-
- /**
- * Returns the natural {@code log} of the Binomial
- * Coefficient, "{@code n choose k}", the number of
- * {@code k}-element subsets that can be selected from an
- * {@code n}-element set.
- *
- * Preconditions:
- *
- * - {@code 0 <= k <= n } (otherwise
- * {@code IllegalArgumentException} is thrown)
- *
- *
- * @param n the size of the set
- * @param k the size of the subsets to be counted
- * @return {@code n choose k}
- * @throws IllegalArgumentException if preconditions are not met.
- */
- public static double binomialCoefficientLog(final int n, final int k) {
- 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) {
- 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 the
* hyperbolic cosine of x.
@@ -387,74 +178,6 @@ public final class MathUtils {
return (x >= ZS) ? PS : NS;
}
- /**
- *
- * Returns the least common multiple of the absolute value of two numbers,
- * using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}.
- *
- * Special cases:
- *
- * - The invocations {@code lcm(Integer.MIN_VALUE, n)} and
- * {@code lcm(n, Integer.MIN_VALUE)}, where {@code abs(n)} is a
- * power of 2, throw an {@code ArithmeticException}, because the result
- * would be 2^31, which is too large for an int value.
- * - The result of {@code lcm(0, x)} and {@code lcm(x, 0)} is
- * {@code 0} for any {@code x}.
- *
- *
- * @param a Number.
- * @param b Number.
- * @return the least common multiple, never negative.
- * @throws MathArithmeticException if the result cannot be represented as
- * a non-negative {@code int} value.
- * @since 1.1
- */
- public static int lcm(int a, int b) {
- if (a == 0 || b == 0){
- return 0;
- }
- int lcm = FastMath.abs(mulAndCheck(a / ArithmeticsUtils.gcd(a, b), b));
- if (lcm == Integer.MIN_VALUE) {
- throw new MathArithmeticException(LocalizedFormats.LCM_OVERFLOW_32_BITS,
- a, b);
- }
- return lcm;
- }
-
- /**
- *
- * Returns the least common multiple of the absolute value of two numbers,
- * using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}.
- *
- * Special cases:
- *
- * - The invocations {@code lcm(Long.MIN_VALUE, n)} and
- * {@code lcm(n, Long.MIN_VALUE)}, where {@code abs(n)} is a
- * power of 2, throw an {@code ArithmeticException}, because the result
- * would be 2^63, which is too large for an int value.
- * - The result of {@code lcm(0L, x)} and {@code lcm(x, 0L)} is
- * {@code 0L} for any {@code x}.
- *
- *
- * @param a Number.
- * @param b Number.
- * @return the least common multiple, never negative.
- * @throws MathArithmeticException if the result cannot be represented
- * as a non-negative {@code long} value.
- * @since 2.1
- */
- public static long lcm(long a, long b) {
- if (a == 0 || b == 0){
- return 0;
- }
- long lcm = FastMath.abs(mulAndCheck(a / ArithmeticsUtils.gcd(a, b), b));
- if (lcm == Long.MIN_VALUE){
- throw new MathArithmeticException(LocalizedFormats.LCM_OVERFLOW_64_BITS,
- a, b);
- }
- return lcm;
- }
-
/**
* Returns the
* logarithm
@@ -475,78 +198,6 @@ public final class MathUtils {
return FastMath.log(x)/FastMath.log(base);
}
- /**
- * 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) {
- 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) {
- 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;
- }
-
/**
* Normalize an angle in a 2&pi wide interval around a center value.
*
This method has three main uses:
diff --git a/src/test/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtilsTest.java b/src/test/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtilsTest.java
index 03acca200..d1707c5ef 100644
--- a/src/test/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtilsTest.java
+++ b/src/test/java/org/apache/commons/math/analysis/polynomials/PolynomialsUtilsTest.java
@@ -18,8 +18,8 @@ package org.apache.commons.math.analysis.polynomials;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.analysis.integration.LegendreGaussIntegrator;
+import org.apache.commons.math.util.ArithmeticsUtils;
import org.apache.commons.math.util.FastMath;
-import org.apache.commons.math.util.MathUtils;
import org.apache.commons.math.util.Precision;
import org.junit.Assert;
import org.junit.Test;
@@ -289,7 +289,7 @@ public class PolynomialsUtilsTest {
for (int w = 0; w < 10; ++w) {
for (int i = 0; i < 10; ++i) {
PolynomialFunction jacobi = PolynomialsUtils.createJacobiPolynomial(i, v, w);
- double binomial = MathUtils.binomialCoefficient(v + i, i);
+ double binomial = ArithmeticsUtils.binomialCoefficient(v + i, i);
Assert.assertTrue(Precision.equals(binomial, jacobi.value(1.0), 1));
}
}
diff --git a/src/test/java/org/apache/commons/math/linear/InverseHilbertMatrix.java b/src/test/java/org/apache/commons/math/linear/InverseHilbertMatrix.java
index 0c036d01c..bd68cab0d 100644
--- a/src/test/java/org/apache/commons/math/linear/InverseHilbertMatrix.java
+++ b/src/test/java/org/apache/commons/math/linear/InverseHilbertMatrix.java
@@ -17,7 +17,7 @@
package org.apache.commons.math.linear;
import org.apache.commons.math.exception.DimensionMismatchException;
-import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.util.ArithmeticsUtils;
/**
* This class implements inverses of Hilbert Matrices as
@@ -54,13 +54,13 @@ public class InverseHilbertMatrix
*/
public long getEntry(final int i, final int j) {
long val = i + j + 1;
- long aux = MathUtils.binomialCoefficient(n + i, n - j - 1);
- val = MathUtils.mulAndCheck(val, aux);
- aux = MathUtils.binomialCoefficient(n + j, n - i - 1);
- val = MathUtils.mulAndCheck(val, aux);
- aux = MathUtils.binomialCoefficient(i + j, i);
- val = MathUtils.mulAndCheck(val, aux);
- val = MathUtils.mulAndCheck(val, aux);
+ long aux = ArithmeticsUtils.binomialCoefficient(n + i, n - j - 1);
+ val = ArithmeticsUtils.mulAndCheck(val, aux);
+ aux = ArithmeticsUtils.binomialCoefficient(n + j, n - i - 1);
+ val = ArithmeticsUtils.mulAndCheck(val, aux);
+ aux = ArithmeticsUtils.binomialCoefficient(i + j, i);
+ val = ArithmeticsUtils.mulAndCheck(val, aux);
+ val = ArithmeticsUtils.mulAndCheck(val, aux);
return ((i + j) & 1) == 0 ? val : -val;
}
diff --git a/src/test/java/org/apache/commons/math/util/ArithmeticsUtilsTest.java b/src/test/java/org/apache/commons/math/util/ArithmeticsUtilsTest.java
index b4eba0fde..1753d7d6d 100644
--- a/src/test/java/org/apache/commons/math/util/ArithmeticsUtilsTest.java
+++ b/src/test/java/org/apache/commons/math/util/ArithmeticsUtilsTest.java
@@ -17,6 +17,9 @@
package org.apache.commons.math.util;
import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
import org.apache.commons.math.exception.MathArithmeticException;
import org.apache.commons.math.exception.MathIllegalArgumentException;
@@ -30,15 +33,16 @@ import org.junit.Test;
* @version $Id$
*/
public class ArithmeticsUtilsTest {
- /**
- * Exact direct multiplication implementation to test against
- */
- private long factorial(int n) {
- long result = 1;
- for (int i = 2; i <= n; i++) {
- result *= i;
- }
- return result;
+
+ /** cached binomial coefficients */
+ private static final List