Felix Breuer's Blog

Visualizing the GCD

There are many beautiful geometric interpretations of the GCD and visualizations of the Euclidean algorithm out there, see for example here, here and here. I have written about this myself some time ago. But I find that some of the pictures become much clearer when lifted into the third dimension.

[Note: I am using Sage Cells to produce the pictures in this post. This has the advantage that you, dear reader, can interact with the pictures by modifying the Sage code. However, this will probably not run everywhere. It certainly won’t work in most RSS readers…]

Here is the graph of the GCD! First a little static preview…

Graph of the GCD

… and now a proper interactive graph. Be sure to rotate the image to get a 3d impression. (Requires Java…)

def mygcd(a,b,depth=0): if a > b: return mygcd(a-b,b,depth+1) if a < b: return mygcd(a,b-a,depth+1) if a == b: return [a, depth] def drawgcd(boxsize=10): M = boxsize s = Set([QQ(p)/QQ(q) for p in xrange(1,M) for q in xrange(1,M)]) lines = [] for x in s: p = x.numerator() q = x.denominator() l = min(M/p,M/q) [gcd,depth] = mygcd(p,q) lines.append(line3d([(p,q,1),(l*p,l*q,l)], rgbcolor=Color(depth/M,1,0.8,space="hsv").rgb())) return reduce(lambda a,b: a+b, lines) drawgcd(10)

What precisely are we seeing here? The gcd is a function $\mathrm{gcd}:\mathbb{N}^2\rightarrow\mathbb{N}$. Its graph is a (countable) set of points in $\mathbb{N}^3$. This graph is precisely the set of integer points contained in one of the (countably many) lines indicated in the above picture. The value of the gcd is given by the vertical axis.

Now, what are those lines? Suppose $a$ and $b$ are two natural numbers that are relatively prime, i.e., $\gcd(a,b)=1$. Then we know that $\gcd(ka,kb)=k$ for any positive integer $k$. So all (positive) integer points on the line through the origin and $(a,b,1)$ are in the graph of the gcd. Conversely, for any numbers $x,y$ with $gcd(x,y)=z$, we know that $\gcd(\frac{x}{z},\frac{y}{z})=1$, that is, the point $(x,y,z)$ lies on the line through the origin and $(\frac{x}{z},\frac{y}{z},1)$. So by drawing all these lines, we can be sure we hit all the integer points in the graph of the gcd.

What I find very beautiful about this is that when looking at this picture from the right angle, you can see the binary tree structure of this set of lines. It is precisely this binary tree that we walk along when running the Euclidean algorithm. The lines in the above picture are color coded by their depth in this “Euclidean” binary tree.

The version of the Euclidean algorithm I like best this one. Suppose you want to find $\gcd(a,b)=1$. If $a>b$, compute $\gcd(a-b,b)$. If $a<b$, compute $\gcd(a,b-a)$. If $a=b$, we are done and the number $a=b$ is the gcd. At each point in the algorithm we make a choice, depending on whether $a>b$ or $a<b$. If we keep track of these choices by writing down a 1 in the former case and a 0 in the latter, then, for any pair (a,b), we get a word of 0s and 1s that describes a path through a binary tree. As we walk from one pair $(a,b)$ to the next in this fashion, we trace out a path in the plane: For each 1 we take a horizontal step and for each 0 we take a vertical step. This path can look as follows.

a = 31 b = 19 def gcd1(a,b): # returns the gcd, a list of choices and a drawing if a > b: [gcd, path, drawing] = gcd1(a-b,b) return [gcd, path+[1], drawing+point2d([a,b])+line2d([[a-b,b],[a,b]])] if a < b: [gcd, path, drawing] = gcd1(a,b-a) return [gcd, path+[0], drawing+point2d([a,b])+line2d([[a,b-a],[a,b]])] if a == b: return [a, [], point2d([a,b])] [gcd,path,drawing] = gcd1(a,b) print "The gcd of %d and %d is %d." % (a,b,gcd) print "During the computation the Euclidean algorithm, makes the following choices:" print path drawing

Now comes the twist. We can turn the Euclidean algorithm upside down. Instead of walking from the point $(a,b)$ to a point of the form $(g,g)$ and thus finding the gcd, we can also start from the point $(1,1)$ and walk to a point of the form $(\frac{a}{g},\frac{b}{g})$ to find the gcd $g$. (Why $(\frac{a}{g},\frac{b}{g})$? Well, $(\frac{a}{g},\frac{b}{g})$ is the primitive lattice vector in the direction of $(a,b)$. In terms of the graph at the top, $(\frac{a}{g},\frac{b}{g},1)$ is the starting point of the ray through $(a,b,\gcd(a,b))$.)

How do we get there? We start out by fixing the lattice basis with “left” basis vector $(0,1)$ and “right” basis vector $(1,0)$. Their sum $(1,1)$, which I call the “center”, is a primitive lattice vector. Now we ask ourselves: is the point $(a,b)$ to the left or to the right of the line through the center? If it is to the left of that line, we continue with the basis $(0,1)$ and $(1,1)$. If it is to the right of that line, we continue with the basis $(1,1)$ and $(1,0)$ and recurse: We take the sum of our two new basis vectors as the new center and ask if $(a,b)$ lies to the left or to the right of the line through the center, and continue accordingly. If $(a,b)$ lies on the line through the center, we are done: the gcd of $a$ and $b$ is the factor by which I have to multiply the center to get to $(a,b)$.

As we follow this algorithm, we draw the path from one center to the next, until we reach $(\frac{a}{g},\frac{b}{g})$. Then we get the following picture. It looks almost like a straight line from $(1,1)$ to our goal, but it is only straight whenever we take two consecutive left turns or two consecutive right turns. Whenever we alternate between right and left, we have a proper angle, it only gets difficult to see this angle the farther out we get. If we write down a sequence of 1s and 0s, corresponding to the left and right turns we take, we get the exact same sequence as with the standard Euclidean algorithm.

a = 31 b = 19 def gcd2(a,b,left=vector([0,1]),right=vector([1,0])): center = right + left if b/a < center[1]/center[0]: [gcd, path, drawing] = gcd2(a,b,center,right) return [gcd, path+[1], drawing+point2d(center)+line2d([center,right+center])] if b/a > center[1]/center[0]: [gcd, path, drawing] = gcd2(a,b,left,center) return [gcd, path+[0], drawing+point2d(center)+line2d([center,left+center])] if b/a == center[1]/center[0]: return [a/center[0], [], point2d(center)] [gcd,path,drawing] = gcd2(a,b) print "The gcd of %d and %d is %d." % (a,b,gcd) print "During the computation the reverse Euclidean algorithm, makes the following choices:" print path drawing

The obvious next step is to draw the complete binary trees given by the above two visualizations of the Euclidean algorithm (up to a certain maximum depth). Drawing the first tree gives a big mess, because many of the edges of the tree overlap. You can clean up the picture by removing the lines and drawing only the points: In the code below comment the first definition of thisdrawing and uncomment the second. This then gives you a picture of all primitive lattice points (pairs of relatively prime numbers) in the plane at depth at most maxdepth in the tree.

maxdepth = 7 def drawtreehelper(a,b,parenta,parentb,depth): thisdrawing = [point2d([a,b]),line2d([[parenta,parentb],[a,b]])] #thisdrawing = [point2d([a,b])] if depth == 0: return thisdrawing else: return thisdrawing + drawtreehelper(a+b,b,a,b,depth-1) + drawtreehelper(a,b+a,a,b,depth-1) def drawtree(depth): return [point2d([1,1])] + drawtreehelper(2,1,1,1,depth) + drawtreehelper(1,2,1,1,depth) reduce(lambda a,b:a+b,drawtree(maxdepth))

The second variant of the Euclidean algorithm gives a much nicer picture, revealing the fractal nature of the set of primitive lattice points. The nodes of the tree, i.e., the points drawn in the picture are the starting points of the rays plotted in the 3d graph of the gcd given at the beginning of this post. In this way, this drawing of the tree gives the precise structure of the binary tree intuitively visible in the three-dimensional graph. One thing you may want to try is to increase the maxdepth of the tree to 10. (Warning: this may take much longer to render!) Note how far out from the origin some of the points at depths 10 are.

maxdepth = 7 unitmatrix = matrix([[1,0],[0,1]]) lefttransform = matrix([[1,0],[1,1]]) righttransform = matrix([[1,1],[0,1]]) thevector = vector([1,1]) def drawtreehelper2(basis,parent,depth): #thisdrawing = [point2d([a,b]),line2d([[parenta,parentb],[a,b]])] thispoint = basis*thevector #print thispoint,parent thisdrawing = [point2d(thispoint),line2d([thispoint,parent])] if depth == 0: return thisdrawing else: left = matrix return thisdrawing + drawtreehelper2(basis*lefttransform,thispoint,depth-1) + drawtreehelper2(basis*righttransform,thispoint,depth-1) def drawtree2(depth): return [point2d(thevector)] + drawtreehelper2(lefttransform,thevector,depth) + drawtreehelper2(righttransform,thevector,depth) reduce(lambda a,b:a+b,drawtree2(maxdepth))

The reader may have observed that the binary tree described in this post gives us a very nice way of enumerating the rational numbers. Calkin and Wilf made this observation in the paper Recounting the Rationals, whence it is sometimes called the Calkin-Wilf tree, even though Euclid tree might be a better name as it is given by the Euclidean algorithm itself. As we have seen above, a look through the lens of lattice geometry tells us how to draw this tree in the plane and how to find it in the graph of the GCD.