The Euclidean algorithm

The Euclidean algorithm is perhaps the oldest algorithm in the world, being attributed to Euclid and appearing in his Elements. It applies on any Euclidean domain:

Integers and polynomials form a Euclidean domain, where the metric for integers is simply the absolute value, and for polynomials it is the degree of the polynomial.

Before we get to Euclid's algorithm, we need one more fact about GCDs: For any integer $x$, $(a,b) = (b,a) = (a,-b) = (a,b+ax)$.

Euclid's algorithm finds the greatest common divisor of a pair elements in a Euclidean domain.

The Euclidean algorithm works by repeated division.



This can be coded in nothing at all in MATLAB:

We can do more, however. It is a fact of number theory that

$$(a,b) = ax + by$$

for some integers $x$ and $y$. We can use the Euclidean algorithm to determine $x$ and $y$. In many of our applications, these will be the numbers that we need. (Another useful fact: if there is an $x$ and $y$ such that $ax + by=1$, then $(a,b)=1$.)

More formally, the algorithm may be stated in the following theorem, which also fixes some useful notation.

Theorem 1 (The Euclidean algorithm)   Let $b$ and $c$ be integers $>0$. Then by repeated application of the division algorithm write

$$\begin{aligned} b &= cq_1 + r_1 \qquad 0 < r_1 < c \\ c &= r_1 q... ..._j + r_j \qquad 0 < r_j < r_{j-1} \\ r_{j-1} &= r_j q_{j+1} + 0. \end{aligned}$$

Then $(b,c) = r_j$, the last nonzero remainder of the division process.

That the theorem stops after a finite number of steps follows since every remainder must be less than the preceding remainder.

If the GCD $g$ and the coefficients $x$ and $y$ are desired such that

$$b x + cy = g,$$

a little more work is required. Observe that the recursion in the Euclidean algorithm may be expressed as

$$q_i = \lfloor r_{i-2}/r_{i-1} \rfloor$$

$$r_{i} = r_{i-2} - r_{i-1}q_i$$ (1)

for $i=1,2,\ldots$ (until termination) with $r_{-1} = b$ and $r_0 = c$. The values for $x$ and $y$ may be obtained by finding intermediate integers $x_i$ and $y_i$ satisfying

$$b x_i + c y_i = r_i$$

In conjunction with (1) we obtain

$$\begin{split}x_i &= x_{i-2} - q_i x_{i-1} \ y_i &= y_{i-2} - q_i y_{i-1}\end{split}$$ (2)

for $i=1,2,\ldots$ (until termination) with

$$\begin{aligned} x_{-1} &= 1\qquad x_0 = 0 \\ y_{-1} &= 0\qquad y_0 = 1 \end{aligned}$$

This algorithm applies as well to polynomials (over a field) as to integers. The following implementation is for polynomials.

As an example, we can try

>> [g,x,y] = gcdpoly([4 10 8 2],[8 14 7 1])
g =
    1.0000    1.5000    0.5000
x =
    0.3333
y =
   -0.1667