Forth for Enlightenment, part four of ten, GCD

GCD implements an interesting way to compute Greatest Common Divisor. You remember fractions? Let's say you need to reduce the fraction

210
---
462

They're both even, so you can divide by 2.

210   105
--- = ---
462   231

Uhm... the digits of 105 add to 6, and the digits of 231 add to 6, so they're both divisible by 3.

210   105   35
--- = --- = --
462   231   77

Then, 35 is 5 * 7, and 77 is 7 * 11, so they're both divisible by 7.

210   105   35   5
--- = --- = -- = -
462   231   77   11

So, effectively, 210 / 42 = 5, and 462 / 42 = 11. 42 is the greatest common divisor between 210 and 464. The way this was solved was to attempt to find prime numbers that are factors. And it turned out that it wasn't that difficult. That's because we were able to divide big numbers by small numbers until all the numbers were pretty small. But you might have felt a looming fear that you'd have to do a bigger divide. There are lots of ways to get the greatest common divisor. Astonishingly, there is one that doesn't involve division. In C, it looks something like this:

int gcd(int a, int b) {
  if (a == b) return a;
  if (a > b) {
    return gcd(a - b, a);
  } else {
    return gcd(a, b - a);
  }
}

Let's use this to figure out the GCD of 210 and 462. Since b > a, the first step is to subtract 462 - 210 = 252, so it's gcd(210, 252). Then, b > a, so gcd(210, 252 - 210 = 42). Then gcd(168, 42). Then gcd(126, 42). Then gcd(84, 42). Then gcd(42, 42). Since 42 = 42, the answer is 42. OK, maybe this isn't so easy to do in your head. But there used to be computer processors like the 8080 that did not know how to divide directly. Such processors could use this sort of thing, and they were really relatively fast at it.

Anyway, the recursive version of this function on the HP-28 looks like this.

« DUP2 IF == THEN
      DROP ; return a
   ELSE
      DUP2 IF > THEN
         SWAP DUP ROT - SWAP GCD ; return gcd(a-b, a)
      ELSE
         SWAP DUP ROT SWAP - GCD ; return gcd(a, b-a)
      END
   END
» 'GCD' STO

This GCD function is called with two numbers, which i'll call a and b on the stack. The first thing the function does is make copies of the two arguments. They are then compared. If they are the same as each other, one of them is dropped from the stack. This is the return value of the function. If a and b aren't the same, the if a > b it returns the value from the recursive call of gcd(a - b, a). Let's look at this with a stack trace.

Level	Start	SWAP	DUP	ROT	-	SWAP
3			b	a
2	a	b	a	a	a	a - b
1	b	a	a	b	a - b	a

And then GCD is called. The stack trace for when b > a follows.

Level	Start	SWAP	DUP	ROT	SWAP	-
3			b	a	a
2	a	b	a	a	b	a
1	b	a	a	b	a	b - a

This program works. However, for large arguments it can run out of memory. That is, despite not putting much in the way of data on the stack, it runs out of recursive call depth. For example, it chokes on gcd(100, 101). I wanted to know exactly how deep it could go, and i tried various number pairs that are separated by 1. And gcd(100, 101) failed. So i wrote this little program to find out when it first fails. It starts at 50, and tries up to 200. It then calls gcd(x + 1, x). And, unexpectedly, it returned that gcd(180,179) works. That was really unexpected. It turns out that gcd(179,180) fails. The order of the arguments is important, for no apparent reason.

« 50 200 FOR x
      x
      x 1 + x GCD
      DROP DROP
   NEXT
» 'G' STO

Anyway, the fact that fairly small arguments lead to a stack crash is not good. The C version, using the full -O3 optimization of the gcc compiler, never runs out of memory. That's because this compiler optimizes the tail recursion, and it never generates much of a call stack. Likewise, the Scheme version never runs out of memory. The Scheme language requires that tail recursion optimization is implemented. Perl 6 is supposed to do some simple tail recursion optimization. All my code is in Perl 5, and it fails pretty quickly. It complains of Deep recursion with gcd(101, 100). By gcd(131, 130) it runs out of stack memory. But could the program be written without recursion? Sure. In fact, it's pretty easy. First we'll see what it might look like in C. Often to eliminate tail recursion, all you need is a jump back to the top.

int gcd(int a, int b) {
  while (a != b) {
    if (a > b) {
      a = a - b;
    } else {
      b = b - a;
    }
  }
  return a;
}

What this does is loop until a and b are the same. It changes it's own arguments each loop instead of calling itself. If this is done on the HP-28, then there's no recursion, and large numbers can be handled, given time. Note that the ≠ in the first line is a single HP-28 character, a shifted '='.

« DUP2 WHILE ≠ REPEAT
      DUP2 IF > THEN
         SWAP DUP ROT - SWAP ; repeat with gcd(a-b, a)
      ELSE
         SWAP DUP ROT SWAP - ; repeat with gcd(a, b-a)
      END
      DUP2 ; for the != in WHILE
   END
   DROP ; drop one of the equal a, b.
» 'GCD2' STO

This version is also a little quicker, because the while loop is quicker than a function call. At first, things seem to be moved around at random compared with the original version. The DROP that removed one of the identical values on the stack moved to after the while loop, just as the C version's return statement moved. Otherwise, the first IF became a WHILE.

Is the iterative version easier or harder to follow than the recursive version? This is an important maintenance question. And it isn't clear that either one are very easy. After all, even having stepped through an example, it's not clear how the subtractions end up giving you what amounts to division. And the begining of the answer is that repeated subtraction can be used in division. Think about how you have to guess the next digit in long division. You don't really have to guess. You can try the digit one, do a subtraction of the divisor, and repeat this (while incrementing the answer digit) until just before it goes negative. But that's just a hint. On the other hand, the nonrecursive version demonstrates that nothing needs to be remembered and that the only time something is returned is when a == b.