# What are the steps involved in finding the Greatest Common Divisor of two polynomials?

Ultimately I'm trying to define all the steps necessary to go from this toy quartic polynomial modulus:

$$x^4 + 21x^3 + 5x^2 + 7x + 1 \equiv 0 \mod 23$$ to: $$x = 18, 19$$

One of the recommended first steps is to use the polynomial version of Euclid's Algorithm like this:

$$\gcd(x^4 + 21 x^3 + 5 x^2 + 7 x + 1,x^{23}-x)$$

From what I've gathered, the solution is: $$x^2+9x+20$$

For this question, can someone guide me through the steps involved in reaching $x^2+9x+20$? Bonus points if steps to $x = -5, -4$ were provided. :)

I've read many wiki pages, watched several youtube videos, and countless stackexchange posts on the matter and I'm more confused now than when I started. In each example, it seems the initial setup introduces forms of polynomials I have no idea how they're are obtained.

Part of my confusion stems from the fact that polynomial long division is often mentioned as one of the steps and illustrations indicate that when using my $x^{23}-x$, that one of the steps is a polynomial with an x of all powers from 19 to 0.

For example:

$$x^{19} + 2 x^{18} + 22 x^{17} + 4 x^{16} + 21 x^{15} + 4 x^{14} + 14 x^{13} + 18 x^{12} + 9 x^{11} + 10 x^{10} + 19 x^{9} + 22 x^{8} + 8 x^{7} + 16 x^{6} + 3 x^{5} + 9 x^{4} + 21 x^{3} + 6 x^{2} + 2 x + 2$$

I need to use a solution that does not do this because my real problem involves a modulus n of $115792089237316195423570985008687907852837564279074904382605163141518161494337$ and I couldn't even write this out in this universe's lifetime.

This post suggests that Square-and-multiply can be used to quickly find this answer. I am familiar with the Square-and-multiply technique using numbers, but am all thumbs where polynomials are involved.

Any clarity would be helpful here. Steps would be great. Thanks!

By the way, the modulus here is prime so the Chinese remainder theorem and Hensel lifting aren't necessary.

Update

Here is what I'm actually trying to solve. Hope it sheds some light into the right approach with my toy example.

$$\gcd(x^4 + 104904789243076764769272258331008861968773673797810778195683332912883890030960x^3 + 0x^2 + 7476413576101162797801345879216150799700321200725070869821865316357935848858x + 115422971207940030049746909572711733051107198341303901253516904650890299203952, x^{115792089237316195423570985008687907852837564279074904382605163141518161494337} - x)$$

• Your modulus is large, but what about the degrees of the polynomials? storing and computing with $2\times 23$ integers of size $2^{256}$ is not that bad. – user580373 Aug 2 '18 at 13:32
• @Levitikon I can explain how to get $(x^2+12x+15)(x^2+9x+20)$ if you want. – Michael Rozenberg Aug 2 '18 at 13:51
• @spiralstotheleft I updated my question to include an example of the large version. Looks like on the left side, the degrees remain small, but the coefficients will be large. The right side however will have a large degree. – Levitikon Aug 2 '18 at 14:00
• @MichaelRozenberg that would be great! I've surrendered to the fact that I'm going to have to break this problem into several questions. This question was about getting to $x^2+9x+20$. I didn't realize there was a second part there. My next question was going to be about how two distinct zeros modulo 23 were determined from this solution. I think that partially answers that :) – Levitikon Aug 2 '18 at 14:03
• For heaven's sake don't do long division in the first step of Euclid when the other polynomial has ridiculously high degree. That is precisely the step where you are to use square-and-multiply instead. – Jyrki Lahtonen Aug 3 '18 at 21:33

Ok. One more annotated run of Cantor-Zassenhaus :-)

Let $f(x)$ be your quartic. The first step is, indeed, to calculate $\gcd(f(x),x^{23}-x)$. The reason for this is that modulo the prime $p=23$ $$x^{23}-x=\prod_{j=0}^{22}(x-j).$$ By uniqueness of factorization in the polynomial ring $\Bbb{F}_p[x]$ it follows that $$g(x):=\gcd(f(x),x^{23}-x)=\prod_{a\in\Bbb{F}_{23},f(a)=0}(x-a).$$ Implying that any zero $a\in\Bbb{F}_p$ of $f(x)$ is a simple zero of $g(x)$. Further implying that the degree of $g(x)$ tells us the number of distinct zeros in the prime field.

The first step in Euclid's algorithm calls for the remainder of $x^{23}-x$ modulo $f(x)$. This I recommend (at least in the not toy example!) you to calculate with square-and-multiply. First, by repeated squaring \begin{aligned} x&\equiv x\pmod{f},\\ x^2&\equiv x^2\pmod{f},\\ x^4&\equiv x^4-f(x)=2 x^3+18 x^2+16 x+22\pmod{f},\\ x^8&\equiv (2 x^3+18 x^2+16 x+22)^2\equiv 4 x^3+6 x^2+10 x+2\pmod{f},\\ x^{16}&\equiv(4 x^3 + 6 x^2 + 10 x + 2)^2\equiv (15 + 14 x + 17 x^2 + 16 x^3)\pmod{f}. \end{aligned} In the last couple steps we had to do a banal squaring, BUT it was a calculation involving low degree polynomials only (here cubics). This is a theme in all of Cantor-Zassenhaus (and also Berlekamp): arithmetic operations involving low degree polynomials are cheap (in terms of complexity).

Next I want to calculate the remainder of $x^{23}$ modulo $f$. All I need to observe is that $$23=16+4+2+1$$ from writing $p=23$ in binary. Therefore $$x^{23}\equiv x^{16}\cdot x^4\cdot x^2\cdot x\pmod {f}.$$ But we already calculated the remainders of those factors on the right. They are all cubics at worst, so it is cheap to calculate their product modulo $f$. Do it two factors at a time! Then you never need to deal with polynomials of degrees higher than six, and reducing those modulo $f$ is cheap. The result of this calculation is that $$x^{23}-x\equiv 21 + 6 x + 16 x^2\pmod {f}.$$ The next step in Euclid's algorithm calls us to calculate the remainder of $f$ modulo that first remainder $r_1(x)=21 + 6 x + 16 x^2$. This is, again, a calculation involving low degree polynomials only. Today is our lucky day, because it turns out the next remainder is zero. In other words, $r_1(x)$ is a factor of $f$. Normally we would have gotten another remainder $r_2(x)$, then calculated the remainder of $r_1(x)$ modulo $r_2(x)$ etc. But that's all "cheap".

For easier manipulation I normalize and make $r_1(x)$ monic by dividing with its leading coefficient: $$g(x)=\gcd(f(x),x^{23}-x)=16^{-1}r_1(x)=13r_1(x)=x^2+9x+20.$$ An important observation here is that if we use a larger prime $p$ the only step where we had to deal with high degree polynomials was the first one of calculating the remainder of $x^p-x$ modulo $f(x)$. And that could be carried out with $\log_2p$ repeated squarings (each "cheap" individually) and the at most $\log_2p$ multiplications involving low degree polynomials. So to summarize: calculating the remainder of $x^p-x$ modulo $f$ takes at worst $2\log_2p$ cheap operations. This is not a bad complexity at all!

But we are not done. Imagining that in place of $23$ we have a huge prime, we cannot find the roots of $g(x)$ either by testing the candidates in sequence. We need to factor $g(x)$ further. But, we already know that $g(x)$ factors into linear terms because $g(x)$ is a factor of $x^p-x$.

The idea in Cantor-Zassenhaus factorization is to use proper "random" factors of $x^{23}-x$ in its place. A common scheme is to use the set $Q_p$ of quadratic residues modulo $p$. From Euler's formula we know that all the elements of $Q_p$ are zeros of the polynomial $x^{(p-1)/2}-1$. In our toy case, of $x^{11}-1$. As above, $$\gcd(g(x),x^{11}-1)=\prod_{a\in Q_p,g(a)=0}(x-a).$$ So here $x^{11}-1$ (really $x^{(p-1)/2}-1$) can be a ridiculously high degree polynomial. But, because it is a binomial, its remainder modulo $g(x)$ can be calculated by square-and-multiply "cheaply". And then we are, again, left with low degree polynomials making the rest of Euclid easy. I spare you the details and simply tell that $$\gcd(x^2+9x+20,x^{11}-1)=x+5.$$ The following things need to be noted:

• This tells us that $x+5$ is a factor of $g(x)$, hence of $f(x)$. And hence that $x=-5=18$ is a zero of $f$.
• By polynomial division (cheap!) $g(x)/(x+5)=(x+4)$, so we have fully factored $g(x)$ at this point, and therefore fully solved the toy example.
• Here we were a bit lucky. It could easily have happened that neither of the two zeros of $g(x)$ were $\in Q_p$, in which case we would have gotten $\gcd(g(x),x^{11}-1)=1$. Equally bad would have been that both of the zeros turned out to be in $Q_p$. For then we would have gotten $g(x)=\gcd(g(x),x^{11}-1)$ and no lower degree factors of $g(x)$ would have come to light.
• To remedy the problem of the previous bullet Cantor-Zassenhaus calls for the calculation of $\gcd(g(x),(x-a)^{11}-1)$ for several randomly chosen elements $a\in\Bbb{F}_p$. This gcd spits out a polynomial that has as its zeros exactly those zeros $x=z$ of $g(x)$ such that $z+a\in Q_p$. Doing this for many distinct choices of $a$ is guaranteed to eventually give a proper factor of $g(x)$! The reason is that the two random events: 1) $z$ is a quadratic residue and 2) $z+a$ is a quadratic residue both occur at close to 50% probability, and there is next to no correlation. See a recent answer of mine for the details of how this follows from Weil's bounds on character sums.
• This is incredible! Thanks @Jyrki for this insightful annotation. I'm sure once I fully digest this the path will be clear. – Levitikon Aug 4 '18 at 15:02
• Could you point me to some info around banal squaring? It might be a typo as I can't seem to find any info on the topic. – Levitikon Aug 4 '18 at 15:04
• @Levitikon "banal" = "high school polynomial multiplication. Actually you can speed that up a little by having precalculated the remainders of $x^4$, $x^5$ and $x^6$ modulo $f(x)$. You see, when you multiply two remainders, both cubic at worst, you get a polynomial of degree six. Instead of performing long division by $f$ you can at that point simply replace the eventual terms of degrees four, five and six, and merrily dance to the next step. In fact, this helps you in the square-and-multiply step also! – Jyrki Lahtonen Aug 4 '18 at 16:11
• More generally, when applying Cantor-Zassenhaus to a polynomial $f(x)$ of degree $n$, you ever need to use polynomials of degrees higher than $2(n-1)$. If you checked the comments under the main question you see how Mathematica couldn't cope when I started increasing the degree of the polynomial as well as using a huge prime. Both square-and-multiply and Euclid become more time consuming when $\deg f$ grows. – Jyrki Lahtonen Aug 4 '18 at 16:13
• A little stuck at $\gcd(x^2+9x+20,x^{11}-1)$. Everything before this point was about solving $\gcd(x^4 + 21 x^3 + 5 x^2 + 7 x + 1,x^{23}-x)$. The same steps applied to this second $gcd$ don't seem to produce your results. I'm sure it has to do with switching from $a\in\Bbb{F}_{23}$ to $a\in Q_p$. Could you briefly point out how any of the gcd steps change as a result of using 23-adic? I'm guessing all the $\mod f$ and $\mod p$ of the coefficients should be omitted, and the monic division should be true division instead of multiplication of modular inverse of the leading coefficient? – Levitikon Aug 7 '18 at 12:49

gp-pari says the same, coprime in $\mathbb Q[x]$ There are algorithms for factoring mod p, such as https://en.wikipedia.org/wiki/Berlekamp%27s_algorithm My guess is that the gcd problem you set is not quite right for Berlekamp. https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields

? factormod(  x^4 + 21 * x^3 + 5 * x^2 + 7 * x + 1 , 23 )
%4 =
[                   Mod(1, 23)*x + Mod(4, 23) 1]

[                   Mod(1, 23)*x + Mod(5, 23) 1]

[Mod(1, 23)*x^2 + Mod(12, 23)*x + Mod(15, 23) 1]

?


=====================================================

? gcd( x^23 - x , x^4 + 21 * x^3 + 5 * x^2 + 7 * x + 1 )
%2 = 1
? bezout( x^23 - x , x^4 + 21 * x^3 + 5 * x^2 + 7 * x + 1 )
%3 = [601562001455448034332304838/601559367722392947539135415*x^3 + 12632734286485513487226173386/601559367722392947539135415*x^2 + 1002598481708396788521486307/200519789240797649179711805*x + 4210915650874735492265878157/601559367722392947539135415,
-601562001455448034332304838/601559367722392947539135415*x^22 + 67744078895233752228212/601559367722392947539135415*x^21 - 1408063504750102699727183/601559367722392947539135415*x^20 + 29248972518676735993385492/601559367722392947539135415*x^19 - 6100312465279544295752063/601559367722392947539135415*x^18 - 8349600368157764610272068/601559367722392947539135415*x^17 + 2508425931723729040502482/601559367722392947539135415*x^16 + 2524272012870587277687167/601559367722392947539135415*x^15 - 1004326886517081466286378/601559367722392947539135415*x^14 - 32168547887068011724661/26154755118364910762571105*x^13 + 380713040221417010130502/601559367722392947539135415*x^12 + 210425355112047121912952/601559367722392947539135415*x^11 - 139034562125043256867703/601559367722392947539135415*x^10 - 57515651081682016589308/601559367722392947539135415*x^9 + 49310957334791769192817/601559367722392947539135415*x^8 + 14864731141038606058352/601559367722392947539135415*x^7 - 17070020938952200196618/601559367722392947539135415*x^6 - 3514266249057193923193/601559367722392947539135415*x^5 + 5785620602900061768862/601559367722392947539135415*x^4 + 698714016011467787837/601559367722392947539135415*x^3 - 1931212668388574729918/601559367722392947539135415*x^2 + 76817984859491930252/601559367722392947539135415*x + 1, 1]


this is from my own program, I get your two polynomials are coprime

$$\left( x^{23} - x \right) \left( \frac{ 601562001455448034332304838 x^{3} + 12632734286485513487226173386 x^{2} + 3007795445125190365564458921 x + 4210915650874735492265878157 }{ 601559367722392947539135415 } \right) - \left( x^{4} + 21 x^{3} + 5 x^{2} + 7 x + 1 \right) \left( \frac{ 601562001455448034332304838 x^{22} - 67744078895233752228212 x^{21} + 1408063504750102699727183 x^{20} - 29248972518676735993385492 x^{19} + 6100312465279544295752063 x^{18} + 8349600368157764610272068 x^{17} - 2508425931723729040502482 x^{16} - 2524272012870587277687167 x^{15} + 1004326886517081466286378 x^{14} + 739876601402564269667203 x^{13} - 380713040221417010130502 x^{12} - 210425355112047121912952 x^{11} + 139034562125043256867703 x^{10} + 57515651081682016589308 x^{9} - 49310957334791769192817 x^{8} - 14864731141038606058352 x^{7} + 17070020938952200196618 x^{6} + 3514266249057193923193 x^{5} - 5785620602900061768862 x^{4} - 698714016011467787837 x^{3} + 1931212668388574729918 x^{2} - 76817984859491930252 x - 601559367722392947539135415 }{ 601559367722392947539135415 } \right) = \left( 1 \right)$$

• Thanks Will! One big thing that stands out is $x^{22} + x^{21} + x^{20} + x^{19} + ...$ all the way to $x^0$. That's my expressed great concern. Unless i'm missing something, like a cutoff somewhere, I can't do that when starting with $x^{115792089237316195423570985008687907852837564279074904382605163141518161494337}$ and finish writing all that out in this universe's lifetime. WolframAlpha can compute that large expression in a fraction of a second so I know they aren't computing every power. – Levitikon Aug 2 '18 at 23:53