commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Phil Steitz <p...@steitz.com>
Subject Re: svn commit: r1381206 - /commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/ArithmeticUtils.java
Date Wed, 05 Sep 2012 16:25:08 GMT
On 9/5/12 7:46 AM, erans@apache.org wrote:
> 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.

The mod operator is used below.  Javadoc should be modified to
reflect the implementation.  The algorithm appears to be a modified
version of what is in the reference.  Modifications should  be
called out if possible.  I haven't fully understood the
modifications, but it looks like there is some preprocessing of
arguments that speeds up the algorithm.  It would be great to
describe how it works or provide an external reference.  If we can't
do that, we should just say it is a "modified version" of binary gcd
and drop the statement about not using mod.  I was going to suggest
the that unmodified binary gcd in the private method below for
positive arguments be exposed for usage when both arguments are
known positive, but IIUC what is going on, even in that case it will
be faster to call the exposed method.  Do I have that right?

Phil
> +     * 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;
>      }
>  
>      /**
>
>
>


-- 


---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org


Mime
View raw message