The Extended Extended Euclidean Algorithm

Recently I've been working on a Rust library implementing some structures and algorithms from commutative algebra. One of the algorithms I implemented is the Euclidean algorithm which computes a greatest common divisor $g$ of two elements $a, b$ in a Euclidean domain. There is also the well-known extended Euclidean algorithm which computes not only the gcd but also Bezout coefficients $x, y$ satisfying $g = xa + yb$. What I hadn't realized until recently (and I never saw it mentioned in the literature) is that the Euclidean algorithm can be further extended to also compute the quotients $d$ and $e$ satisfying $a = d g$ and $b = e g$.

The Classical Euclidean Algorithm

Let $(A, \delta, {/})$ be a Euclidean domain, i.e. an integral domain $A$ together with a function $\delta \colon A\setminus\{0\} \to \IN$ and an operation ${/} \colon A \times A \setminus \{0\} \to A$ such that for every $a, b \in A$ with $b \neq 0$ we have either $b \mid a$ or $\delta(a - b \cdot (a / b)) < \delta(b)$.

The Euclidean algorithm is based on the observation that for $a, b \in A$, $b \neq 0$ the two ideals $\langle a, b\rangle$ and $\langle b, a - b \cdot (a / b)\rangle$ are equal (it is easy to see that the two generators of either ideal are contained in the other ideal). If we define inductively the pairs $(a_i, b_i)$ by $(a_0, b_0) := (a, b)$ and $(a_{i + 1}, b_{i + 1}) := (b_i, a_i - b_i (a_i / b_i))$ whenever $b_i \neq 0$, then we necessarily reach an $N \in \IN$ where $b_N = 0$, because we have $\delta(b_{i + 1}) < \delta(b_i)$ as long as $b_{i + 1} \neq 0$.

Consequently we have $\langle a, b\rangle = \langle a_0, b_0\rangle = \dots = \langle a_N, b_N\rangle = \langle a_N\rangle$. It is easy to see that this implies that an element of $A$ divides $a_N$ if and only if it divides both $a$ and $b$, thereby establishing that $g := a_N$ is a gcd of $a$ and $b$.

The Extended Euclidean Algorithm

The equality of ideals $\langle a, b\rangle = \langle g\rangle$ tells us already that there exist $x, y \in A$ satisfying $xa + yb = g$ (as well as $d, e \in A$ satisfying $a = dg$ and $b = eg$ for that matter). However we need to track down the inducitve construction of $g$ in order to get a computation of $x$ and $y$ as well. It is useful to phrase this in the language of matrices.

We can rephrase the construction of the sequence $(a_i, b_i)$ as follows:

We have $$ \begin{pmatrix}a_0 \\ b_0\end{pmatrix} = \begin{pmatrix}a \\ b\end{pmatrix}, \qquad \begin{pmatrix}a_{i + 1} \\ b_{i + 1}\end{pmatrix} = \underbrace{\begin{pmatrix}0 & 1 \\ 1 & - a_i / b_i\end{pmatrix}}_{M_i}\cdot \begin{pmatrix}a_i \\ b_i\end{pmatrix}. $$ We thus have $$ \begin{pmatrix}a_i \\ b_i\end{pmatrix} = N_i \cdot \begin{pmatrix}a \\ b\end{pmatrix} $$ where $N_i = M_i \cdot \dots \cdot M_1$ (in the case $i = 0$ this is the $2 \times 2$ unit matrix). We want to find inductive formulas for the coefficients $s_i, t_i, u_i, v_i$ satisfying $$ N_i = \begin{pmatrix}s_i & t_i \\ u_i & v_i\end{pmatrix}. $$ We have $$ \begin{pmatrix}s_0 & t_0 \\ u_0 & v_0\end{pmatrix} = N_0 = \begin{pmatrix}1 & 0 \\ 0 & 1\end{pmatrix} $$ and

$$ \begin{split} \begin{pmatrix}s_{i + 1} & t_{i + 1} \\ u_{i + 1} & v_{i + 1}\end{pmatrix} &= M_{i + 1} \cdot N_i \\ &= \begin{pmatrix}0 & 1 \\ 1 & -a_i / b_i\end{pmatrix} \cdot \begin{pmatrix}s_i & t_i \\ u_i & v_i\end{pmatrix} \\ &= \begin{pmatrix} u_i & v_i \\ s_i - u_i \cdot (a_i / b_i) & t_i - v_i \cdot (a_i / b_i) \end{pmatrix}. \end{split} $$

In particular, we have $$ \begin{pmatrix}a_N \\ b_N\end{pmatrix} = N_i \cdot \begin{pmatrix}a_0 \\ b_0\end{pmatrix} = \begin{pmatrix}s_N & t_N \\ u_N & v_N\end{pmatrix} \cdot \begin{pmatrix}a \\ b\end{pmatrix}, $$ i.e. $g = a_N = s_N a + t_N b$.

Thus, setting $x := s_N$ and $y := t_N$ gives us expressions for the Bezout coefficients which can be inductively computed.

The Extended Extended Euclidean Algorithm

The matrices $$M_i = \begin{pmatrix}0 & 1 \\ 1 & -a_i / b_i\end{pmatrix}$$ are obviously invertible with determinant equal to $-1$. Hence, $N_N := M_1 \cdots M_N$ is invertible with determinant $(-1)^N$. Computing the inverse in terms of the determinant and the adjugate matrix, we find that $$ N_N^{-1} = (-1)^N \begin{pmatrix}v_N & -t_N \\ -u_N & s_N \end{pmatrix}. $$ From $$ \begin{pmatrix}g \\ 0\end{pmatrix} = \begin{pmatrix}a_N \\ b_N\end{pmatrix} = N_N \cdot \begin{pmatrix}a \\ b\end{pmatrix} $$ we get $$ \begin{pmatrix}a \\ b\end{pmatrix} = N_N^{-1} \begin{pmatrix}g \\ 0\end{pmatrix} = (-1)^N \begin{pmatrix}v_N & -t_N \\ -u_N & s_N \end{pmatrix} \begin{pmatrix}g \\ 0\end{pmatrix}. $$ In particular, setting $d := (-1)^N v_N$ and $e := (-1)^{N + 1} u_N$ we have $a = dg$ and $b = eg$.

Is It Worth It?

Well, probably not really. In most standard examples of Euclidean domains, like the integers $\mathbb{Z}$ or the polynomial ring over a field $K[T]$, it is easy to compute the quotients $d = \frac{a}{g}$, $e = \frac{b}{g}$ (when $g \neq 0$). However the computation of these quotients is not part of the structure of a Euclidean domain, so for my Rust implementation which is supposed to operate on arbitrary Euclidean domains, I didn't have it available which was my motivation to figure out how to read the off from the Euclidean algorithm.

In fact it might be the case that the computation of quotients where the numerator is a multiple of the denominator actually is part of the structure of a Euclidean domain: At least in the examples of $\mathbb{Z}$ and $K[T]$ it coincides with the Euclidean division function. But I am not sure if this is the case for all Euclidean domains. It certainly is the case when you require some more properties for the Euclidean norm function, like multiplicativity and positivity away from $0$ (as in the case of $\mathbb{Z}$).

So I want to end this post with a question:

Question. Let $(A, \delta, /)$ be a Euclidean domain and let $a, b \in A$ where $b \neq 0$ and $b \mid a$. Is it necessarily the case that $a = b \cdot (a / b)$? If not, which additional properties on $(A, \delta, /)$ are necessary to make this conclusion?