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