Solution: Testing a system of polynomials for roots over $\mathbb{F}_p$

This is a well-studied problem. An algorithm with the right asymptotics is given in this paper by Ming-Deh Huang and Yiu-Chung Wong. Unfortunately the paper is behind a paywall (it was published in 1999, just before the use of preprint servers like arXiv and ePrint became widespread).

The paper uses some heavy algebraic geometry; the full details are beyond the level of this blog. The best we can do here is to give an impressionistic sketch of the main ideas.

The follow-up problem

Is there a solution that runs in polynomial time in $n$, $m$, $d$, and $\log p$?

No.

This is easy. Polynomial satisfiability is known to be NP-complete (if you haven’t seen this before, it’s a good exercise!), so (assuming P does not equal NP) there can’t be a polynomially fast algorithm for the problem.

Warmup: Solving polynomials in one variable

Suppose you have a polynomial equation $f(x) = 0$ in one variable over a finite field $\mathbb{F}_p$. For this post, we’ll assume it’s a prime field, meaning that $p$ is prime.

For intuition, you should assume that $p$ is large. In cryptography we often work with primes of size about $2^{256}$. The degree of $f$ will be much smaller – imagine it’s somewhere between $2$ and $1000$.

In this setting there is a standard algorithm, called Berlekamp’s algorithm, to find all the roots over $\mathbb{F}_p$.

First off, a bit of algebra. In the finite field $\mathbb{F}_p$ the polynomial $x^p - x$ plays a magical role: it has a root at every single element of the field. (This fact is known as Fermat’s Little Theorem.) In other words, there is a factorization

$$x^p - x = x (x-1) (x-2) \cdots (x - p + 1).$$

Suppose we can compute

$$g(x) = \gcd(x^p - x, f(x)).$$

Then $g$ and $f$ will have all the same roots in $\mathbb{F}_p$, but $g$ will have no more roots. In other words, if $f$ has exactly $k$ roots in $\mathbb{F}_p$, then $g$ will have degree $k$. So just computing $g$ tells you how many roots $f$ has in $\mathbb{F}_p$.

How to compute $g$? The standard method is to use the Euclidean algorithm. In the first step of the Euclidean algorithm, you have to compute $(x^p - x)$ modulo $f(x)$. This seems hard at first glance (remember, $p$ is really big!), but in fact it can be done by repeated squaring. The idea is to compute $x, x^2, x^4, x^8...$ modulo $f(x)$, squaring and reducing modulo $f(x)$ at each step. So you can finish the whole calculation in only about $2 \log p$ calculations (multiplications and reductions) with polynomials of about the same degree as $f$.

OK so now we have this polynomial $g$ with the same roots as $f$, and in fact the degree $k$ of $g$ tells us how many roots $f$ has. But how can we find them?

The idea is to recursively factor $g$. Suppose we can write $g$ as $g(x) = g_1(x) g_2(x)$, where $g_1$ and $g_2$ have smaller degree than $g$. If we keep doing this, eventually we will get down to linear polynomials, which are easy to solve.

For the factorization, we’ll use a randomized algorithm. We’ll choose a random $a \in \mathbb{F}_p$ and split the roots $r$ of $g$ according to whether or not $(r - a)$ is a square modulo $p$. If $(r - a)$ is a square, $r$ will be a root of $g_1$; if not, it will be a root of $g_2$. Since about half of elements modulo $p$ are squares, we expect that $g_1$ and $g_2$ will each have about half the degree of $g$.

In the event that the $(r - a)$’s are either all squares or all not squares, one of $g_1$ or $g_2$ is equal to $g$, and we don’t reduce the degree at all. But that’s okay! If this happens, just choose another random $a$ and try again. One can prove that this bad event will happen with probability at most $1/2$, so if you try a handful of $a$’s the algorithm will break $g$ down with very high probability.

OK, so how do we find $g_1$ and $g_2$? Simple! The polynomial $x^{(p-1)/2} + 1$ has roots at all non-squares and no squares. So just take

$$g_1(x) = \gcd(x^{(p-1)/2} + 1, g(x))$$

and

$$g_2(x) = g(x) / g_1(x).$$

You compute the gcd $g_1(x)$ by repeated squaring, as before.

This is a randomized algorithm; it runs in polynomial time in $\log p$, and it’s fast enough to use in practice.

Back to the problem: One irreducible polynomial

Suppose you just have a single polynomial

$$f(x_1, \ldots, x_n) = 0.$$

To start with, let’s assume that $f$ is what they call absolutely irreducible. This means that you can’t factor $f$ as a product of polynomials

$$f(x_1, \ldots, x_n) = f_1(x_1, \ldots, x_n) f_2(x_1, \ldots, x_n),$$

even if you are allowed to take coefficients in an extension field.

In this case there is a theorem (assuming $p$ is very large, but if not we can just do a brute-force search anyway) which guarantees that $f$ has a large number of roots. In other words: Just plug in random values $(x_1, \ldots, x_{n-1})$ for all but one of the variables, and solve for the last one using Berlekamp’s algorithm.

The theorem gives a guarantee that this will succeed with some probability (roughly $1 / \deg f$). So… just keep trying values until you find a root.

One reducible polynomial

Now suppose you have a polynomial $f$ that has a factorization: In other words, you want to solve

$$f_1(x_1, \ldots, x_n) f_2(x_1, \ldots, x_n) = 0.$$

Well, this is easy!
You just have to find a root for either $f_1$ or $f_2$, and you’re done.

Well we’re sort of hiding a difficulty here. How do you test whether a polynomial is reducible? And if it is, how do you find a factorization? You might remember that factoring integers into prime factors is very hard! Fortunately for us, it’s easier with polynomials, and there are known algorithms to solve the problem. But we won’t get into it here.

Several polynomials

What happens if you have a system of several polynomials? This is where the heavy algebraic geometry comes in.

As an example, imagine a system of 2 equations in 3 variables. Unless the equations are somehow redundant, 2 constraints leave only 1 degree of freedom: the solution set will be a curve in three dimensions.

The trick is to eliminate one of the variables from the two equations, reducing it to a single equation in two variables. Geometrically, this is like projecting the curve down to a plane. Once you do the elimination, you’re back in the setting above, with one equation in two variables, and you’re good to go.

There are standard techniques for this sort of elimination, based on the idea of a polynomial resultant.

The paper does a lot of technical work to avoid unlikely edge cases. (For example: what if the projection causes the solution space to collapse in dimension? What if the solution space has several pieces, some are surfaces, some are curves, and some are points? What if…?) The details are far too intricate for a blog post.

Conclusion

We found an algorithm in the literature with the right asymptotics. It seems that this algorithm isn’t fast enough to be helpful with the concrete problem at hand; we’re continuing to investigate algorithms that will be useful for the concrete systems we’re interested in (roughly, systems of polynomials constraints coming from SNARKs).