Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!usc!wuarchive!psuvax1!psuvm!cxt105
From: CXT105@psuvm.psu.edu (Christopher Tate)
Newsgroups: comp.lang.postscript
Subject: Re: GCD algorithm
Message-ID: <90305.105824CXT105@psuvm.psu.edu>
Date: 1 Nov 90 15:58:24 GMT
Organization: Penn State University
Lines: 83
Supersedes: <90305.103320CXT105@psuvm.psu.edu>

In article <1990Oct31.182705.3184@light.uucp>, bvs@light.uucp (Bakul Shah) says:

>- To justify posting this to a PostScript group, here is a faster
>  version of basically the same algorithm:
>
>    /gcd {              % given m, n on the stack
>      dup 3 1 roll mod  % leaves n, m mod n on the stack
>      dup 0 eq {
>        pop             % leaves n on the stack when done
>      } {
>        gcd             % find gcd of the new pair
>      } ifelse
>    } bind def
>
>  Not too surprisingly, this recursive version is faster than
>  the one below which uses an explicit loop:
>
>    /gcd {
>      { dup 3 1 roll mod dup 0 eq { pop exit } if } loop
>    } bind def
>
>                     [stuff deleted]
>
>  So....  Who is right?  How did two different alg. came to be
>  known as Euclid's?  Does the faster alg. also have such a nice
>  generalization for gcd of N numbers?  The slower alg. can be
>  modified slightly to simulteneously compute LCM.  Can the same
>  thing be done with the faster one?  Time to visit a library!

As it happens, I'm takinga computational number theory course even as
we speak :-), so here's what my text has to say on the subject of Euclid's
GCD algorithm....

Essentially, the whole thing works because of this odd property of the GCD;
namely, that  gcd(a, b) = m*a + n*b  for some integers m and n (not
necessarily unique).  From this, we can do some fiddling with the properties
of division and so on to show that  gcd(a, b) = gcd(b, a - m*b)  for ANY
integer m.  In Euclid's algorithm, we essentially take m to be as large as
possible, i.e. a - m*b is as close to zero as possible.  Since we can write
r = a mod b = a - m*b (for m as large as possible), we simply use a mod b
for the next iteration, and try to calculate gcd(b, a mod b).

I once ran across a version of this algorithm (in C) which looked like
this:

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

Essentially, taking the smallest steps down (using a-b instead of a mod b)
possible at a time, but avoiding any division.  Still used an awful lot of
stack space sometimes, though....

I'm very suprised that the recursive version is faster than the iterative
one; I suspect it has to do with the way recursion is implemented in
PostScript.  Is it still faster if you don't bind the gcd procedures?  How
about if you give it two large, almost-equal primes?  In most programming
languages (i.e. native code compilation, etc.), given an iterative algorithm
and an equivalent recursive algorithm, the iterative algorithm will be faster
if only because it won't have to go through all the overhead of the recursive
function calls....

There's a version of this algorithm known as the binary GCD algorithm,
which works on the binary representation of the numbers (indirectly) instead
of on the numbers themselves.  This allows it to take advantage of the fact
that just finding a modulus is computationally icky (as it requires division),
while dividing an integer by two is what computers do best.  It's not usually
faster than the above, though, unless implemented directly in machine code.

And, for everyone who's been following this, the version of the GCD algorithm
shown above IS much faster than mine -- use it instead!  Mine was a quick-
and-dirty rewrite of a Rexx version of the algorithm, and so was not at all
suited to the real stack-handling capabilities of PostScript.

-------
Christopher Tate                   | "Pardon me, but is this your bar
                                   |    of soap?"
cxt105@psuvm.psu.edu               | "Why, yes, I suppose it is...."
{...}!psuvax1!psuvm.bitnet!cxt105  |
cxt105@psuvm.bitnet                | "So do we."