(*^ ::[ frontEndVersion = "Microsoft Windows Mathematica Notebook Front End Version 2.2"; microsoftWindowsStandardFontEncoding; fontset = title, "Times", 24, L0, center, nohscroll, bold; fontset = subtitle, "Times", 18, L0, center, nohscroll, bold; fontset = subsubtitle, "Times", 14, L0, center, nohscroll, italic; fontset = section, "Times", 18, L0, nohscroll, bold, grayBox; fontset = subsection, "Times", 14, L0, nohscroll, bold, blackBox; fontset = subsubsection, "Times", 12, L0, nohscroll, bold, whiteBox; fontset = text, "Times", 12, L0, nohscroll; fontset = smalltext, "Times", 10, L0, nohscroll; fontset = input, "Courier", 12, L-5, nowordwrap, bold; fontset = output, "Courier", 12, L-5, nowordwrap; fontset = message, "Courier", 12, L-5, nowordwrap, R65535; fontset = print, "Courier", 12, L-5, nowordwrap; fontset = info, "Courier", 12, L-5, nowordwrap, B65535; fontset = postscript, "Courier", 12, L0, nowordwrap; fontset = name, "Geneva", 10, L0, nohscroll, italic; fontset = header, "Times", 12, L0; fontset = footer, "Times", 12, L0, center; fontset = help, "Times", 10, L0, nohscroll; fontset = clipboard, "Times", 12, L0, nohscroll; fontset = completions, "Times", 12, L0, nohscroll; fontset = graphics, "Courier New", 10, L0, nowordwrap, nohscroll; fontset = special1, "Times", 12, L0, nohscroll; fontset = special2, "Times", 12, L0, nohscroll; fontset = special3, "Times", 12, L0, nohscroll; fontset = special4, "Times", 12, L0, nohscroll; fontset = special5, "Times", 12, L0, nohscroll; fontset = leftheader, "Times", 12, L2; fontset = leftfooter, "Times", 12, L2; fontset = reserved1, "Courier New", 10, L0, nowordwrap, nohscroll;] :[font = title; inactive; preserveAspect; startGroup; nohscroll; center; ] GCD Examples :[font = text; inactive; preserveAspect; nohscroll; ] 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. :[font = section; inactive; preserveAspect; startGroup; Cclosed; nohscroll; ] Simple GCD :[font = text; inactive; preserveAspect; nohscroll; ] 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. :[font = input; initialization; preserveAspect; startGroup; nowordwrap; ] *) gcd[a_,b_]:= If[ a==b, (* then *) Abs[a], (*otherwise *) If[ a> b, (*then*) gcd[a-b,b], (*otherwise*) gcd[b-a,a] ] ] (* (* :[font = message; inactive; preserveAspect; endGroup; nowordwrap; ] General::spell1: Possible spelling error: new symbol name "gcd" is similar to existing symbol "GCD". :[font = input; preserveAspect; startGroup; nowordwrap; ] gcd[13,39] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] 13 ;[o] 13 :[font = input; preserveAspect; startGroup; nowordwrap; ] gcd[ -729, 90] :[font = message; inactive; preserveAspect; nowordwrap; ] $IterationLimit::itlim: Iteration limit of 4096 exceeded. :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] Hold[If[995175 == -729, Abs[995175], If[995175 > -729, gcd[995175 - -729, -729], gcd[-729 - 995175, 995175]]]] ;[o] Hold[If[995175 == -729, Abs[995175], If[995175 > -729, gcd[995175 - -729, -729], gcd[-729 - 995175, 995175]]]] :[font = text; inactive; preserveAspect; nohscroll; ] 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. :[font = input; initialization; preserveAspect; nowordwrap; ] *) 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]] ] ] (* (* :[font = input; preserveAspect; startGroup; nowordwrap; ] newgcd[2,0] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] 2 ;[o] 2 :[font = input; preserveAspect; startGroup; nowordwrap; ] newgcd[0,0] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] undefined ;[o] undefined :[font = input; preserveAspect; startGroup; nowordwrap; ] newgcd[-729,90] :[font = output; inactive; formatted; output; preserveAspect; endGroup; endGroup; nowordwrap; ] 9 ;[o] 9 :[font = section; inactive; preserveAspect; startGroup; startGroup; nohscroll; ] Euclidean Algorithm :[font = text; inactive; preserveAspect; nohscroll; ] 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). :[font = input; initialization; preserveAspect; nowordwrap; ] *) 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]]] (* (* :[font = input; preserveAspect; startGroup; Cclosed; nowordwrap; ] EuclideanGCD[12,36] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] 12 ;[o] 12 :[font = input; preserveAspect; startGroup; nowordwrap; ] EuclideanGCD[13,0] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] 13 ;[o] 13 :[font = input; preserveAspect; startGroup; nowordwrap; ] EuclideanGCD[-13, -39] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] 13 ;[o] 13 :[font = input; preserveAspect; startGroup; nowordwrap; ] EuclideanGCD[-13,65] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] 13 ;[o] 13 :[font = input; preserveAspect; startGroup; nowordwrap; ] EuclideanGCD[0,0] :[font = output; inactive; formatted; output; preserveAspect; endGroup; endGroup; nowordwrap; ] undefined ;[o] undefined :[font = text; inactive; endGroup; nohscroll; ] 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. :[font = section; inactive; preserveAspect; startGroup; Cclosed; nohscroll; ] Complete GCD :[font = text; inactive; preserveAspect; nohscroll; ] 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. :[font = input; initialization; preserveAspect; nowordwrap; ] *) NewMod[a_,b_]:= If[b>0, Mod[a,b], Mod[a,b] -b] (* :[font = input; preserveAspect; startGroup; nowordwrap; ] Mod[27,-13] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] -12 ;[o] -12 :[font = input; preserveAspect; startGroup; nowordwrap; ] NewMod[27,-13] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] 1 ;[o] 1 :[font = input; initialization; preserveAspect; nowordwrap; ] *) NewQuotient[a_,b_]:= If[ b>0, Floor[a/b], Ceiling[a/b]] (* :[font = input; preserveAspect; startGroup; nowordwrap; ] NewQuotient[27,-13] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] -2 ;[o] -2 :[font = text; inactive; preserveAspect; nohscroll; ] 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. :[font = input; initialization; preserveAspect; nowordwrap; ] *) 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]} ] ] (* :[font = input; preserveAspect; startGroup; nowordwrap; ] CompleteGCD[678,21] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] {{32, 6}, {3, 3}, {2, 0}} ;[o] {{32, 6}, {3, 3}, {2, 0}} :[font = input; preserveAspect; startGroup; nowordwrap; ] CompleteGCD[751,22] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] {{34, 3}, {7, 1}, {11, 1}, {20, 1}, {6, 1}, {5, 1}, {89, 1}, {4, 1}, {3, 1}, {2, 1}, {1, 1}, {16, 1}, {-8, 1}, {3, 0}} ;[o] {{34, 3}, {7, 1}, {11, 1}, {20, 1}, {6, 1}, {5, 1}, {89, 1}, {4, 1}, {3, 1}, {2, 1}, {1, 1}, {16, 1}, {-8, 1}, {3, 0}} :[font = text; inactive; preserveAspect; nohscroll; ] 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. :[font = text; inactive; preserveAspect; nohscroll; backColorRed = 65535; backColorGreen = 65535; backColorBlue = 0; fontColorRed = 0; fontColorGreen = 0; fontColorBlue = 0; plain; fontName = "Times"; fontSize = 12; ] Exercise: The completegcd program should give the reader an idea to rewrite the EuclideanGCD program in a slightly better way using the NewMod function. :[font = text; inactive; preserveAspect; nohscroll; ] The quotients are the first elements of each pair. This can extracted using the Transpose and First commands. :[font = input; initialization; preserveAspect; nowordwrap; ] *) QuotientsInGCD[a_,b_]:= First[Transpose[CompleteGCD[a,b]]] (* :[font = input; preserveAspect; startGroup; nowordwrap; ] QuotientsInGCD[751,22] :[font = output; inactive; formatted; output; preserveAspect; endGroup; endGroup; nowordwrap; ] {34, 7, 11, 20, 6, 5, 89, 4, 3, 2, 1, 16, -8, 3} ;[o] {34, 7, 11, 20, 6, 5, 89, 4, 3, 2, 1, 16, -8, 3} :[font = section; inactive; preserveAspect; startGroup; Cclosed; nohscroll; ] An Experiment :[font = text; inactive; preserveAspect; nohscroll; ] 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. :[font = input; preserveAspect; startGroup; nowordwrap; ] pairs= Table[{t=Random[Integer, 25000], Random[Integer,t]}, {i,1,10}] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] {{4507, 494}, {17792, 8221}, {1331, 506}, {21453, 18208}, {9466, 2824}, {12375, 11471}, {2800, 567}, {15682, 14599}, {19859, 3489}, {19390, 6666}} ;[o] {{4507, 494}, {17792, 8221}, {1331, 506}, {21453, 18208}, {9466, 2824}, {12375, 11471}, {2800, 567}, {15682, 14599}, {19859, 3489}, {19390, 6666}} :[font = text; inactive; preserveAspect; nohscroll; ] Then apply the QuotientsInGCD procedure to each of the pair using the MapThread command. :[font = input; preserveAspect; startGroup; nowordwrap; ] allquotients = MapThread[QuotientsInGCD, Transpose[pairs]] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] {{9, 8, 10, 6, 28, 1, 2, 5, 6}, {2, 6, 11, 6, 2, 14, 1, 1, 5, 2, 12, 2, 10, 3, 5, 1, 2}, {2, 1, 1, 1, 2, 1, 2, 2}, {1, 5, 1, 1, 1, 1, 3, 180}, {3, 2, 1, 5, 3, 2, 3, 1, 3}, {1, 12, 1, 2, 4, 1, 1, 1, 1, 2, 5, 1, 2}, {4, 1, 15, 5}, {1, 13, 2, 12, 10, 3, 1, 29, 1, 16, -8, 3}, {5, 1, 2, 4, 13, 1, 3, 8, 2, 10, 3, 5, 1, 2}, {2, 1, 9, 1, 26, 4, 1, 1, 1, 3}} ;[o] {{9, 8, 10, 6, 28, 1, 2, 5, 6}, {2, 6, 11, 6, 2, 14, 1, 1, 5, 2, 12, 2, 10, 3, 5, 1, 2}, {2, 1, 1, 1, 2, 1, 2, 2}, {1, 5, 1, 1, 1, 1, 3, 180}, {3, 2, 1, 5, 3, 2, 3, 1, 3}, {1, 12, 1, 2, 4, 1, 1, 1, 1, 2, 5, 1, 2}, {4, 1, 15, 5}, {1, 13, 2, 12, 10, 3, 1, 29, 1, 16, -8, 3}, {5, 1, 2, 4, 13, 1, 3, 8, 2, 10, 3, 5, 1, 2}, {2, 1, 9, 1, 26, 4, 1, 1, 1, 3}} :[font = text; inactive; preserveAspect; nohscroll; ] Next, we count the relative proportion of 1's in all the quotients we obtained. :[font = input; preserveAspect; startGroup; nowordwrap; ] quotients=Flatten[allquotients]; N[Count[quotients,1]/Length[quotients]] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] 0.3269230769230769 ;[o] 0.326923 :[font = text; inactive; preserveAspect; nohscroll; ] 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. :[font = input; preserveAspect; startGroup; nowordwrap; ] 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]] :[font = output; inactive; formatted; output; preserveAspect; endGroup; nowordwrap; ] 0.3245227606461087 ;[o] 0.324523 :[font = text; inactive; preserveAspect; endGroup; endGroup; nohscroll; backColorRed = 65535; backColorGreen = 65535; backColorBlue = 0; fontColorRed = 0; fontColorGreen = 0; fontColorBlue = 0; plain; fontName = "Times"; fontSize = 12; ] 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. ^*)