Discovering Number Theory

Jeffrey J. Holt
University of Virginia

John W. Jones
Arizona State University

Copyright © 2001 by W. H. Freeman and Company

The Euclidean Algorithm and
Linear Diophantine Equations

Initializations

\(Off[General::spell]; 
Off[General::spell1]; \)

The Division Algorithm Revisited

Below is a statement of the division algorithm, as given in the Prelab section of this chapter:

Theorem (Division Algorithm)  Let a be an integer and b be a positive integer.  Then there exist unique integers q and r such that

a  =  b q  +  r	and	0 ≤ r < b.

Recall that at the end of the last chapter, you used a Mathematica function called divalg to compute q and r for a given a and b.  Here is the definition for divalg:

divalg[a_, b_] := 
      Module[{q, r}, q = Floor[a/b]; r = a - b q; Return[{q, r}]];

Here's how it works for a = 123456789 and b = 369:

divalg[123456789, 369]

Thus, for this choice of a and b, we have q = 334571 and r = 90.  In Prelab exercise 2, you were asked to compute gcd(a, b) and gcd(b, r) for different choices of a and b.  If we do this computation in this case, we get

\(a\  = \ 123456789; \ 
b\  = \ 369; 
r\  = \ 90; 
    
{GCD[a, b], GCD[b, r]}\)

We can save ourselves a little typing by combining the computation required to compute r with the computation to compute gcd(a, b) and gcd(b, r).  We do this by using the fact that divalg[a,b][[2]] will return just r.

a = 123456789;
b = 369;
{GCD[a, b], GCD[b, divalg[a, b][[2]]]}

Mathematica Note

In general, if you have a list such as mylist = {1, 2, 4, 8, 16}, you can select a particular element from the list by using double brackets: 

   mylist[[4]] = 8
   
For those with computer programming experience, this is the analogue of working with an array, where we would specify the nth entry of the array mylist with mylist[[n]].

Research Question 1

Compute the values of gcd(a, b) and gcd(b, r) as above for different pairs a and b of your choosing until you have enough data to form a conjecture concerning the relationship between gcd(a, b) and gcd(b, r).  As always, once you have formed your conjecture, try to prove it!

The Euclidean Algorithm 1

Now that you have completed Research Question 1, you know that for a pair of integers a and b, we have

gcd(a, b) = gcd(b, r),

where r is the remainder when a is divided by b.  (Sure, this gives away the conjecture for Research Question 1, but you already knew the conjecture anyway, right?)  Thus, for example, if  a = 123456789 and b = 369 (as in the preceding section), then r = 90, and we see that

gcd(123456789, 369) = gcd(369, 90).

It's clear that computing gcd(369, 90) will be easier than computing gcd(123456789, 369).  (This will be true using any of the methods discussed  in the previous chapter.)  But rather than trying to compute gcd(369, 90), why don't we make the problem even easier by repeating what we did above?  If now we think of a = 369 and b = 90, then we have r equal to

divalg[369, 90][[2]]

Thus it follows that gcd(369, 90) = gcd(90, 9).  Repeating the above, we  compute

divalg[90, 9][[2]]

It shouldn't be surprising that the remainder in this case is equal to 0.   After all, it's easy to see that 90/9 = 10.  Thus 9 is a divisor of 90, so that we have gcd(90, 9) = 9.   Finally, we just string all of this together to arrive at

\(TraditionalForm
    \`gcd(123456789, \ 369)\  = \ 
      \(gcd(369, \ 90)\  = \ \(gcd(90, \ 9)\  = \ 9. \)\)\)

The process illustrated is called the Euclidean Algorithm, and it is very efficient for computing gcds.  In fact, this is probably the method that the Mathematica command GCD uses for computing greatest common divisors.

Exercise 1

Repeat the above procedure to compute gcd(7920, 4536).

Linear Diophantine Equations 1

As stated in the Prelab section, one goal of this chapter is to develop a systematic method for finding all integer solutions x and y (when there are any) to the linear diophantine equation

a x  +  b y  =  c,

where a, b, and c are integer constants.  As you discovered when working on the Prelab exercises, for a given choice of a and b, equation (1) may have several solutions or possibly no solutions; it depends on the value of c.  If we have specific values for a and b, how can we figure out the values of c for which equation (1) will have a solution?  As a first step, let's look at the situation for a specific choice of a and b, say a = 6 and b = 4.  In this case, the above equation becomes

6  x  +  4  y  =  c,

and our question is: For what values of c is there a solution to this equation?  One way to approach this is to try plugging a bunch of different values for x and y into the left-hand side of (2), and see what we get out.  After all, any value that comes out must be a suitable value for c.  (Do you see why?)  Below is the Mathematica code to compute the values of 6x + 4y for each choice of x and y satisfying -3 ≤ x ≤ 5 and -3 ≤ y ≤ 5:

Flatten[Table[6 x  +  4 y,  {x, -3, 5}, {y, -3, 5}]]

Mathematica Note

The Flatten is there to remove extra sets of unnecessary braces.  If you really want to see what it does, remove the Flatten command and re-execute the cell.

On the basis of this list, it looks like c must be an even integer in order for equation (2) to have a solution.  We can also see that some values appear several times; this indicates that for a given value of c there may be lots of solutions.

If we think about it, if x = m and y = n are solutions to 6  x  +  4  y  =  c then x = -m and y = -n are solutions to 6  x  +  4  y  =  -c.  Therefore, if we know the positive values of c for which there are solutions to equation (2), then we also know the story for negative values of c.  (Of course, c = 0 is easy.  Quick, name choices for x and y such that 6  x  +  4  y  =  0.)  Thus, we can restrict ourselves to c > 0.

Below is a Mathematica function that will compute a x  +  b y for specified values of a and b by plugging in x and y satisfying -n ≤ x ≤ n and -n ≤ y ≤ n and then removing the repeated and nonpositive values:

\(\(findc[a_, b_, n_]\  := \ 
	
      Module[{clist}, 
						
        clist\  = \ 
          Table[a\ \ x\  + \ b\ \ y, \ {x, \(-n\), n}, {y, \(-n\), n}]; 

        						clist\  = \ Rest[Union[Abs[Flatten[clist]]]]; 

        						Return[clist]]; \)\)

Mathematica Note

The Union command behaves just as you would expect with sets: it removes duplicate entries.  Along the way, it also sorts the list.

Here we also see Abs, which takes the absolute value of a number.  In this case, it is applied to a list.  The end result is to take the absolute value of every entry in the list.  Mathematica tries to do something sensible whenever an operation is applied to a list.  The best way to understand what it will do is by  experimenting.

Finally, we encounter the Rest command.  It removes the first entry of a list.  Can you guess why we wanted to do that?  If not, make a copy of findc in a new cell, give it a new name, and remove Rest from the definition. See if you can spot the difference between the original function findc and the new one you created.

Here is what we get when we try findc out on our equation, which corresponds to a = 6 and b = 4:

findc[6, 4, 10]

As we observed above, it looks like c must be a multiple of 2 in order for equation (2) to have a solution.   Let's try a different equation, say 9  x  +  12  y  =  c.  What do we get in this case?

findc[9, 12, 10]

Hmm, this time it looks like we get solutions only if c is a multiple of 3.  How about 5  x  +  8  y  =  c ?

findc[5, 8, 10]

Research Question 2

Evaluate findc[a,b] for values of a and b of your choosing until you have enough data to fill in the blank at the end of the following conjecture:

"In order for a x + b y  =  c to have solutions, c must be of the form ______."

The Euclidean Algorithm 2

Finding One Solution

Your conjecture in Research Question 2 gives conditions that c must satisfy in order for

a x  +  b y  =  c

to have solutions.  Suppose that c does have the form given in your conjecture; in this case, are we guaranteed that there will be solutions to our equation?  This is a pretty general question.  Let's look at one special case.  Suppose that a = 408, b = 126, and c = gcd(408, 126) = 6.  Are there any solutions to the equation

408  x  +  126  y  =  6 ?

It turns out that the answer to this question is yes.  To find a solution, we start by going through the steps of the Euclidean Algorithm to show that gcd(408, 126) = 6:

divalg[408, 126]

Therefore, 408 = 3 · 126 + 30.  Remember that

gcd(408, 126) = gcd(126, 30).

Now we just repeat:

divalg[126, 30]

Therefore, 126 = 4 · 30 + 6.

divalg[30, 6]

Therefore, 30 = 5 · 6 + 0.  Since clearly 30 is divisible by 6, it follows that gcd(30, 6) = 6.  Thus it must be that gcd(408, 126) = 6.   Now, what does this have to do with solutions to

408  x  +  126  y  =  6 ?

Good question.  Let's take a look at the steps in the Euclidean Algorithm again:

Reorganizing the equation in step (b), we have

6  =  126 - 4 · 30.

From step (a), we see that 30 = 408  -  3 · 126.  Substituting into equation (3), we get

\(TraditionalForm
    \`6\ \  = \ 
      \(126\  - \ 4\ ·\ 30
	 = \ 
        \(126\  - \ 4\ ·\ \((408\  - \ 3\ ·\ 126)\)

          	 = \ \(\(-4\)\ ·\ 408\  + \ 
              \((1\  + \ 12)\)\ ·\ 126
	 = \ 
            \(-4\)\ ·\ 408\  + \ 13\ ·\ 126\)\)\)\)

Therefore, we see that x  =  -4, y = 13 is a solution to

408  x  +  126  y  =  6.

Just to be on the safe side, let's check using Mathematica by executing the code below.

408*(-4)  +  126*13  ==  6

Mathematica Note

The "==" means that this code is treated as an equation rather than a definition.  When you execute this code, Mathematica compares the values on both sides of the equation, and returns True if both sides are equal, and False otherwise.

The True means that the left and right sides of the equation are equal.  Let's look at another example.  Suppose that a = 1232 and b = 573, and we want to find a solution to

1232 x + 573 y = d,

where d = gcd(1232, 573).  First we compute d using the Euclidean Algorithm:

divalg[1232, 573]

divalg[573, 86]

divalg[86, 57]

divalg[57, 29]

divalg[29, 28]

divalg[28, 1]

We see that d = gcd(1232, 573) = 1, and so we are looking for a solution to

1232 x + 573 y = 1.

Here's a summary of the above steps:

Now we work our way back up the chain:

\(TraditionalForm
    \`1\ \  = \ \ 
      \(29\  - \ 1\ ·\ 28
	 = \ \ 
        \(29\  - \ 1\ ·\ \((57\  - \ 1\ ·\ 29)\)\ \ 

          	 = \ \ \(\(-\ 1\)\ ·\ 57\  + \ 2\ ·\ 29

            	 = \ \ \(\(-\ 1\)\ ·\ 57\  + \ 
                2\ ·\ \((86\  - \ 1\ ·\ 57)\)
	 = \ \ 
              \(2\ ·\ 86\  - \ 3\ ·\ 57
	 = \ \ 
                \(2\ ·\ 86\  - \ 
                    3\ ·\ \((573\  - \ 6\ ·\ 86)\)

                  	 = \ \ 
                  \(\(-\ 3\)\ ·\ 573\  + \ 20\ ·\ 86

                    	 = \ \ 
                    \(\(-\ 3\)\ ·\ 573\  + \ 
                        20\ ·\ 
                          \((1232\  - \ 2\ ·\ 573)\)
	 = \ \ 
                      20\ ·\ 1232\  - \ 
                        43\ ·\ 573. \)\)\)\)\)\)\)\)\)

So, x = 20 and y = -43 should be a solution to 1232  x  +  573  y  =  1.  As before, let's check:

20*1232  -  43*573   ==  1

Solving the linear equation

a x  +  b y  =  gcd(a,  b)

is useful in a variety of places in number theory.  Indeed, the fact that this equation always has a solution is so handy that we give it a special  name: The GCD Trick.  We'll have several opportunities to exploit the GCD Trick as the course progresses.  Mathematica has a built-in command that implements the procedure illustrated above.  Given a and b, the command ExtendedGCD[a, b] computes both gcd(a, b) and a pair (x, y) which satisfies the above equation.  For example, if we start with a = 1232 and b = 573 as above, we get:

ExtendedGCD[1232,  573]

The output comes in the form {gcd(a, b), {x, y}}.  Note that it gives the same answer as was found above.

Finding All Solutions

So far, so good.  The method that we used in the two previous examples is quite general, and will work to find a solution to any equation of the form

a x  +  b y  =  d,

where d = gcd(a, b).  Although we have made some progress towards our original goal of  finding all solutions to a x  +  b y  =  c, we still have a long way to go.  The next question that we shall  consider is the following: Are there any solutions to equation (4) other than the one guaranteed by the GCD Trick and found using the reverse Euclidean Algorithm method?  Let's go with a simple example,

7  x  +  2  y  =  1.

It's clear that gcd(7, 2) = 1, and we get the solution (Subscript[x, 0])  =  1, (Subscript[y, 0])  =  -3 by reversing the steps of the Euclidean Algorithm.  One way to look for other solutions is to use Mathematica to do a simple search.  Solving for y in terms of x, we get the equation

y = (1  -  7  x)/2.

If we plug in an integer value of x, and the corresponding value of y is also an integer, then the pair (x,  y) is a solution to our equation.  (Be sure that you understand why before moving forward.)  Here's the Mathematica code to do the work.  In this case, we are just using values of x satisfying -10 ≤ x ≤ 10.

Do[If[IntegerQ[(1 - 7 x)/2], Print[{x, (1 - 7 x)/2}]], {x, 
        -10, 10}]

Mathematica Note

IntegerQ[t] returns True if t is an integer, and False otherwise.

As you can see, our solution (Subscript[x, 0]) = 1, (Subscript[y, 0])  =  -3 is among those found.  Let's try a different equation, say

8  x  +  3  y  =  1.

In this case, it's easy to see that gcd(8, 3) = 1.  The Euclidean Algorithm yields the solution (Subscript[x, 0])  =  -1, (Subscript[y, 0])  =  3.  Solving for y in terms of x, we have y  =  (1  -  8  x)/3.  Here's the test:

Do[If[IntegerQ[ (1 - 8 x)/3], Print[{x, (1 - 8 x)/3}]], {x, 
        -10, 10}]

As you might guess, we are working towards forming a conjecture, and it will be handy to have a function that will automatically look for solutions.  Here's one provided for your convenience:

sols[a_, b_, d_, n_] := 
      Do[If[IntegerQ[(d - a x)/b], Print[{x, (d - a x)/b}]], {x, 
          -n, n}]

The function sols will search for solutions to a x  +  b y  =  d, trying values of x satisfying -n ≤ x ≤ n.  Let's check this function against the preceding output:

sols[8, 3, 1, 10]

Looks fine.  Here's one more equation: 11  x  +  4  y  =  1.  As we can see, gcd(11, 4) = 1, and the Euclidean Algorithm produces the solution (Subscript[x, 0])  =  -1, (Subscript[y, 0])  =  3.  Here's the result of a search:

sols[11, 4, 1, 20]

Research Question 3

Suppose that gcd(a, b) = 1, and that ((Subscript[x, 0]), (Subscript[y, 0])) is the solution to 

				a x  +  b y  =  1

produced by the Euclidean Algorithm.  Find a general formula for all solutions (x,  y) to the above equation, giving x in terms of (Subscript[x, 0]) and y in terms of (Subscript[y, 0]).  Be sure to prove that your formula works and that it accounts for all solutions.

Research Question 4

Now suppose that gcd(a, b) = d > 1, and that ((Subscript[x, 0]), (Subscript[y, 0])) is the solution to 

				a x  +  b y  =  d			

produced by the Euclidean Algorithm.  Find a general formula for all solutions (x,  y) to the above equation, giving x in terms of (Subscript[x, 0]) and y in terms of (Subscript[y, 0]).  

Hint:  The equations a x  +  b y  = d and (a/d) x + (b/d) y = 1 have the same set of solutions.

Linear Diophantine Equations 2

Now we're getting somewhere.  With the research questions in the preceding section complete, we now know all of the solutions to the equation

a x  +  b y  =  d,

where d = gcd(a, b).  Therefore, it remains to determine all of the solutions to the  equation

a x  +  b y  =  k d,

where k is an integer.  To get a feel for what is going on, let's look at an  example.  In the preceding section, we looked at the equation 7  x  +  2  y  =  1.  Let's change this a bit, say to

7  x  +  2  y  =  5.

Using the Euclidean Algorithm, we found that x  =  1, y  =  -3 is a solution to 7  x  +  2  y  =  1.  Thus, it is easy to see that x = 1·5 = 5 and  y = -3·5 = -15 is a solution to the equation 7  x  +  2  y  =  5.  What about the other solutions?  Let's put sols to work to find some others:

sols[7, 2, 5, 10]

Research Question 5

Suppose that gcd(a, b) = d and that ((Subscript[x, 0]), (Subscript[y, 0])) is a solution to a x + b y  =  d .   Find the general form of all solutions (x, y) to

				a x + b y = k d,		

giving x in terms of (Subscript[x, 0]) and y in terms of (Subscript[y, 0]).