# 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 – 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$