Here's a variant using the even-odd speedup:
(Feel free to add my S-o-b if you use it.)
/*
* This implements Brent & Kung's "even-odd" variant of the binary GCD
* algorithm. (Often attributed to Stein, but as Knith has noted, appears
* a first-century Chinese math text.)
*
* First, find the common powers of 2 from a and b. Here, they are not
* divided by that common power, but rather we form a bitmap r which
* covers the least-significant bit set in either.
* Let s = (r+1)/2 be the common factor of 2.
*
* Then, shift down both a and b until both are odd multiples of s.
* Assuming a > b, then gcd(a, b) = gcd(a-b, b). But a-b is an even
* multiple of s and we can divide it by two at least once.
*
* Further, either a-b or a+b is a multiple of 4s, so by choosing between
* the two we can divde by 4. This is done by (in the case that a > b):
* 1. Compute (a - b) / 2
* 2. If this is an odd multiple of s (AND with r is non-zero), then
* add b to make an even multiiple of s. (This cannot overflow,
* as it equals (a + b) / 2, which must be less than max(a, b).)
* 3. Divide by 2 one or more more times until we are left with an odd
* multiple of r.
* 4. The result replaces a, and we go back to finding the larger.
*
* While the algorothm is equally correct without step2 and the first
* divide by 2, the benefit of adding it it can be evaluated with a CMOV
* rather than a full conditional branch.
*/
unsigned long bgcd2(unsigned long a, unsigned long b)
{
unsigned long r = a | b; /* To find common powers of 2 */
if (!a || !b)
return r;
r ^= r-1; /* Only the msbit of r (r/2 + 1) really matters */
if (!(a & r))
goto a_even;
if (!(b & r))
goto b_even;
while (a != b) {
/* (a & r) == (b & r) == (r >> 1) != 0 */
if (a > b) {
a = (a - b) >> 1;
if (a & r)
a += b;
a_even: do
a >>= 1;
while (!(a & r));
} else {
b = (b - a) >> 1;
if (b & r)
b += a;
b_even: do
b >>= 1;
while (!(b & r));
}
}
return b;
}
--
To unsubscribe from this list: send the line "unsubscribe linux-kernel" in
the body of a message to majord...@vger.kernel.org
More majordomo info at http://vger.kernel.org/majordomo-info.html
Please read the FAQ at http://www.tux.org/lkml/
