36

First of all, I hope this is the right place to post this question. I was looking through one of my dad's old programming books from the 80s, and at the back it has a list of utility functions, one of which is gcd. Its implementation is as follows:

gcd(u, v)
int u, v;
{
    int k, t, f;
    if ((u <= 0) || (v <= 0)) return 1;
    k = 0; f = 1;
    while ((0 == (u % 2)) && (0 == (v % 2))) {
        k++; u >>= 1; v >>= 1; f *= 2;
    }
    if (u & 1) {
        t = -v; goto B4;
    } else {
        t = u;
    }
B3: if ( t > 0) { t >>= 1; } else { t = -((-t) >> 1); };
B4: if (0 == (t%2)) goto B3;
    if (t > 0) u = t; else v = -t;
    if (0 != (t = u - v)) goto B3;
    return (u * f);
}

It's clearly meant to be compatible with an old version of C, given it assumes an int return type and the location it declares the type of the parameters. And from my quick testing (and conversion to newer C), it seems to work. But why is it so complicated? Euclid's Algorithm has been known for literally thousands of years, and is just

gcd(int u, int v) {
    if (v == 0) return u;
    return gcd(v, u % v);
}

Assuming recursion would be a problem, this could also easily be converted to an iterative version. Are there any technical reasons that the above algorithm would be preferred?

14
  • 8
    Perhaps "it's so complicated" because the programmer wasn't very good at programming? The layout is definitely spaghetti. We knew better in the 1970s. Dec 22 '20 at 0:16
  • 16
    You assume they had a reasonably fast % operator.
    – user253751
    Dec 22 '20 at 13:43
  • 6
    On un-powerful CPUs, there's a good chance there's no division instruction at all, and your compiler uses a division algorithm someone has written in software in the standard library - u%v would actually do something like __divide_int32(u,v).remainder!
    – user253751
    Dec 22 '20 at 15:16
  • 7
    It looks like it may have been assembly or Fortran code translated to C.
    – Barmar
    Dec 22 '20 at 15:56
  • 4
    Why just say "an old programming book" and not name the book and the author? It might lead to better answers. And it's only polite, anyway.
    – TRiG
    Dec 23 '20 at 21:02
38

This is (or at least appears to be) the binary euclidean GCD. https://en.wikipedia.org/wiki/Binary_GCD_algorithm has a slighty different version (but the version you post has the nice feature that there's no recursion (which I seem to recall earlier C compilers wanted to limit). As I noted in my comment that I deleted when I realized I do have access to Wikipedia even if I don't have access to my crypto text book that covered various GCD algorithms, I didn't think 80's C compilers would do the "optimization" that modulus by a power of 2 can be done using a bitwise and of that number -1 (e.g. N % 2 == N & 1, N % 4 == N & 3), but that would be the main benefit (since division is much slower than bitwise operators, especially in 80's / 90's hardware)

Edit fwiw, I was hoping I could run a test on an Apple IIGS emulator but alas while the Orca C compiler is open sourced, to actually get the binary requires buying it and I wasn't that interested at this point. I instead ran it on QEMU running 386 emulator using OpenWatcom with the -0 flag which is supposed to generate code that will run on 8086+ ... and lo and behold, the recursive one was about 50% faster to calculate gcd(i,j) for i and j between 1 and 1000 (and I ran some code with openwatcom that did foo % 2 vs foo & 1 and they both took the same amount of time so while I don't know how to get openwatcom to dump assembly to compare, I rather expect it to be doing the optimization (I did try to go through quickly the source of Orca C to see what that would spit out for a 6502 / 65816, but wasn't able to find what I was looking for quickly enough)

8
  • 13
    @CalvinGodfrey Yes, possibly. Many 16 and 32 bit machines did not have a hardware division instruction and would have to run a little subroutine for every division. Even for CPUs that did have hardware division, it could take hundreds of cycles while rotates would only take a few. Today most general-purpose CPUs take a single cycle for both rotates and division.
    – RETRAC
    Dec 21 '20 at 23:04
  • 16
    On hardware without integer divide (or modulus) instructions it would probably be much faster. Note that the only modulus calculations are "x % 2" which were presumably written that way for readability and would be optimized by the compiler. Otherwise they could all have been rewritten as "x & 1" tests, which was actually done in one line of the code.
    – alephzero
    Dec 21 '20 at 23:13
  • 5
    @RETRAC I don't believe modern architectures have single-cycle latency division. They might have pipelines with a throughput of one division per cycle, but that's no reason to stop avoiding divisions.
    – Nobody
    Dec 22 '20 at 9:21
  • 2
    @Nobody Yes, they don't have single cycle latency. Nothing has single cycle latency due to the pipelining.
    – RETRAC
    Dec 22 '20 at 15:47
  • 11
    Shifting takes 1 (amortized) cycle, but on the latest Intel I9s, integer division takes 10 cycles. A few years ago it was ~30 cycles. It has never been 1. Source (page 286) Dec 22 '20 at 17:59
20

The reason is that it's quicker.

You just need to look up GCD in Knuth, that is Volume 2 of TAOCP, to get answers and analysis. There you have comparisons of "Original Euclidean algorithm", "Modern Euclidean algorithm" and "Binary gcd algorithm". The last one is introduced with the remark that "Since Euclid's patriarchal algorithm has been used for so many centuries, it is a rather surprising fact that it may not always be the best method for finding the greatest common divisor after all." This new algorithm was discovered by J. Stein in 1961, and has the advantage that it has no division, but only (i) subtraction, (ii) testing whether a number is even or odd, and (iii) shifting the binary representation of an even number to the right (halving)".

The code you are reading follows the how "Algorithm B" is described there (Vol. 2, p. 321), and the labels B3 and B4 reference those steps in the algorithm. I think it's sensible to use labels that directly correspond to the steps as named by Knuth, since it makes it easy to compare the program with something known to work and thoroughly analyzed. Using that instead of writing your own (maybe naive) GCD isn't that far off from using a standard library function in a higher-level language, like math.gcd in Python. It is there, and you can assume it will be done in an efficient way.

With the make-up processor Knuth analyzes these algorithms for, the program that implements this binary algorithm is with some assumptions on average 8.8n + 5 cycles where the average value for Euclid is 11.1n + 7 (with random inputs between 1 and 2^n).

That is made with a variant of MIX, which is Knuth's original make-up processor, similar to real processors from the 1960s. How it fares with more modern processors (for example Knuth's replacement MMIX) I don't know. But at the Wikipedia page about the binary algorithm there is a reference from 2000 proving that proved that binary GCD uses about 60% fewer bit operations than the Euclidean algorithm.

1
  • 1
    The B3 and B4 labels immediately made me think that it's a literal translation of an algorithm as presented in TAOCP.
    – texdr.aft
    Dec 25 '20 at 0:01
13

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;
  }
}
2
  • 2
    How does you hybrid fare on the benchmark Stein and Euclid scored 2 s and 8 s on?
    – equaeghe
    Dec 24 '20 at 14:06
  • 1
    @equaeghe Great question. When I tried it, the hybrid was around 4s. Both numbers being between 10,000 and 20,000 is an unnaturally good fit for Stein, though. When I tried larger numbers, the hybrid was more competitive... until I got above 2^32, when the slow speed of 64-bit divide made it suddenly much less competitive.
    – jorgbrown
    Jan 4 at 10:58
7

It avoids recursion, so the resources required (in particular stack depth) are independent of the parameters. OP might find that the book didn't assume the availability of virtual memory, or targeted CP/M or real-mode x86 systems with a stack limit of 64K.

How does Knuth deal with this in his 1970s books?

4
  • 7
    While it avoids recursion, that is not the main point of the algorithm. The main point is already explained in an another answer, it avoids the use of dividing numbers to get modulo of decimal numbers by using binary arithmetic operations and bit shifts.
    – Justme
    Dec 22 '20 at 8:19
  • 1
    I don't believe you can say "not the main point" without more info from OP. I was highlighting the most glaring difference between the "by the book" code and his recursive alternative, and particularly since I explicitly mentioned CP/M etc. assumed that the lack of division and other complex operations would be assumed. Dec 22 '20 at 8:36
  • 2
    What I gather is that the OP did cover that a standard recursive gcd could be converted to iterative gcd if recursive execution was a problem, and was still wondering about the algorithm itself in the last sentence, instead of wondering the algorithm being recursive or iterative.
    – Justme
    Dec 22 '20 at 8:42
  • 1
    Knuth wrote on it extensively in TAOCP, at least in modern editions. As noted above, it is a binary GCD, there is also a Lehmer GCD, for example. See here for some implementations. Dec 23 '20 at 19:58
6

As others have pointed out, this GCD algorithm avoids the use of division and modulus (except for division and modulus by 2, which can be implemented very easily in binary arithmetic). There might be good reasons to avoid division and modulus in some environments.

Let's take a closer look at the code. The algorithm is basically in two pieces.

Top part: find the highest power of 2 which divides both u and v. So you divide both variables by 2 until one of them is odd, and count the number of times you did that.

At the end of the top part, you know that one of the two variables is odd (the other one might still be even). The code keeps track of the one that we know for sure is odd. And then, for the one that might be even, we might as well divide out all its remaining factors of 2 to make it odd also.

Bottom part: repeatedly subtract either u from v or vice-versa. Naturally you want to know which variable is which, so that you can update the right one. And also, when you subtract an odd number from an odd number, the result is even. So for that variable, you might as well divide out the factors of 2 again!

Comment: in both parts, the code keeps track of which variable is which, by setting t to be either positive or negative. This is a confusing way to do it.

Comment: I am not convinced that there is any compelling reason to remember which variable is (known to be) odd, in either part. You can always just divide out the 2's in both variables. The one that's odd won't take long.

Comment: They replace all the /2 operations with >>1, but they leave %2 everywhere instead of &1. Maybe this makes sense depending on your compiler.

Comment: As noted in the comments, k is never used, so it could be removed. Alternatively, you could remove all instances of f and write return (u << k); at the end.

Comment: As noted in the comments, this gives the wrong answer if either input is 0. Also, I'd say it gives the wrong answer for negative inputs. Personally, I would say that the greatest common divisor of 6 and -16 is 2. It definitely isn't 1.

Comment: So, a lot of the confusion comes from this "use the sign of t to encode a boolean" idea. And also, I don't think that boolean value is all that important to begin with. And also, we could probably do without t altogether.

Bottom line: I don't think the code had to be as confusing as it was.

For reference, here's a Python version of this algorithm, with some changes as described above. It should be straightforward to port to C, just put in enough brackets and semicolons, etc.

def gcd(u, v):
    if u < 0: u = -u
    if v < 0: v = -v
    # handle negative input before zero input,
    # to ensure output is nonnegative in cases like gcd(0, -1)
    if u == 0: return v
    if v == 0: return u
    # Top part: find exponent on 2 in output
    k = 0
    while ((0 == u&1) and (0 == v&1)):
        k += 1; u >>= 1; v >>= 1
    # Bottom part: subtract u from v or vice-versa
    while True:
        while u&1 == 0:
            u >>= 1
        while v&1 == 0:
            v >>= 1
        if u > v:
            u -= v
        elif u < v:
            v -= u
        else: # u == v, DONE.
            return u << k
3

Other answers/comments have touched on two very good reasons: Avoiding recursion (stack overflow) and avoiding division/modulo which may take orders of magnitude more CPU cycles than the messy-looking alternative.

If the book is 1980s you've got to consider it was likely targeting very early microprocessors/minicomputers - 1970s hardware - where you have incredibly limited resources by modern standards.

Recursive code can eat resources very quickly and the cost in time and memory could be relatively huge by the standards of the time.

Also remember that a lot of optimisation kung-fu that's well understood (or even done behind your back by the compiler) these days was hardly thought of back then, C was a relatively new thing and a lot of coding was done in more arcane languages or assembler.

The fact the code contains gotos makes it pretty darn old-school considering it was considered harmful since 1968.

The book you've got may be 1980s but is the 1st printing of the book earlier? It could easily be a 10-year-old text.

9
  • 1
    Please remember that "considered harmful" was not Dijkstra's choice of title: it was Wirth's and judging by Wirth's early code was more than a little tongue in cheek. Dec 22 '20 at 17:23
  • True, but he does make a good point and it's held true for ~50 years so I don't think it was too extreme. Looking at the state of modern code I can only imagine how bad it would be if we were throwing goto in everywhere with wild abandon. Anyway, the title isn't important here - it's the effect - and the fact that goto rapidly fell from favour helps date the text of the book.
    – John U
    Dec 23 '20 at 8:59
  • I agree up to a point, but would remind you that Dijkstra himself became unhappy about the extent to which avoidance of goto became treated as a religious tenet (see the archive of his papers). Dec 23 '20 at 9:04
  • As I said, for the sole purpose of dating a piece of code containing goto it doesn't matter what anyone thinks of it.
    – John U
    Dec 23 '20 at 9:12
  • Note that if you know that both numbers are positive, then you can replace the modulo with a pair of subtractions. en.wikipedia.org/wiki/…
    – CSM
    Dec 23 '20 at 16:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.