Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

I just asked this question on stackoverflow.com and had it closed before I could get any reasonable help, with the suggestion to move it to a math site. I don't understand the math, and don't speak the math speak, and need an algorithm, not simply math jargon for the answer.

With a, b, and c as known values, % being the typical computer modulus operator, and with a and b relatively prime, how do I compute a legitimate x? I know that x exists, that's how c came into being... I just need to solve for x.

Someone on stackoverflow pointed me to wikipedia on Linear Diophantine equations for a solution to ax + by = c where x and y are integers (which is almost a re-write of ax % b = c, I need ax - by = c), but I don't know what half this stuff means, or how to implement it in code.

"Let g be the greatest common divisor of a and b. Both terms in ax + by are divisible by g; therefore, c must also be divisible by g, or the equation has no solutions. By dividing both sides by c/g, the equation can be reduced to Bezout's identity sa + tb = g where s and t can be found by the extended Euclidean algorithm."

Ok, that's great in theory... but. How do I use such a thing to find s and t? And even further, how would I use s to find x?

Thanks.

share|improve this question
    
Google "extended Euclidean algorithm", e.g. Wikipedia –  Bill Dubuque Apr 16 '11 at 20:28
    
I have googled this problem over and over. I DO NOT understand the page on extended Euclidean algorithm, or its implications to possibly solving my problem. Once again, I need an algorithm please. –  Steven Apr 16 '11 at 20:38
    
@Steven, the algorithms are on that page. What more could we possibly give you? –  Raskolnikov Apr 16 '11 at 20:47
    
$\rm a\ x\ mod\ b\ =\ c\ $ implies that $\rm\ a\ x + b\ y\ =\ c\ $ for some integer $\rm\:y\:.\:$ Use the extended Euclidean algorithm to solve for $\rm\:x\:$, then return $\rm\ x\ mod\ b\:.$ –  Bill Dubuque Apr 16 '11 at 20:53
    
I do not see how to leverage the Euclidean algorithm to solve ax - by = c. Using the Euclidean Algorithm gives me the gcd(m, n), and the extended version gives me ax + by = 1. I happen to know that gcd(a, b) = 1 (a and b are relatively prime). How do I connect my equation of ax - by = c to the extended algorithm of ax + by = 1? –  Steven Apr 16 '11 at 21:34
show 1 more comment

3 Answers

up vote 4 down vote accepted

Here it is as C-style code.

You can see the Euclidean algorithm at work: p and q get reduced while we simultaneously keep track of both p and q as linear combinations of a and b. In the end this gives us exactly the numbers we need. The comments should allow you to understand it enough to debug it in whatever your favorite language is. The assert statements are basically asserting that c indeed comes from an a*x%b calculation with the given a and b.

I will assume a>=0 and b>0.

if (a == b) {
  assert(c==0);
  return 0; (* any value of x satisfies a*x % b == c in this case *)
}
int pa = 1, pb = 0, p = a, qa = 0, qb = 1, q = b; (* p always == pa*a+pb*b *)
while (p > 0 && q > 0) {
  if (q > p) { (* might be false on first pass, but then always true *)
    m = q / p;  (* integer arithmetic, assuming positive p,q *)
    q -= m * p;
    qa -= m * pa;
    qb -= m * pb;
  }
  if (q > 0) { (* q must be less than p *)
    m = p / q;
    p -= m * q;
    pa -= m * qa;
    pb -= m * qb;
  }
}
(* now one of p or q is zero, and the other is the gcd *)
assert(c % (p+q) == 0); (* c should be a multiple of the gcd, p+q *)
m = c / (p+q);
if (p > 0) (* p==pa*a+pb*b, and m*p==c, so c==m*pa*a+m*pb*b, so x=m*pa *)
  return pa > 0 ? m*pa % b : b - (m*-pa) % b;
else
  return qa > 0 ? m*qa % b : b - (m*-qa) % b;
share|improve this answer
    
This is awesome Matt! You saved my bacon. I have NO idea what is going on here sadly... I'm not sure what we are doing with p and q. I follow the code fine, but don't know why any of it is there. It works perfectly though. Thank you! If someone can explain in "stupidese" what is going on, I will certainly try to follow (and would enjoy the exercise also). –  Steven Apr 16 '11 at 22:01
    
I wish I had enough reputation to upvote this response. This is perfect. I keep looking it over trying to figure out what is going on. I know the gcd for a and b is 1. What I don't know is what you are computing along the way that is making pa the value we need. –  Steven Apr 16 '11 at 22:11
    
Steven: If you desire to understand the algorithm rather than treat the code as a black box then see the example in my post here. Studying obfuscated code in random programming lanuages is a very poor way to understand mathematical concepts. The language of mathematics is much higher-level and much better suited for such purposes. –  Bill Dubuque Apr 16 '11 at 23:46
    
@Steven Even if you can't upvote it, you certainly can check the checkmark next to it to accept it as the answer to your question. –  Matthew Frederick Apr 17 '11 at 3:00
add comment

Here, for example, is the extended Euclidean algorithm in Python. Apply that to your numbers $a$ and $b$, to get what they call $u_1$, $u_2$, $u_3$ corresponding to your $s$, $t$, $g$. Then $x=s \cdot (c/g)$ (where $c/g$ must be an integer since you knew that there was a solution).

share|improve this answer
    
I'm all kinds of confused. I implemented the code in the example you provided, and passed in my a and b into the function, which yielded u1, u2, and u3 - where u3 was 1. I don't know how this can work without c being taken into account somewhere before the final multiplication. The u1 result I got was negative, and u1 * c was an absolutely huge number that when substituted into the original equation did not yield a solution. I wonder if the parameters to that function need to be something other than a and b. I wondered about a/c or b/c, but a/c is not an integer. –  Steven Apr 16 '11 at 21:25
    
@Steven: Perhaps you could tell us your numbers a, b, c? –  Hans Lundmark Apr 16 '11 at 21:58
    
a, b, and c will vary depending on the situation. All I know is that they are known, and that x was also known at one time (used to compute c). Although they are known, the algorithm must work for any values of a, b, and c where the equation holds true. –  Steven Apr 16 '11 at 22:03
    
@Steven: I meant the values that you used for the trials in your previous comment. –  Hans Lundmark Apr 16 '11 at 22:08
    
@Steven: I tried a few numbers here, and it works perfectly, so I wanted to see your specific example. Regarding your comment "I don't know how this can work without c being taken into account somewhere before the final multiplication": Once you know numbers such that $a u_1 + b u_2 = u_3$, you can multiply that whole equation by $k$ to get numbers such that $a (u_1 k) + b (u_2 k) = (u_3 k)$, so if you can find $k$ such that $u_3 k=c$, then $x=u_1 k$ will be the number that you're looking for! –  Hans Lundmark Apr 16 '11 at 22:18
show 2 more comments

Below is an example of the extended Euclidean algorithm from one of my prior posts. It computes the Bezout representation for $\rm\:gcd(80,62) = 2\ $ viz. $\ 7\cdot 80 - 9\cdot 62\ = 2\:.\:$ Now, to solve $\rm\ 80\ x\ mod\ 62\ =\ c\ $ simply multiply the Bezout equation by $\rm\:c/2\:$ if $\rm\:c\:$ is even (if $\rm\:c\:$ is odd there is no solution). For example, for $\rm\:c = 6\:$ we scale Bezout by $\rm\:3\:$ yielding $\ 21\cdot 80 - 27\cdot 62\ =\ 6\ $ therefore $\rm\ 21\cdot 80\ mod\ 62\ =\ 6\:.$ This Bezout scaling method works generally.

For example, to solve  mx + ny = gcd(x,y) one begins with
two rows  [m   1    0], [n   0    1], representing the two
equations  m = 1m + 0n,  n = 0m + 1n. Then one executes
the Euclidean algorithm on the numbers in the first column,
doing the same operations in parallel on the other columns,

Here is an example:  d =  x(80) + y(62)  proceeds as:

                      in equation form   | in row form
                    ---------------------+------------
                    80 =   1(80) + 0(62) | 80   1   0
                    62 =   0(80) + 1(62) | 62   0   1
 row1 -   row2  ->  18 =   1(80) - 1(62) | 18   1  -1
 row2 - 3 row3  ->   8 =  -3(80) + 4(62) |  8  -3   4
 row3 - 2 row4  ->   2 =   7(80) - 9(62) |  2   7  -9
 row4 - 4 row5  ->   0 = -31(80) -40(62) |  0 -31  40

Above the row operations are those resulting from applying
the Euclidean algorithm to the numbers in the first column,

        row1 row2 row3 row4 row5
namely:  80,  62,  18,   8,   2  = Euclidean remainder sequence
               |    |
for example   62-3(18) = 8, the 2nd step in Euclidean algorithm

becomes:   row2 -3 row3 = row4  on the identity-augmented matrix.

In effect we have row-reduced the first two rows to the last two.
The matrix effecting the reduction is in the bottom right corner.
It starts as the identity, and is multiplied by each elementary
row operation matrix, hence it accumulates the product of all
the row operations, namely:

       [  7 -9] [ 80  1  0]  =  [2   7  -9]
       [-31 40] [ 62  0  1]     [0 -31  40]

The 1st row is the particular  solution: 2 =   7(80) -  9(62)
The 2nd row is the homogeneous solution: 0 = -31(80) + 40(62),
so the general solution is any linear combination of the two:

       n row1 + m row2  ->  2n = (7n-31m) 80 + (40m-9n) 62

The same row/column reduction techniques tackle arbitrary
systems of linear Diophantine equations. Such techniques
generalize easily to similar coefficient rings possessing a
Euclidean algorithm, e.g. polynomial rings F[x] over a field, 
Gaussian integers Z[i]. There are many analogous interesting
methods, e.g. search on keywords: Hermite / Smith normal form, 
invariant factors, lattice basis reduction, continued fractions,
Farey fractions / mediants, Stern-Brocot tree / diatomic sequence.
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.