As others have said, this is Stein's algorithm, and is remarkable for having taken so long to be discovered, considering the importance of GCD to number theory. While it was discovered at least as early as 1962, it was not published until 1967.
One thing that hasn't been mentioned is that it scales much better than division. As you deal with larger and larger numbers, division time grows quadratically. Shift and subtraction operations, however, do not. On an AMD Zen processor, a 64-bit integer divide takes 13-44 cycles, whereas every Stein's algorithm op takes one cycle, with the exception of the bit-scan op, which takes two.
Also note that in the example you cite, many things are done in loops that we now have instructions for - specifically the loop that continually divides the number by two until it's odd: that's just two instructions now.
Also note that while on x86, dividion and modulo are computed and returned simultaneously, on other architectures, a separate multiplication and subtraction operation is needed in order to compute the modulo, further slowing down the traditional GCD calculation.
Modern Stein's algorithm code looks like this. (Note __builtin_ctzll, count trailing zeroes, is the gcc built-in to tell you how many times you need to divide a number by two until it's odd.)
uint64_t gcd(uint64_t x, uint64_t y) {
if (x == 0) return y;
if (y == 0) return x;
auto xs = __builtin_ctzll(x);
x >>= xs;
auto ys = __builtin_ctzll(y);
y >>= ys;
auto topshift = std::min(xs, ys);
// x and y are both odd and we don't know which is bigger.
for (;;) {
if (x < y) std::swap(x, y);
x -= y;
if (x == 0) break;
x >>= __builtin_ctzll(x);
}
return y << topshift;
}
See https://godbolt.org/z/7hKa34 for the generated code, and note how complex the "simple" divide algorithm is, at least on x86-64.
On my MacBook Pro just now, computing GCD of all pairs of numbers between 10,000 and 20,000 takes a little over 2 seconds using Stein, and about 8 seconds using Euclid.
Finally, note that Euclid is worst when the incoming numbers have a ratio close to the golden ratio; the first several divisions produce simply 1, and the so the first few modulo operations are merely extremely expensive subtraction operations, while the number becomes smaller very slowly.
By contrast, Euclid is fastest, and Stein the worst, when the numbers have largely differing magnitudes; Euclid will immediately bring the larger number into the same range as the small one, while Stein only clears away one or two bits per iteration.
You can alleviate these problem areas by combining the algorithms like so:
// The "Brown-Euclid-Stein Hybrid", if no one has thought of it.
uint64_t besh(uint64_t x, uint64_t y) {
if (x == 0) return y;
if (y == 0) return x;
auto xs = __builtin_ctzll(x);
x >>= xs;
auto ys = __builtin_ctzll(y);
y >>= ys;
auto topshift = std::min(xs, ys);
if (x > y) std::swap(x, y);
for (;;) {
y %= x;
y >>= __builtin_ctzll(y);
if (y <= 1) return (y ? 1 : x) << topshift;
x %= y;
x >>= __builtin_ctzll(x);
if (x <= 1) return (x ? 1 : y) << topshift;
}
}
u%vwould actually do something like__divide_int32(u,v).remainder!