Author: erans
Date: Wed Sep 5 14:46:59 2012
New Revision: 1381206
URL: http://svn.apache.org/viewvc?rev=1381206&view=rev
Log:
MATH-841
Performance improvement in method "gcd(int, int)" (~2 to 4 times faster than
the previous implementation). Thanks to Sebastien Riou.
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java?rev=1381206&r1=1381205&r2=1381206&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java
Wed Sep 5 14:46:59 2012
@@ -358,25 +358,25 @@ public final class ArithmeticUtils {
}
/**
- * <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>
+ * Computes 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.
+ * 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>The invocation {@code gcd(0, 0)} is the only one which returns
- * {@code 0}.</li>
+ * <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.
@@ -386,62 +386,118 @@ public final class ArithmeticUtils {
* a non-negative {@code int} value.
* @since 1.1
*/
- public static int gcd(final int p, final int q) {
- int u = p;
- int v = q;
- if ((u == 0) || (v == 0)) {
- if ((u == Integer.MIN_VALUE) || (v == Integer.MIN_VALUE)) {
+ 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(u) + FastMath.abs(v);
+ return FastMath.abs(a + b);
}
- // keep u and v negative, as negative integers range down to
- // -2^31, while positive numbers can only be as large as 2^31-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 < 31) { // while u and v are
- // both even...
- u /= 2;
- v /= 2;
- k++; // cast out twos.
+
+ 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 (k == 31) {
- throw new
MathArithmeticException(LocalizedFormats.GCD_OVERFLOW_32_BITS,
- p, q);
+ if (b < 0) {
+ if (Integer.MIN_VALUE == b) {
+ useLong = true;
+ } else {
+ b = -b;
+ }
+ bl = -bl;
}
- // B2. Initialize: u and v have been divided by 2^k and at least
- // one is odd.
- int 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
+ if (useLong) {
+ if(al == bl) {
+ throw new
MathArithmeticException(LocalizedFormats.GCD_OVERFLOW_32_BITS,
+ p, q);
}
- // B5 [reset max(u,v)]
- if (t > 0) {
- u = -t;
- } else {
- v = t;
+ 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;
}
- // 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 * (1 << k); // gcd is u*2^k
+ blbu = bl;
+
+ // Now "al" and "bl" fit in an "int".
+ b = (int) al;
+ a = (int) (blbu % al);
+ }
+
+ return gcdPositive(a, b);
+ }
+
+ /**
+ * 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;
}
/**
