{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# A Modular Euclidean Algorithm\n",
"\n",
"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)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"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\n",
"5*x^9 + 7*x^8 + 6*x^7 + 3*x^6 + 4*x^5 + 7*x^4 + x^3 + 4*x^2\n",
"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\n"
]
}
],
"source": [
"# An example demonstrating the coefficient growth that occurs in the Euclidean algorithm in Q[x]\n",
"F. = QQ[]\n",
"a = F(random_vector(ZZ, 10, 10).list())\n",
"b = F(random_vector(ZZ, 10, 10).list())\n",
"g, s, t = xgcd(a,b)\n",
"print(a)\n",
"print(b)\n",
"print(s)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Additionally, the modular approach also works in $\\renewcommand{\\Z}{\\mathbb{Z}}\\Z[x]$ (not just $\\Q[x]$)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## GCDs in $\\Z[x]$\n",
"\n",
"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\\}$).\n",
"\n",
"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.\n",
"\n",
"### Irreducible polynomials\n",
"\n",
"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$).\n",
"\n",
"For example, $x^2-1$ is not irreducible, since it factors as $(x-1)(x+1)$.\n",
"\n",
"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$.\n",
"\n",
"### Gauss' lemma\n",
"\n",
"A polynomial is called _primitive_ if the greatest common divisor of its coefficients is $1$.\n",
"\n",
"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$.\n",
"\n",
"A property of integer polynomials proven by Gauss is that the product of two primitive polynomials is also a primitive polynomial.\n",
"\n",
"Furthermore, a nonconstant polynomial $f$ is irreducible (over $\\Z$) if and only if $f$ is primitive and $f$ is irreducible (over $\\Q$).\n",
"\n",
"In other words, for nonconstant primitive polynomials irreducibility over $\\Z$ and irreducibility over $\\Q$ correspond exactly.\n",
"\n",
"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.\n",
"\n",
"### Simplifying assumption\n",
"\n",
"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.\n",
"\n",
"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)$.\n",
"\n",
"### Computing GCDs in $\\Z[x]$ via GCDs in $\\Q[x]$\n",
"\n",
"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.\n",
"\n",
"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$.\n",
"\n",
"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\n",
"\n",
"\\begin{equation*}\\renewcommand{\\lc}{\\operatorname{lc}}\n",
"h / \\lc(h) = v \\qquad\\text{where $\\lc(h)$ is the leading coefficient of $h$} .\n",
"\\end{equation*}\n",
"\n",
"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)$.\n",
"\n",
"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.\n",
"\n",
"In summary, when $f$ and $g$ are primitive we have\n",
"\n",
"\\begin{equation*}\n",
"\\gcd\\nolimits_{\\Z[x]}(f,g) = \\pp\\bigl(\\gcd(f_n,g_m)\\cdot\\gcd\\nolimits_{\\Q[x]}(f,g)\\bigr) .\n",
"\\end{equation*}\n",
"\n",
"### Example\n",
"\n",
"Suppose $\\tilde f:=30x^3-10x^2+30x-10$ and $\\tilde g:=6x^2-14x+4$.\n",
"\n",
"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)$.\n",
"\n",
"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$:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"x - 1/3"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"R. = QQ[]\n",
"f = 3*x^3-x^2+3*x-1\n",
"g = 3*x^2-7*x+2\n",
"gcd(f,g)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Furthermore, the leading coefficients of $f$ and $g$ is $f_3=g_2=3$, so $\\gcd(f_3,g_2)=3$.\n",
"\n",
"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$."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6*x - 2"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"R. = ZZ[]\n",
"ftilde = 30*x^3-10*x^2+30*x-10\n",
"gtilde = 6*x^2-14*x+4\n",
"gcd(ftilde,gtilde)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reducing modulo $p$\n",
"\n",
"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.\n",
"\n",
"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)$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"gcd_F3[x](f, g) = 1\n",
"gcd_F5[x](f, g) = x^2 + x + 4\n",
"gcd_F7[x](f, g) = x + 2\n"
]
}
],
"source": [
"for p in [3, 5, 7]:\n",
" F. = GF(p)[]\n",
" f = 3*x^3-x^2+3*x-1\n",
" g = 3*x^2-7*x+2\n",
" print(\"gcd_F{}[x](f, g) = {}\".format(p, gcd(f,g)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have the following:\n",
"\\begin{align*}\n",
"\\gcd\\nolimits_{\\F_3[x]}(f,g) &= 1 \\\\\n",
"\\gcd\\nolimits_{\\F_5[x]}(f,g) &= x^2 + x - 1 \\\\\n",
"\\gcd\\nolimits_{\\F_7[x]}(f,g) &= x + 2\n",
"\\end{align*}\n",
"\n",
"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$.\n",
"\n",
"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?\n",
"\n",
"### A criteria for nontrivial GCDs\n",
"\n",
"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$.\n",
"\n",
"If $h\\neq 1$ then there is a nontrivial solution to the equation\n",
"\n",
"\\begin{equation*}\n",
"\\text{$sf+tg=0$ with $\\deg(s)<\\deg(g)$ and $\\deg(t)<\\deg(f)$.}\\tag{$*$}\n",
"\\end{equation*}\n",
"\n",
"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).\n",
"\n",
"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$.\n",
"\n",
"Note that ($*$) can equivalently be written as the following matrix-vector product equation:\n",
"\n",
"\\begin{equation*}\n",
"\\begin{bmatrix}\n",
"f_n & & & & g_m \\\\\n",
"f_{n-1} & f_n & & & \\vdots & g_m \\\\\n",
"\\vdots & f_{n-1} & \\ddots & & g_0 & \\vdots & g_m \\\\\n",
"f_1 & \\vdots & & f_n & & g_0 & \\vdots & g_m \\\\\n",
"f_0 & f_1 & & f_{n-1} & & & g_0 & \\vdots & \\ddots \\\\\n",
" & f_0 & & \\vdots & & & & g_0 & & g_m \\\\\n",
" & & \\ddots & f_1 & & & & & \\ddots & \\vdots \\\\\n",
" & & & f_0 & & & & & & g_0\n",
"\\end{bmatrix}\n",
"\\begin{bmatrix}\n",
"s_{m-1} \\\\\n",
"s_{m-2} \\\\\n",
"\\vdots \\\\\n",
"s_0 \\\\\n",
"t_{n-1} \\\\\n",
"t_{n-2} \\\\\n",
"\\vdots \\\\\n",
"t_0\n",
"\\end{bmatrix} = \n",
"\\begin{bmatrix}\n",
"0 \\\\\n",
"0 \\\\\n",
"\\vdots \\\\\n",
"0\n",
"\\end{bmatrix} \\in F^{n+m}\n",
"\\end{equation*}\n",
"\n",
"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$.\n",
"\n",
"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).\n",
"\n",
"The determinant of the Sylvester matrix of $f$ and $g$ is known as the _resultant_ $\\renewcommand{\\res}{\\operatorname{res}}\\res(f,g)$.\n",
"A matrix is singular if and only if its determinant is $0$, so we can equivalently state this as\n",
"\\begin{equation*}\n",
"\\gcd(f,g)\\neq1 \\Longleftrightarrow \\res(f,g)=0.\n",
"\\end{equation*}\n",
"\n",
"The above takes place over a field $F$ but due to Gauss' theorem it can also be modified to work over $\\Z$:\n",
"\n",
"\\begin{equation*}\n",
"\\gcd\\nolimits_{\\Z[x]}(f,g)\\text{ is nonconstant} \\Longleftrightarrow \\res(f,g)=0.\n",
"\\end{equation*}\n",
"\n",
"### A criteria for choosing primes\n",
"\n",
"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$.\n",
"\n",
"#### Theorem (6.26, Modern Computer Algebra)\n",
"\n",
"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)$.\n",
"\n",
"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$.)\n",
"\n",
"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)$.\n",
"\n",
"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.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The modular algorithm for GCDs\n",
"\n",
"So the prime $p$ must satisfy the following:\n",
"\n",
"1. $p$ does not divide $\\gcd(f_n,g_m)$\n",
"2. $p$ does not divide $\\res(f/h,g/h)$\n",
"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\n",
"\n",
"If $p$ satisfies all three conditions then we can compute $h$, the $\\gcd(f,g)$ over $\\Z$, by:\n",
"\n",
"* Using the Euclidean algorithm to compute $\\gcd(f,g)$ over $\\F_p$\n",
"* Multiplying this computed GCD by $\\gcd(f_n,g_m)$ and reduce the coefficients to be in the symmetric range modulo $p$\n",
"* Return the primitive part of the above polynomial"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Example\n",
"\n",
"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.\n",
"\n",
"First, $p$ must not divide $\\gcd(f_3,g_2)=\\gcd(3,3)=3$. Thus $p\\neq3$\n",
"\n",
"Recall that $f/h=x^2+1$ and $g/h=x-2$, and the Sylvester matrix of these two polynomials is\n",
"\n",
"\\begin{equation*}\n",
"\\begin{bmatrix}\n",
"1 & 1 & 0 \\\\\n",
"0 & -2 & 1 \\\\\n",
"1 & 0 & -2\n",
"\\end{bmatrix}\n",
"\\end{equation*}\n",
"\n",
"which has determinant $(-2)^2+1=5$. Thus $p\\neq5$.\n",
"\n",
"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$.\n",
"\n",
"Thus, the simplest selection is $p=7$.\n",
"\n",
"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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Caveats\n",
"\n",
"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!_\n",
"\n",
"How can we get around this?\n",
"\n",
"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.\n",
"\n",
"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$.\n",
"\n",
"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.\n",
"\n",
"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$."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "SageMath 9.5",
"language": "sage",
"name": "sagemath"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
},
"latex_metadata": {
"author": "Curtis Bright",
"date": "December 7, 2022",
"title": "Computational Mathematics: Handout 11"
}
},
"nbformat": 4,
"nbformat_minor": 2
}