GCD Examples

The following programs implement the Euclidean algorithm and some of its variations to compute the greatest common divisor of two integers a and b.   The last example includes some
 code to investigate the frequency of occurrence of a given number as a quotient
 in the Euclidean algorithm.

Simple GCD

Let (a,b) denote the greatest common divisor of  a and b. It has the following properties.
 
   (a,a)= |a|
   
   (-a,b)= (a,b) 
   
   (a,b)= (a-b,b)
   
   These two are sufficient to  write a  program to compute GCD.  If the two numbers in the input are not  equal, then subtract the smaller number from the larger and compute the GCD of the two smaller numbers.

gcd[a_,b_]:= If[ a==b, 
               (* then *) Abs[a],
               (*otherwise *) If[ a> b, 
                               (*then*) gcd[a-b,b], 
                               (*otherwise*) gcd[b-a,a]
                                ]
               ]
(*

gcd[13,39]

gcd[ -729, 90]

The program  has a few problems. First, it will fail if a is negative and b is positive, because the values keep increasing (Try it! )  One solution is to take the Absolute value in the inner gcd computation, or check different cases depending on the signs of a and b. Another defect is that if one of the values is 0, then the  program will not terminate.  Both of these are rectified in the following program.  The first part gives a definition of the gcd when one of the input values is 0.

newgcd[a_,0] :=  If[ a==0, (*then*)undefined, 
                           (*else*)Abs[a]];
newgcd[0,a_] :=  If[ a==0, undefined, Abs[a]];

newgcd[a_,b_]:= If[ a==b, 
               (* then *) Abs[a],
               (*otherwise *) If[ a> b, 
                               (*then*) gcd[a-b, Abs[b]], 
                               (*otherwise*) gcd[b-a,Abs[a]]
                                ]
               ]
(*

newgcd[2,0]

newgcd[0,0]

newgcd[-729,90]

Euclidean Algorithm

The Euclidean algorithm  is based on the following property.

 (a,b) = (b,r), where a= bq +r and r= a mod b.
 
 This  is used to give the following definition of gcd.  If a>b, then we compute the remainder, otherwise we  can interchange the variables  using (a,b)= (b,a).

EuclideanGCD[a_,b_]:= If[ b==0, (*then*) If[a==0, undefined, Abs[a]],
                          (*otherwise*)  If[a>b, EuclideanGCD[Mod[a,b],Abs[b]],
                                                 EuclideanGCD[b,a]]]
(*

EuclideanGCD[12,36]

EuclideanGCD[13,0]

EuclideanGCD[-13, -39]

EuclideanGCD[-13,65]

EuclideanGCD[0,0]

Exercise: Implement the binary gcd algorithm in Mathematica.  The binary gcd method 
uses only division by two and  subtraction, so should be more efficient than the 
Euclidean method when implemented in a language such as C, where it is possible to make use of these machine level operations.

Complete GCD

We can modify the procedure to print quotient and the  remainder in the intermediate steps or to return  a list of all the quotients and remainders in each step.  First, we use the convention that the remainder  r= a mod b  satisfies  0<= r < |b|. This is not the convention adopted by Mathematica, its Mod function gives a remainder which is the same sign as b.  This is easy to correct with a new Mod function follows. This also requires a modification of the quotient function.

NewMod[a_,b_]:= If[b>0, Mod[a,b], 
                      Mod[a,b] -b]

Mod[27,-13]

NewMod[27,-13]

NewQuotient[a_,b_]:= If[ b>0, Floor[a/b],
                          Ceiling[a/b]]

NewQuotient[27,-13]

The Euclidean algorithm applied to a and b takes the form

 a= bq1 + r1,
 b= r1q2 + r2
 r1=r2 q3 + r3
   .
   .
   .
   The first step should give us { q1, r1} and then we apply the method to { b, r1}. The procedure is the following.

CompleteGCD[ a_,b_]:=CompleteGCD[a,b]=
        If[ NewMod[a,b]==0, {{NewQuotient[a,b],0}},
              PrependTo[CompleteGCD[b,NewMod[a,b]],
                         {NewQuotient[a,b], NewMod[a,b]}
                       ]
           ]

CompleteGCD[678,21]

CompleteGCD[751,22]

Some explanation is in order regarding the program. The program  constructs a list of pairs {q,r} which occur in the Euclidean  algorithm.  At the last step,  b divides and a and the procedure stops. Otherwise we compute  q and r and add it to the front of the list obtained by  applying the algorithm to   b and r.   In the definition of the program, the  additional statement CompleteGCD[a,b] in the first line tells Mathematica  remember the values of the function.

Exercise: The completegcd program should give the reader an idea to rewrite  the EuclideanGCD program in a slightly better way using the NewMod function.

The quotients are the first elements of each pair. This can extracted
 using the Transpose and First commands.

QuotientsInGCD[a_,b_]:= First[Transpose[CompleteGCD[a,b]]]

QuotientsInGCD[751,22]

An Experiment

If a and b are any two random integers, one would expect that any quotient can occur with equal likelihood in the Euclidean algorithm.  Let us conduct an experiment. We select    a few random pairs and  compute the quotients in the Euclidean algorithm applied to all these pairs.   We count the  relative proportion of 1's occuring in these quotients. The experiment is performed several times to see the  average  value of the percentage of 1's occuring as quotients in the Euclidean algorithm.  

 First, we make a list of 10 random pairs.

pairs= Table[{t=Random[Integer, 25000], Random[Integer,t]},
                {i,1,10}]

Then apply the QuotientsInGCD procedure to each of the pair  using the MapThread command.

allquotients = MapThread[QuotientsInGCD, Transpose[pairs]]

Next, we count the relative proportion of 1's in  all the quotients we obtained.

quotients=Flatten[allquotients];
 N[Count[quotients,1]/Length[quotients]]

Remarkably,  the proportion of 1's in the quotient is quite large.  We can
 combine all the operations to do the calculations for more pairs. The semicolon's at the end of the intermediate statements suppress all the intermediate computations and only the final  probability is shown.

pairs= Table[{t=Random[Integer, {10000,100000}], Random[Integer,{10000,t}]},
                {i,1,100}];
allquotients = MapThread[QuotientsInGCD, Transpose[pairs]];
 quotients=Flatten[allquotients];
 N[Count[quotients,1]/Length[quotients]]

Exercise: Repeat this exercise several times. What is the probability that 1 is the quotient. Similarly, find the probability the a quotient is 2 in the Euclidean algorithm.