commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject svn commit: r1381206 - /commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java
Date Wed, 05 Sep 2012 14:46:59 GMT
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;
     }
 
     /**



Mime
View raw message