The Extended Euclidean Algorithm
--------------------------------

There are two prevailing versions of expositions of the extended
Euclidean algorithm.  All students I have talked to find both
confusing.  Certainly I felt that there ought to be a better way
after I managed to understand them.  Here is my take.

Given two numbers a and b, the algorithm should end with this equation:

  a*m + b*n = d

where d = gcd(a,b).

The algorithm accomplishes this by working with a pair of equations
instead.  At the beginning, they are:

  a*1 + b*0 = a
  a*0 + b*1 = b

Then we go through iterations of modifying these two equations.
In each iteration, the pair is modified to keep this form:

  a*u + b*v = w    (1)
  a*x + b*y = z    (2)

with w and z satisfying:

  gcd(w,z) = gcd(a,b)    (*)

(Note that the initial pair satisfies (*) already.)

In other words, u,v,w,x,y,z are initially 1,0,a,0,1,b.  Then we keep
updating the equations, changing the values of u,v,w,x,y,z
accordingly, but we have to ensure that gcd(w,z) always happens to be
gcd(a,b).  (This despite the fact that we do not yet know the value of
gcd(a,b) until the end.)

How to do this?  This is where we extend the Euclidean algorithm.
Initially, w,z is a,b, so gcd(w,z) = gcd(a,b) at the beginning for
free.  Now in each iteration, if gcd(w,z) is gcd(a,b), then certainly
gcd(w mod z, z) is still gcd(a,b).  Thus we can indeed maintain the
condition (*) without even knowing the value of gcd(a,b).

So we replace w by w mod z.  Assuming w>=z, this is not a waste of
time.  And u,v,x,y need to be updated in sync.  The new pair of
equations is:

  a*(u-qx) + b*(v-qy) = w-qz
  a*x + b*y = z

where q = w div z (and so w-qz = w mod z, as promised).
I.e., compute q, then replace (1) by (1)-q*(2).

(In the case that w<z, you can do the symmetric thing, or just
swap the two equations and do the above.)

What's good of doing this?  Because this process is just the plain old
Euclidean algorithm when you just look at the values of w,z,
eventually one of them, say w, will become 0, and z will become
gcd(a,b):

  a*u + b*v = 0
  a*x + b*y = gcd(a,b)

Then you take the second equation and you are done!