A Modular Euclidean Algorithm

In this handout we cover a modular version of the Euclidean algorithm. This provides a way to control the coefficient growth of the Euclidean algorithm of polynomials over coefficient fields like $\renewcommand{\Q}{\mathbb{Q}}\Q$. Note that applying the usual Euclidean algorithm on polynomials with coefficients in $\Q$ typically causes a great increase in the size of the numerators and denominators of the intermediate coefficients used in the algorithm (and in the coefficients of the $s,t\in\Q[x]$ provided by the extended Euclidean algorithm, which typically explode in size even when run on coprime $a,b\in\Q[x]$ with small integer coefficients).

In [1]:
# An example demonstrating the coefficient growth that occurs in the Euclidean algorithm in Q[x]
F.<x> = QQ[]
a = F(random_vector(ZZ, 10, 10).list())
b = F(random_vector(ZZ, 10, 10).list())
g, s, t = xgcd(a,b)
print(a)
print(b)
print(s)
3*x^9 + 7*x^8 + 4*x^7 + 8*x^6 + 2*x^5 + 3*x^4 + 4*x^3 + 6*x^2 + 9*x
5*x^9 + 7*x^8 + 6*x^7 + 3*x^6 + 4*x^5 + 7*x^4 + x^3 + 4*x^2
6443632968160/118131505340139*x^7 - 34866263779/6217447649481*x^6 - 958350531281/39377168446713*x^5 - 214215554689/39377168446713*x^4 + 5207481762998/118131505340139*x^3 + 794474335838/118131505340139*x^2 - 5285189195002/118131505340139*x + 1/9

Additionally, the modular approach also works in $\renewcommand{\Z}{\mathbb{Z}}\Z[x]$ (not just $\Q[x]$).

GCDs in $\Z[x]$

We've seen that the Euclidean algorithm does not work in $\Z[x]$ since $\Z$ is not a field. A priori it is not even clear if the concept of GCD makes sense in $\Z[x]$ as not every ring has unique factorization. An example of a ring that does not have unique factorization (and therefore does not have GCDs) is the polynomial ring $\Z[x]$ with arithmetic performed modulo $x^2-3$ (typically denoted $\Z[x]/\langle x^2-3\rangle=\{a+b\sqrt3:a,b\in\Z\}$).

Disregarding this, a theorem of Gauss implies that GCDs do in fact exist in $\Z[x]$ and we will develop an algorithm to compute them.

Irreducible polynomials

A polynomial $f$ in $\Z[x]$ is called irreducible if it cannot be factored any further in $\Z[x]$, i.e., the decomposition $f=gh$ must be trivial (one of $g,h\in\Z[x]$ is invertible and thus $\pm1$).

For example, $x^2-1$ is not irreducible, since it factors as $(x-1)(x+1)$.

Note that the irreducibility of a polynomial can depend on its coefficient ring. For example, $x^2-2$ is irreducible over $\Z$ but not over $\renewcommand{\R}{\mathbb{R}}\R$. Conversely, $2x+2$ is not irreducible over $\Z$ as it factors as $2\cdot(x+1)$ which is nontrivial in $\Z$ (neither factor is invertible). However, $2x+2$ is irreducible over $\R$, since the factorization $2(x+1)$ is trivial over $\R$ as $2$ is invertible in $\R$.

Gauss' lemma

A polynomial is called primitive if the greatest common divisor of its coefficients is $1$.

For example, $6x^2+2x+3$ is primitive as $\gcd(6,2,3)=1$ but $6x+3$ is not primitive as $\gcd(6,3)=3$.

A property of integer polynomials proven by Gauss is that the product of two primitive polynomials is also a primitive polynomial.

Furthermore, a nonconstant polynomial $f$ is irreducible (over $\Z$) if and only if $f$ is primitive and $f$ is irreducible (over $\Q$).

In other words, for nonconstant primitive polynomials irreducibility over $\Z$ and irreducibility over $\Q$ correspond exactly.

These properties are known as Gauss' lemmas and using them it follows that $\Z[x]$ has unique factorization because $\Q[x]$ has unique factorization. More generally, if $R$ has unique factorization then $R[x]$ also has unique factorization.

Simplifying assumption

Say $f,g\in\Z[x]$ and we want to compute $\gcd(f,g)$ over $\Z$. It is not a restrictive assumption to assume that $f$ and $g$ are primitive, because if they were not it is easy to compute their "primitive parts" by dividing through by the greatest common divisor of their coefficients first.

Let $\renewcommand{\pp}{\operatorname{pp}}\pp(f)$ be defined to be $f/\gcd(f_0,f_1,\dotsc,f_n)$. In order to compute the GCD of $f$ and $g$ it sufficies to compute the GCD of the "non-primitive" parts (i.e., $\gcd(f_0,f_1,\dotsc,f_n,g_0,g_1,\dotsc,g_m)$) and the GCD of the primitive parts $\pp(f)$ and $\pp(g)$. Thus, from now on we will assume that $f$ and $g$ are primitive. By Gauss' lemma this also implies their product is primitive and $\pp(fg)=\pp(f)\cdot\pp(g)$.

Computing GCDs in $\Z[x]$ via GCDs in $\Q[x]$

As stated above, we assume that $f,g\in\Z[x]$ are primitive and we want to compute their GCD over $\Z$. We already know how to compute their GCD over $\Q$ using the Euclidean algorithm.

Let $v:=\gcd_{\Q[x]}(f,g)$ be the result of applying Euclid's algorithm. As we previously saw, by construction $v$ will be monic, i.e., have a leading coefficient of $1$. However, its other coefficients will very likely be over $\Q$ and not over $\Z$; thus it is not acceptable as a GCD over $\Z$.

Corollary 6.10 in Modern Computer Algebra states that if $h$ is the GCD of $f$ and $g$ over $\Z$ then $h$ is primitive and

\begin{equation*}\renewcommand{\lc}{\operatorname{lc}} h / \lc(h) = v \qquad\text{where $\lc(h)$ is the leading coefficient of $h$} . \end{equation*}

Thus, we need to multiply $v$ by $\lc(h)$ in order to compute $h$. Of course, we don't know $\lc(h)$ since we don't know $h$. However, we can find a multiple of $\lc(h)$. Because $h$ divides $f$ and $g$ (by definition it is the largest divisor) it also follows that $\lc(h)$ divides $\lc(f)=f_n$ and $\lc(g)=g_m$ and thus also $\gcd(f_n,g_m)$.

It follows $\gcd(f_n,g_m)\cdot v$ is an integer polynomial which is a constant multiple of $h$. It may be a nontrivial multiple (introducing a nonprimitive part) but in such a case we can just take its primitive part as $h$ must be primitive.

In summary, when $f$ and $g$ are primitive we have

\begin{equation*} \gcd\nolimits_{\Z[x]}(f,g) = \pp\bigl(\gcd(f_n,g_m)\cdot\gcd\nolimits_{\Q[x]}(f,g)\bigr) . \end{equation*}

Example

Suppose $\tilde f:=30x^3-10x^2+30x-10$ and $\tilde g:=6x^2-14x+4$.

Since $\gcd(30,-10,30,-10)=10$ and $\gcd(6,-14,4)=2$ we can divide $\tilde f$ by $10$ and $\tilde g$ by $2$ to obtain their primitive parts and take $\gcd_{\Z[x]}(\tilde f,\tilde g)=2\gcd_{\Z[x]}(\tilde f/10,\tilde g/2)$.

Now suppose $f:=\tilde f/10 = 3x^3-x^2+3x-1$ and $g:=\tilde g/2=3x^2-7x+2$. We can compute $\gcd(f,g)$ over $\Q$ as $x-1/3$:

In [2]:
R.<x> = QQ[]
f = 3*x^3-x^2+3*x-1
g = 3*x^2-7*x+2
gcd(f,g)
Out[2]:
x - 1/3

Furthermore, the leading coefficients of $f$ and $g$ is $f_3=g_2=3$, so $\gcd(f_3,g_2)=3$.

It follows that $\gcd_{\Z[x]}(f,g)=\pp(3\cdot(x-1/3))=\pp(3x-1)=3x-1$ and $\gcd_{\Z[x]}(\tilde f,\tilde g)=2(3x-1)=6x-2$.

In [3]:
R.<x> = ZZ[]
ftilde = 30*x^3-10*x^2+30*x-10
gtilde = 6*x^2-14*x+4
gcd(ftilde,gtilde)
Out[3]:
6*x - 2

Reducing modulo $p$

The idea behind the modular GCD algorithm is that will reduce the coefficients of $f$ and $g$ modulo a prime $p$, perform Euclid's algorithm on $\renewcommand{\F}{\mathbb{F}}f,g$ (as elements of $\F_p[x]$), and recover $h:=\gcd_{\Z[x]}(f,g)$ from $\gcd_{\F_p[x]}(f,g)$. In order for the recovery to work correctly $p$ must be large enough so that all of the true (non-reduced) coefficients of $h$ lie in the range $\bigl\{-\frac{p-1}{2},\dotsc,\frac{p-1}{2}\bigr\}$. This is the "symmetric" representation of $\F_p$ and it is used instead of the standard representation (that is, $\{0,\dotsc,p-1\}$) because $h$ may have negative coefficients.

However, some primes $p$ cause problems with this approach. For example, consider $p=3,5,7$ and computing $\gcd_{\F_p[x]}(f,g)$ for the above primitive polynomials $f:=3x^3-x^2+3x-1=(x^2+1)(3x-1)$ and $g:=3x^2-7x+2=(x-2)(3x-1)$.

In [4]:
for p in [3, 5, 7]:
    F.<x> = GF(p)[]
    f = 3*x^3-x^2+3*x-1
    g = 3*x^2-7*x+2
    print("gcd_F{}[x](f, g) = {}".format(p, gcd(f,g)))
gcd_F3[x](f, g) = 1
gcd_F5[x](f, g) = x^2 + x + 4
gcd_F7[x](f, g) = x + 2

We have the following:

\begin{align*} \gcd\nolimits_{\F_3[x]}(f,g) &= 1 \\ \gcd\nolimits_{\F_5[x]}(f,g) &= x^2 + x - 1 \\ \gcd\nolimits_{\F_7[x]}(f,g) &= x + 2 \end{align*}

Note that in the last case ($p=7$) the algorithm works correctly: $\gcd(\lc(f),\lc(g))\cdot\gcd\nolimits_{\F_7[x]}(f,g)\equiv3(x+2)\equiv3x-1\pmod{7}$ is the true GCD of $f$ and $g$ over $\Z$.

However, in the first two cases ($p=3,5$) the algorithm does not work correctly, as the degree of $\gcd_{\F_p[x]}(f,g)$ is not correct (too small when $p=3$ and too large when $p=5$). What is going on here?

A criteria for nontrivial GCDs

Suppose $F$ is a field and $f,g\in F[x]$ and $\gcd(f,g)=h$ over $F$. Recall that Euclid's algorithm allows us to find $s,t\in F[x]$ with $sf+tg=h$.

If $h\neq 1$ then there is a nontrivial solution to the equation

\begin{equation*} \text{$sf+tg=0$ with $\deg(s)<\deg(g)$ and $\deg(t)<\deg(f)$.}\tag{$*$} \end{equation*}

Namely, one can take $s:=g/h$ and $t:=-f/h$. In fact, the existence of such a $(s,t)$ provide a certificate that $\gcd(f,g)$ is nontrivial (see lemma 6.13 in Modern Computer Algebra).

Thus, equation ($*$) can be used to determine if $\gcd(f,g)$ is trivial or nontrivial; if ($*$) has a solution then $\gcd(f,g)\neq1$ and if ($*$) has no solution then $\gcd(f,g)=1$.

Note that ($*$) can equivalently be written as the following matrix-vector product equation:

\begin{equation*} \begin{bmatrix} f_n & & & & g_m \\ f_{n-1} & f_n & & & \vdots & g_m \\ \vdots & f_{n-1} & \ddots & & g_0 & \vdots & g_m \\ f_1 & \vdots & & f_n & & g_0 & \vdots & g_m \\ f_0 & f_1 & & f_{n-1} & & & g_0 & \vdots & \ddots \\ & f_0 & & \vdots & & & & g_0 & & g_m \\ & & \ddots & f_1 & & & & & \ddots & \vdots \\ & & & f_0 & & & & & & g_0 \end{bmatrix} \begin{bmatrix} s_{m-1} \\ s_{m-2} \\ \vdots \\ s_0 \\ t_{n-1} \\ t_{n-2} \\ \vdots \\ t_0 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix} \in F^{n+m} \end{equation*}

Note $\deg(f)=n$ and $\deg(g)=m$ and the $i$th row of this expression corresponds to the coefficient of the $x^{n+m-i}$ and hence why the right-hand side contains all zeros, as there are no terms $x^{n+m-i}$ on the right-hand side of ($*$). The matrix in this expression is known as the Sylvester matrix of $f$ and $g$.

Linear algebra then tells us that $\gcd(f,g)\neq1$ if and only if the Sylvester matrix of $f$ and $g$ is singular (i.e., there is a nontrivial solution of this matrix-vector equation).

The determinant of the Sylvester matrix of $f$ and $g$ is known as the resultant $\renewcommand{\res}{\operatorname{res}}\res(f,g)$. A matrix is singular if and only if its determinant is $0$, so we can equivalently state this as

\begin{equation*} \gcd(f,g)\neq1 \Longleftrightarrow \res(f,g)=0. \end{equation*}

The above takes place over a field $F$ but due to Gauss' theorem it can also be modified to work over $\Z$:

\begin{equation*} \gcd\nolimits_{\Z[x]}(f,g)\text{ is nonconstant} \Longleftrightarrow \res(f,g)=0. \end{equation*}

A criteria for choosing primes

The reason we introduced the resultant is because it allows an easy specification of the primes $p$ for which the GCD over $\F_p$ can be used to find the GCD over $\Z$.

Theorem (6.26, Modern Computer Algebra)

Let $f,g\in\Z[x]$ be nonzero and of degrees $n$ and $m$, let $h=\gcd(f,g)$ over $\Z$, and let $p$ be a prime that does not divide $\gcd(f_n,g_m)$.

Then $\deg\gcd_{\F_p[x]}(f,g)\geq\deg h$ (i.e., polynomials might split "deeper" modulo $p$, causing $\gcd(f,g)$ over $\F_p$ to be larger than the $\gcd(f,g)$ over $\Z$.)

Moreover, the degree of $\gcd(f,g)$ over $\F_p$ will be equal to the degree of $\gcd(f,g)$ over $\Z$ if and only if $p$ does not divide $\res(f/h,g/h)$.

Furthermore, $p$ does not divide $\res(f/h,g/h)$ exactly when $\gcd\nolimits_{\F_p[x]}(f,g) \equiv h/\lc(h) \pmod{p}$. (The inverse of $\lc(h)\bmod p$ exists since $\lc(h)$ divides $\gcd(f_n,g_m)$ which does not have $p$ as a divisor.)

The modular algorithm for GCDs

So the prime $p$ must satisfy the following:

  1. $p$ does not divide $\gcd(f_n,g_m)$
  2. $p$ does not divide $\res(f/h,g/h)$
  3. The coefficients of $\gcd(f_n,g_m)\cdot h/\lc(h)$ have absolute value at most $(p-1)/2$ so they fit in the symmetric range

If $p$ satisfies all three conditions then we can compute $h$, the $\gcd(f,g)$ over $\Z$, by:

  • Using the Euclidean algorithm to compute $\gcd(f,g)$ over $\F_p$
  • Multiplying this computed GCD by $\gcd(f_n,g_m)$ and reduce the coefficients to be in the symmetric range modulo $p$
  • Return the primitive part of the above polynomial

Example

Let's compute the integer GCD of $f:=3x^3-x^2+3x-1=(x^2+1)(3x-1)$ and $g:=3x^2-7x+2=(x-2)(3x-1)$ using this approach.

First, $p$ must not divide $\gcd(f_3,g_2)=\gcd(3,3)=3$. Thus $p\neq3$

Recall that $f/h=x^2+1$ and $g/h=x-2$, and the Sylvester matrix of these two polynomials is

\begin{equation*} \begin{bmatrix} 1 & 1 & 0 \\ 0 & -2 & 1 \\ 1 & 0 & -2 \end{bmatrix} \end{equation*}

which has determinant $(-2)^2+1=5$. Thus $p\neq5$.

The coefficients of $\gcd(f_3,g_2)\cdot(3x-1)/3=3x-1$ have absolute value at most 3, so we must have $(p-1)/2\geq3$, i.e., $p\geq7$.

Thus, the simplest selection is $p=7$.

As we saw above, Euclid's algorithm computes $\gcd\nolimits_{\F_7[x]}(f,g) = x + 2$. We multiply this by $\gcd(f_3,g_2)=3$ to obtain $3x+6$ which when reduced to the symmetric range is $3x-1$ which is already primitive.

Caveats

One unrealistic part of this example: the conditions on $p$ involved $h$ so in order to properly select $p$ we are required to know $h=\gcd(f,g)$ over $\Z$. But that's the very thing we are trying to compute!

How can we get around this?

We could derive an upper bound on $\res(f/h,g/h)$ and then select $p$ to be larger than this. However, this is very wasteful in practice and tends to use a prime $p$ much larger than necessary. So we will ignore the resultant condition for now.

What about the sizes of the coefficients of $h$? It can be shown that the maximum coefficient of $h$ has absolute value at most $\sqrt{n+1}\cdot2^nA$ where $A$ is an upper bound on the coefficients of $f$ and $g$. Thus if we choose a prime larger than $B:=2\gcd(f_n,g_m)\sqrt{n+1}\cdot2^nA$ then we can guarantee that all coefficients of $h$ will be bounded in absolute value by $(p-1)/2$.

It can also be shown that if you choose a random prime between $B$ and $2B$ then $p$ will not divide $\gcd(f/h,g/h)$ with probability at least $1/2$. In other words, it shouldn't be hard to find a prime that works.

How can you tell if a prime works? The simplest approach is simply to verify that the purported GCD is actually a divisor of both $f$ and $g$. If so, it follows that $p$ does not divide $\res(f/h,g/h)$. Why? Because if $p$ did divide $\res(f/h,g/h)$ then by Thm 6.26 the degree of $\gcd_{\F_p[x]}(f,g)$ will be strictly larger than the true GCD $h$. In such a case $\gcd_{\F_p[x]}(f,g)$ cannot possibly divide both $f$ and $g$ (over $\Z$) because then it would also have to divide their GCD $h$ which is nonsensical given that $\gcd_{\F_p[x]}(f,g)$ has a larger degree than $h$.