



Problem You need to find the greatest common divisor (GCD) or the lowest common multiple (LCM) of two numbers. The GCD of two numbers a and b is the greatest number that divides evenly into both a and b. For example the GCD of 24 and 60 is 12.
Solution Both of these problems and many others can be solved with Euclid’s algorithm. We will begin by explaining how to obtain GCD and then show that LCM is a simple extension. To compute the GCD naively we could start from the smallest of the two numbers and work our way downwards until we find a number that divides into both of them:
for (int i=Math.min(a,b); i>=1; i)
if (a%i==0 && b%i==0)
return i;
Although this method is fast enough for most applications, a much faster approach is to use Euclid's algorithm. Euclid's algorithm iterates over the two numbers until a remainder of 0 is found. The algorithm is based on the principle that the GCD of two numbers does not change if the smaller number is subtracted from the larger number. For example, suppose we want to find the GCD of 2336 and 1314. We begin by expressing the larger number (2336) in terms of the smaller number (1314) plus a remainder:
2336 = 1314 x 1 + 1022
We now do the same with 1314 and 1022:
1314 = 1022 x 1 + 292
We continue this process until we reach a remainder of 0:
1022 = 292 x 3 + 146
292 = 146 x 2 + 0
The last nonzero remainder is the GCD. So the GCD of 2336 and 1314 is 146. This algorithm can be easily coded as a recursive function:
//assume that a and b are positive
public long GCD(long a, long b)
{
if (b==0) return a;
return GCD(b,a%b);
}
Using this algorithm we can also find the lowest common multiple (LCM) of two numbers. For example the LCM of 6 and 9 is 18 since 18 is the smallest number that divides both 6 and 9. Notice that 18=(6*9)/3. In general we can express LCM(a,b) as (a*b)/GCD(a,b). To prevent potential overflow we perform division before multiplication in the implementation:
public long LCM(long a, long b)
{
return b/GCD(a,b)*a;
}
Discussion Euclid’s algorithm can compute the GCD of large numbers efficiently. In 1844 Gabriel Lamé showed that it never requires more iterations than five times the number of digits h (in base 10) of the smaller number. Since the computational expense of each iteration has order h, the overall complexity grows as h^2.
The algorithm has many theoretical and practical applications. It is an important element of the RSA algorithm. It can also be used to solve linear Diophantine equations: equations with integer coefficients of the form ax + by = c. To solve such equations we need a variation of the Euclid's Algorithm, where we record the coefficients floor(a/b) at each step. We can then use these coefficients to find integers c1 and d1 such that b*c1  a*d1 = 1. If c does not divide evenly by GCD(a,b) then there are no integer solutions; otherwise there are an infinite number of integer solutions and one of these is x0 = (+/)c*d1, y0 = (+/)c*c1. Note that the signs depend on the sign of b*c1  a*d1, as well as, the signs of a, b and c. Once we have a single solution (x0,y0) then the general solution takes the form x = x0 + b/GCD(a,b)*t, y = y0  a/GCD(a,b)*t, where t is an arbitrary integer. The following code implements these idea:
// Find solutions to a linear diophantine equation ax+by=c
// where a, b and c are integers
public void solveDiophantine(long a, long b, long c)
{
System.out.println("Attempting to solve: "+a+"x + "+b+"y = "+c);
//special case
if (a==0 && b==0 && c==0)
{
System.out.println("x and y can take any value");
return;
}
long a1=Math.abs(a);
long b1=Math.abs(b);
long c1=0, c2=1;
long d1=1, d2=0;
while(b1>0)
{
long mult=a1/b1;
long temp=c2;
c2=c2*mult+c1;
c1=temp;
temp=d2;
d2=d2*mult+d1;
d1=temp;
temp=b1;
b1=a1%b1;
a1=temp;
}
long gcd=a1;
if (gcd==0  Math.abs(c)%gcd!=0)
{
System.out.println("There are no integer solutions!");
return;
}
long det=c1*d2c2*d1;
long signa=a>0 ? 1: 1;
long signb=b>0 ? 1: 1;
long x=d1*c*det*signa/gcd;
long y=c1*c*det*signb/gcd;
System.out.println("One solution is x="+x+", y="+y);
System.out.println("A general solution is x="+x+" + "+b/gcd+"t, y="+y+"  "+a/gcd+"t");
}
For a more detailed description of the above algorithm, see the following: http://mathworld.wolfram.com/DiophantineEquation.html http://mathforum.org/library/drmath/view/51595.html http://www.wikihow.com/SolveaLinearDiophantineEquation
Now let's try solving CommonMultiples (Div I 250 SRM 346). We are asked to count how many numbers between lower and upper inclusive are multiples of all numbers in the set a. We begin by finding the LCM of all numbers in a. It can be shown that LCM of many numbers can be expressed as the LCM of individual pairs, eg., LCM(a, b, c, d) = LCM(a, LCM(b, LCM(c,d))). We find the LCM of all numbers in the first for loop and store it in lcm. We have to be careful, because lcm can overflow an int. To avoid overflow we modify our GCD and LCM functions to return longs. Now, the count of numbers not greater than upper having the required property is upper/lcm (integer division). Similarly, the count of numbers less than lower having the required property is (lower1)/lcm. Taking the difference of these two results gives us the final answer:
public class CommonMultiples
{
public int countCommMult(int[] a, int lower, int upper)
{
long lcm=1;
for (int i=0; i<a.length; i++) lcm=LCM(lcm,a[i]);
return (int)(upper/lcm(lower1)/lcm);
}
public long GCD(long a, long b)
{
if (b==0) return a;
return GCD(b,a%b);
}
public long LCM(long a, long b)
{
return b/GCD(a,b)*a;
}
}
Here are some more practice problems to try:
SRM 358 Div I 500 BalanceScale SRM 427 Div II 500 DesignCalendar
EDIT2: Added a comment to the gcd() that a,b>0. Thanks to syg96 EDIT3: Fixed some bugs in the diophantine code. Should be correct now. 

This is pretty much based on my Mathematics for Topcoders article, with a few minor modifications. I am sure there are many other problems that use this algorithm, but I couldn't find them. 

Good. Would be even better to expand "Discussion" with the actual algorithm of solving Diophantine equations and examples of two problems solved using GCD. 

I can certainly give an algorithm for solving Diophantine equations, but I am not sure if I can do a good job at explaining the mathematics behind it. I don't want the cookbook to be too mathematical. What do you think? 

It's Cookbook on programming, not on math, so I think it's fine just to describe the algorithm without its math background (maybe add a "See also" link to some place where it's explained in detail). 

First of all, gcd can be written in one line:
int gcd(int a, int b) { return (b ? gcd(b, a%b) : abs(a)); }
This code works well for any signed integers except the pair (0, 0).
Second, I feel that this statement is wrong: "Since the computational expense of each iteration has order h, the overall complexity grows as h^2." Actually, each iteration requires division of hdigit numbers so it requires O(h^2) time instead of O(h). In this way you can only prove the overall O(h^3) complexity. If the division is implemented in a proper way, it can be proved that the algorithm requires O(h^2) time. But the proof is rather difficult and not that obvious. So let's change this statement just to "With proper division implementation the computational complexity for large numbers is O(h^2)" or something like that.
If we speak of large number gcd, we should include binary algorithm as it is much easier to implement for large numbers. Here is the code for ints:
int gcd(int a, int b) {
if (a == 0) return b;
if (b == 0) return a;
if (a%2==0 && b%2==0) return gcd(a/2, b/2) * 2;
else if (a%2 == 0) return gcd(a/2, b);
else if (b%2 == 0) return gcd(a, b/2);
else return (a>=b ? gcd(ab, b) : gcd(a, ba));
}
Works for nonnegative numbers only (except zero pair as always). Time is O(h^2).
I think we should make the extended Euclid algorithm clear:
int Euclid(int a, int b, int &x, int &y) {
if (b == 0) {
x = (a<0 ? 1 : 1);
y = 0;
return a*x;
}
int q = a/b;
int r = a%b;
int d = Euclid(b, r, y, x);
y = q*x;
return d;
}
...
int x, y;
g = Euclid(a, b, x, y);
This code returns g = gcd(a, b) and also finds such x and y that: a*x + b*y = g It works for any pair of signed numbers (except zero pair). Moreover, the solution it finds is almost minimal in absolute value: x <= b/g y <= a/g This experimental fact often allows us not to be afraid of overflow.
The solution for general diophantine equation includes a lot of cases such as a==0 or b==0 of c==0 plus signs of a and b etc. Do we really need to include it? The extended Euclid algorithm code is sufficient to find reciprocal element modulo M which is the most common usage. 

This code works well for any signed integers except the pair (0, 0). What's wrong with the pair (0,0)? gcd(0,0) is 0. 

Yes, the algo will give zero. If you think it is correct, then the algorithm works for zero pair too=) Speaking strictly mathematically, gcd(0, 0) is not defined. That's why I said it doesn't work for zero pair. 

Thank you for your input. Here are some comments:
* I avoided the oneline gcd() definition for clarity purposes. * For the complexity I used the information from Wikipedia. Perhaps Wikipedia is wrong, or I misinterpreted. * Thanks for the large number gcd. It will be nice to include it. * I don't understand your extended euclid's algorithm, but I am happy to include it. * I think my solution to diophantine equations should handle +/ of a and b. I am not sure that it deals with all 0 cases correctly. Please let me know if you find a bug in it. Personally, I feel that it is not necessary to include it, because it appears quite infrequently in problems. 

Well, your gcd can return negative number on negative input. It can be a bad surprise for someone=) Better add "return abs(a);" or write a,b>=0 as input constraint or at least mention that your gcd can return negative result.
Wikipedia is right on complexity. And the proof about O(h^2) is correct there: it assumes that division is O(h(l+1)) and calculates the sum...
Well, I should have written the explanation of extended Euclid. We want to find any solution (x, y) of equation: a*x + b*y = gcd(a, b);
Two cases: 1. Let b == 0. Then x = sgn(a); y = 0; is a solution: a*sgn(a) + b*0 = a*sgn(a) = abs(a) = gcd(a, 0); Moreover, we choose positive gcd here. 2. Let b != 0. Then Divide a/b with reminder: int q = a / b; (quotient) int r = a % b; (reminder) Then it is guaranteed (even for signed a,b !) that: a = q*b + r; Ok, now let's call Euclid(b, r, y, x) to solve the smaller problem recursively. It returns gcd(b, r) = gcd(a, b) and fills numbers y and x such that: b*y + r*x = d (note that x and y are swapped) Now use r = (aq*b): b*y + (aq*b)*x = d Now group everything around a and b: d = b*y + (aq*b)*x = b*y + a*x  q*b*x = a*x + b*(yq*x) This equation tells us that: nx = x; ny = y  q*x; is a solution for problem a*nx + b*ny = gcd(a, b); So the problem is solved.
In case I forget the formulas, I just write the equation for recursive call Euclid(b, aq*b, recx, recy) on paper, group everything around a and b and get new solution. The swapping of x and y is not necessary. It is just a hack which make code a bit shorter=) So as soon as you understand this algorithm it is relatively easy to code.
About checking Diophantine code for all sorts of weird cases...=) I made a training with such a problem. The absolute values of a, b, d are up to 10^9 there. You have to output any signed 32bit integer solution. If you want to check yourself, I can send you tests + checker over email (stgatilov@gmail.com)


syg96 was very kind and provided me with his code, as well as, some input/output data. Big thanks for that! I used this to find some small bugs in my code: overflowing int, not handling Math.abs(c)%gcd correctly when gcd==0. I fixed these bugs and now the code passes all tests and now I am 99% sure that it is correct. BTW I am now handling a==0, b==0 and c==0 as a special case. 
