# How does Mathematica solve $f(x)\equiv 0\pmod p$?

I input

Solve[35345 x^47 + 1345326 x^2 + 134636246 x + 13451345 == 0, x, Modulus ->2037035976334486086268445688409378161051468393665936250636140449354381299763336706183399101]

to Mathematica and the output is

x-> 341367685451396522792304545319361914942948733456457424375640931820820238770112149829924341

I wonder how Mathematica gets this result? Are there special algorithms to solve this problem?

• It's probably off topic to try to guess what algorithm software is using. Someone might be able to tell known techniques that can solve such equations, but the nature of algorithms is that two different algorithms should give the same result... Aug 10, 2014 at 15:13
• @Thomas Andrews Thank you! I'm not going to guess what algorithm Mathematica is using, but to learn a good algorithm to solve this porblem. I guess Mathematica would use the fast one.. Aug 10, 2014 at 15:17
• Yeah, if the GCD is not linear, you've still got a problem. Heck, $f(x)$ might properly divide $x^p-x$, which just means that there are $47$ distinct roots. :) @JyrkiLahtonen Aug 10, 2014 at 15:26
• I checked a few things. Mathematica cannot factor $p-1$ for the given prime, but yet it can list all the solutions of $x^{100}=1$ (clearly $100\mid (p-1)$) in less than five seconds. This is not a true test, because those can be found non-deterministically by raising random integers modularly to the power $(p-1)/100$ by square-and-multiply. Also the solution set is a multiplicative subgroup. Aug 10, 2014 at 15:28
• Yeah, I edited the comment. Equations of the type $x^m+x=a$ give the following results. With $m=300$ about 20 seconds, with $m=3000$ I gave up after 3 minutes had passed. Aug 10, 2014 at 15:39

The factorization of a polynomial over a finite field is a well-understood problem.

The key idea is that $$x^p-x = \prod_{\xi\in\mathbb{F}_p}(x-\xi),$$ hence the greatest common divisor between $x^p-x$ and $f(x)$ gives all the roots of $f$ in $\mathbb{F}_p$, and there are many efficient methods to compute the $\gcd$ of two polynomials, the classical ones relying on the theory of (sub-)resultants. The complexity essentially depends on a fixed power of the degree of $f$ and on a fixed power of $\log(p)$, hence the problem is very affordable even if the field has a huge characteristic.

Addendum. To recover the roots from $g(x)=\gcd(x^p-x,f(x))$ you can for first get rid of multiple roots by considering $h(x)=g(x)/\gcd(g(x),g'(x))$.

Then you can use a "divide-and-conquer" strategy base on the fact that $x^{\frac{p+1}{2}}-x$ is the product of all the linear factors $(x-\xi)$ with $\xi$ being a quadratic residue, hence you can "split" $h(x)$ into two polynomials, the first one having quadratic residues as roots, the second one having non-quadratic residues as roots. If you know how $p-1$ factors, you can continue this way by separating (now I am assuming that $3\mid(p-1)$) the roots in any factor according to the value of $$\xi^{\frac{p-1}{3}}\in\{1,\omega,\omega^2\}\subset\mathbb{F}_p,$$ till you get the complete list of linear factors.

The probabilistic Cantor-Zassenhaus algorithm do essentially the same without the knowledge of the factorization of $(p-1)$.

• $\gcd(f(x),x^p-x)$ only gives you a minimal polynomial with the same roots, it doesn't give "all the roots." For example, if $\gcd(f(x),x^p-x)=f(x)$, all you know is that $f(x)$ has $\deg f$ distinct roots. Aug 10, 2014 at 15:52
• @ThomasAndrews: I hope my last addendum solves this issue, too. Aug 10, 2014 at 16:05
• Addendum on the addendum: the Cantor-Zassenhaus' algorithm. Aug 10, 2014 at 16:14
• That's interesting. Thanks, Jack! Aug 10, 2014 at 16:16
• They are comparable, but for practical purposes the probabilistic algorithm performs better. A similar "struggle" happens for primality testing, for which the (probabilistic) Miller-Rabin's algorithm outperforms the (deterministic) AKS-algorithm. Aug 10, 2014 at 16:48

As I did a little bit of homework on Cantor-Zassenhaus algorithm I will write up a toy example that can be done with paper and pencil (if you have enough patience and paper). Explaining bits and pieces while I go.

Task: Find the zeros of $f(x)=x^4+13x^2+6x+5$ in the field $F=\Bbb{F}_{17}$.

The first step is to calculate the greatest common divisor $\gcd(f(x),x^{17}-x)$ in the ring $F[x]$. This is fast enough, and it turns out (I set it up that way) that $f(x)\mid x^{17}-x$. This has the following two important consequences.

• All the irreducible factors of $f(x)$ in $F[x]$ are linear. This is because we know from basic finite field theory that $x^{17}-x$ is the product of the monic linear polynomials.
• All the irreducible factors of $f(x)$ are simple. No repeated factors are possible, because $x^{17}-x$ does not have any.

Consequently $f(x)=(x-a_1)(x-a_2)(x-a_3)(x-a_4)$ for some pairwise distinct $a_1,a_2,a_3,a_4\in F$. We could easily find them by trial and error, but let's not. I want to describe the C-Z idea that also works for primes $q$ much larger than $p=17$. At this point we know (from Chinese remainder theorem) that \begin{aligned} R=F[x]/\langle f(x)\rangle&\cong F[x]/\langle x-a_1\rangle\oplus\cdots F[x]/\langle x-a_4\rangle\\ &\cong F\oplus F\oplus F\oplus F. \end{aligned} The idea is that if we take a random non-zero element $r\in R$, we can easily calculate $r^{(p-1)/2}-1=r^8-1$. If in the above isomorphism the image of $r$ is $(r_1,r_2,r_3,r_4)\in F^4$ is such that some of the (unknown) components $r_i$, $i=1,2,3,4,$ are non-zero quadratic residues, and some are not, then some of the components of $r^8-1,$ i.e.$r_i^8-1, i=1,2,3,4,$ are non-zero and some are not. This happens with probability at least $\approx 1/2$ (for larger primes). When this happens it means that some but not all of the factors $(x-a_i)$ are also factors $r^{(p-1)/2}-1$. Calculating the greatest common divisor of $f(x)$ and the lowest degree polynomial in the coset of $r^8-1$ modulo $f(x)$ then reveals a non-trivial factor. Rinse. Repeat.

In our example case let us try $r=x$ for simplicity. Then it turns out that $$x^8-1\equiv 12+15x+9x^2+3x^3\pmod{f(x)}.$$ We see (Euclid's algorithm)) that with the remainder $3f_1, f_1(x)=x^3+3x^2+5x+4$, we actually have $f_1(x)\mid f(x)$. A long division shows that $$f(x)/f_1(x)=x+14.$$ This shows that $a_4=3$ is one of the zeros. This step was already explained in Jack's excellent answer. We can conclude that $3$ is the only zero that is NOT a quadratic residue modulo $17$. After all, if $a$ is a square in $F$, then $x-a\mid x^8-1$.

But we are not done, so let's use another $r$. Let's try $r=x+1+\langle f(x)\rangle$. This time $r^8-1\equiv 6+3x+7x^2+16x^3\pmod{f(x)}$. Another run of Euclid's algorithm gives $$f_2(x)=\gcd(6+3x+7x^2+16x^3,f(x))=7+6x+x^2.$$ This is a quadratic factor, so our luck is in. Yet another run of Euclid's algorithm shows that $$\gcd(f_2(x),f_1(x))=x+9,$$ so we see that $a_3=8$ is one of the zeros. By polynomial division (or using Vieta relations) we see that $a_4=3$ is the other zero of $f_2(x)$, so two zeros remain unknown.

We could keep going, but it is worth pointing out that it is better (=computationally simpler) to start working with $$f_3(x)=f(x)/f_2(x)=8+11x+x^2$$ instead of $f(x)$.

Let's pick $r=(x+2)$. A calculation shows that $$(x+2)^8-1\equiv 2+16x=2-x\pmod{f_3(x)}.$$ As $\gcd(x-2,f_3(x))=x-2$, this reveals the zero $a_2=2$. Vieta relations then reveal the final zero $a_1=4$.

Observe that powers in rings such as $R$ can be computed fast with the aid of square-and-multiply. The run time of Euclid's algorithm depends on the degree of the polynomials. Here only in the first step of calculating $\gcd(f(x),x^p-x)$ is there a high degree polynomial. But because it is a binomial, its remainder can be computed (as pointed out by Thomas Andrews) very quickly, again by square-and-multiply.

The non-deterministic element here is in the choice of $r$. If at first you don't succeed, pick another. The factors I used in the above toy example don't look random, but I swear :-) I didn't design the example so that they would work (but in retrospect I could have - I only understood this afterwards).

• It may be worth pointing out that the choice $r=x-a$, $a$ a random constant, produces a factor that has as zeros those roots $a_i$ with the property that $a_i+a$ is a square. With the aid of character sums involving the Legendre character it follows easily that the subsets of roots we get in this way are independent from each other in the sense that two such tests have correlation that $\to0$ when $p\to\infty$. Aug 17, 2014 at 6:27
• @Levitikon The first step in the Euclidean algorithm is to compute the remainder of $x^{23}$ (or more generally of $x^p$ for some large prime $p$) modulo your quartic. This is easy with the good ole square-and-multiply. Notice that in all the iteration there is just one round of Euclid involving a very high degree polynomial. And that high degree polynomial has very few terms, so it can be managed quickly (not necessarily by paper and pencil, but still). Jul 31, 2018 at 15:32
• Note, if your modulus is not a prime, then you need other techniques also. Look up Hensel lifting and Chinese remainder theorem. Jul 31, 2018 at 15:33
• Luckily in ECC the modulus is always prime. So I've crash coursed on the square-and-multiply technique and wrote a scripting function that works great for integers. Now, how does it work with polynomials? Can you help me take the first step and using this $gcd(x^4+21x^3+5x^2+7x+1,x^23−x)$, can you show me what the first couple operations are? What is the coefficient and what is the exponent of the first step and what do you do with the result? Aug 1, 2018 at 23:09
• The idea is to calculate $\gcd(f(x),r^8-1)$. If $r=x-a$ then the gcd you get has as zeros those roots $x=b$ of $f(x)$ such that $a+b$ is a quadratic residue. A point is that for two different choices of $a$ we get independent information about the zeros. It can be shown that if $b+a_1$ is a quadratic residue then $b+a_2$ is a quadratic residue with probability close to $1/2$. When you get different random factors of $f(x)$ you quickly get it down to linear factors. And then you have found all the zeros. Aug 3, 2018 at 20:27