Tag / extended

by Carsten König

extended Euclidean algorithm in F#

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) d of two (normally natural) numbers a and b. As the name suggests this is the greatest of all divisors of both a and b. There is always such a number because obviously 1 divides every number and no number bigger than the minimum of a and b can be a common divisor as such a number won’t divide min\left\{a,b\right\}. That leaves you with an nonempty, bounded set of natural numbers that must have a maximum – namely d.

One way to find it is based on the observation that

gcd(a, b) = gcd (b, r)

where r is the remainder of dividing a by b. Without going into every messy detail: You can see this by writing a = q b+r where 0 \leq r < q (division with remainder) and noting that d will divide both a and b IIF d will divide a-q b and r.

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

gcd(a,0) = a

and if we now iterate the step till one remainder gets 0 the gcd 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 a = r_{-2}, b = r_{-1} and the step as gcd(r_{i-2}, r_{i-1}) = gcd (r_{i-1}, r_{i}) with r_{i-2} = q_{i}r_{i-1} + r_{i}. As F# might have some issues with subscript i will use r’’ for r_{i-2}, r’ for r_{i-1} 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 gcd but two additional integers s and t so that

gcd(a,b) = s a + t b

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 Zwinkerndes Smiley – 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

r_{i-2} = s_{i-2}a + t_{i-2}b

and

r_{i-1} = s_{i-1}a + t_{i-1}b

then it follows from r_{i-2}=q_{i}r_{i-1}+r_{i} that

r_{i}=a (s_{i-2}-q_{i}s_{i-1})+b(t_{i-2}-q_{i}t_{i-1})

and we see that s_{i} = s_{i-2}-q_{i} s_{i-1} and t_{i} = t_{i-2}-q_{i} t_{i-1} which gives us a way to recursively define the s_i and t_i values.

As a last step we only have to notice that the first of those values are just s_{-2} = 1, s_{-1} = 0, t_{-2} = 0 and t_{-1} = 1 because – of course – a = 1a + 0b and b = 0a + 1b. 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 gcd(a,b) = d = s a + t b.

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 gcd(76, 16) = 4 = -1*76+5*16