I wanted to write something about this for some time. Yes there are very good articles out there – for example in Wikipedia (see here and here). But why not another with some explanations and code in F#?

The basic algorithm

As you can read in those superb articles I linked above, the Euclidean algorithm is used to find the GCD (greatest common divisor) of two (normally natural) numbers and . As the name suggests this is the greatest of all divisors of both and . There is always such a number because obviously divides every number and no number bigger than the minimum of and can be a common divisor as such a number won’t divide . That leaves you with an nonempty, bounded set of natural numbers that must have a maximum – namely .

One way to find it is based on the observation that

where is the remainder of dividing by . Without going into every messy detail: You can see this by writing where (division with remainder) and noting that will divide both and IIF will divide and .

Based on this we can already write down the Euclidean algorithm: Because every number divides 0 we see that

and if we now iterate the step till one remainder gets 0 the will be just the remainder we calculated before:

let gcd (a : bigint) (b : bigint) = let rec inner r'' r' = let step () = let q = r'' / r' let r = r'' - q*r' r if r' = 0I then r'' else inner r' (step()) inner a b

here I just wrote the step based on the latest two remainders (r’’ and r’) – generalizing what was said above by writing , and the step as with . As F# might have some issues with subscript will use r’’ for , r’ for and similar notation for other values than r.

The extended algorithm

The extended algorithm takes this idea and extends it to find not only the but two additional integers and so that

Normally you will write down the series of steps the Euclidean algorithm will be taking and then move from bottom to top always replacing the the numbers by what you found (see the wiki articles) – this is not only boring but rather error prone and of course should be taken care by the computer – I did not go into details because in my opinion it’s easier to see the algorithm directly by taking the general step.

So let’s say we already know that

and

then it follows from that

and we see that and which gives us a way to recursively define the and values.

As a last step we only have to notice that the first of those values are just , , and because – of course – and . Finally we get our algorithm in the form:

let extGCD (a : bigint) (b : bigint) = let rec inner (r'', s'', t'') (r', s', t') = let step () = let q = r'' / r' let r = r'' - q*r' let s = s'' - q*s' let t = t'' - q*t' (r, s, t) if r' = 0I then (r'', s'', t'') else inner (r', s', t') (step()) inner (a, 1I, 0I) (b, 0I, 1I)

That takes two bigint and results in a triple (d, s, t) with .

Here is a quick example for this:

> extGCD 76I 16I;; val it : System.Numerics.BigInteger * System.Numerics.BigInteger * System.Numerics.BigInteger = (4 {IsEven = true; IsOne = false; IsPowerOfTwo = true; IsZero = false; Sign = 1;}, -1 {IsEven = false; IsOne = false; IsPowerOfTwo = false; IsZero = false; Sign = -1;}, 5 {IsEven = false; IsOne = false; IsPowerOfTwo = false; IsZero = false; Sign = 1;}) >

stating that