MATH-841
Performance improvement in method "gcd(int, int)" (~2 to 4 times faster than the previous implementation). Thanks to Sebastien Riou. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1381206 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Gilles Sadowski
parent
350f726ce5
commit
68fe06032a
|
|
@ -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.
|
||||
}
|
||||
if (k == 31) {
|
||||
throw new MathArithmeticException(LocalizedFormats.GCD_OVERFLOW_32_BITS,
|
||||
p, q);
|
||||
}
|
||||
// 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
|
||||
}
|
||||
// B5 [reset max(u,v)]
|
||||
if (t > 0) {
|
||||
u = -t;
|
||||
|
||||
long al = a;
|
||||
long bl = b;
|
||||
boolean useLong = false;
|
||||
if (a < 0) {
|
||||
if(Integer.MIN_VALUE == a) {
|
||||
useLong = true;
|
||||
} else {
|
||||
v = t;
|
||||
a = -a;
|
||||
}
|
||||
// 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
|
||||
al = -al;
|
||||
}
|
||||
if (b < 0) {
|
||||
if (Integer.MIN_VALUE == b) {
|
||||
useLong = true;
|
||||
} else {
|
||||
b = -b;
|
||||
}
|
||||
bl = -bl;
|
||||
}
|
||||
if (useLong) {
|
||||
if(al == bl) {
|
||||
throw new MathArithmeticException(LocalizedFormats.GCD_OVERFLOW_32_BITS,
|
||||
p, q);
|
||||
}
|
||||
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;
|
||||
}
|
||||
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;
|
||||
}
|
||||
|
||||
/**
|
||||
|
|
|
|||
Loading…
Reference in New Issue