1. History

The method was proposed by Elwyn Berlekamp in his 1970 work^[1] on polynomial factorization over finite fields. His original work lacked a formal correctness proof^[2] and was later refined and modified for arbitrary finite fields by Michael Rabin.^[2] In 1986 René Peralta proposed a similar algorithm^[3] for finding square roots in ${\displaystyle {\mathbb{Z}}_p}$.^[4] In 2000 Peralta's method was generalized for cubic equations.^[5]

2. Statement of Problem

Let ${\displaystyle p}$ be an odd prime number. Consider the polynomial ${\displaystyle f(x)=a_0+a_{1}x+\cdots +a_n{x}^n}$ over the field ${\displaystyle {\mathbb{Z}}_p}$ of remainders modulo ${\displaystyle p}$. The algorithm should find all ${\displaystyle \lambda }$ in ${\displaystyle {\mathbb{Z}}_p}$ such that ${\displaystyle f(\lambda )=0}$ in ${\displaystyle {\mathbb{Z}}_p}$.^[2][6]

3. Algorithm

3.1. Randomization

Let ${\displaystyle f(x)=(x-{\lambda }_1)(x-{\lambda }_2)\cdots (x-{\lambda }_n)}$. Finding all roots of this polynomial is equivalent to finding its factorization into linear factors. To find such factorization it is sufficient to split the polynomial into any two non-trivial divisors and factorize them recursively. To do this, consider the polynomial ${\displaystyle f_z(x)=f(x-z)=(x-{\lambda }_1-z)(x-{\lambda }_2-z)\cdots (x-{\lambda }_n-z)}$ where ${\displaystyle z}$ is some any element of ${\displaystyle {\mathbb{Z}}_p}$. If one can represent this polynomial as the product ${\displaystyle f_z(x)=p_0(x)p_1(x)}$ then in terms of the initial polynomial it means that ${\displaystyle f(x)=p_0(x+z)p_1(x+z)}$, which provides needed factorization of ${\displaystyle f(x)}$.^[1][6]

3.2. Classification of ${\displaystyle {\mathbb{Z}}_p}$ Elements

Due to Euler's criterion, for every monomial ${\displaystyle (x-\lambda )}$ exactly one of following properties holds:^[1]

1. The monomial is equal to ${\displaystyle x}$ if ${\displaystyle \lambda =0}$,
2. The monomial divides ${\displaystyle g_0(x)=(x^{(p-1)/2}-1)}$ if ${\displaystyle \lambda }$ is quadratic residue modulo ${\displaystyle p}$,
3. The monomial divides ${\displaystyle g_1(x)=(x^{(p-1)/2}+1)}$ if ${\displaystyle \lambda }$ is quadratic non-residual modulo ${\displaystyle p}$.

Thus if ${\displaystyle f_z(x)}$ is not divisible by ${\displaystyle x}$, which may be checked separately, then ${\displaystyle f_z(x)}$ is equal to the product of greatest common divisors ${\displaystyle gcd(f_z(x);g_0(x))}$ and ${\displaystyle gcd(f_z(x);g_1(x))}$.^[6]

3.3. Berlekamp's Method

The property above leads to the following algorithm:^[1]

1. Explicitly calculate coefficients of ${\displaystyle f_z(x)=f(x-z)}$,
2. Calculate remainders of ${\displaystyle x,x^2,x^{2^2},x^{2^3},x^{2^4},\dots ,x^{2^{\lfloor {\log }_{2}p\rfloor }}}$ modulo ${\displaystyle f_z(x)}$ by squaring the current polynomial and taking remainder modulo ${\displaystyle f_z(x)}$,
3. Using exponentiation by squaring and polynomials calculated on the previous steps calculate the remainder of ${\displaystyle x^{(p-1)/2}}$ modulo ${\displaystyle f_z(x)}$,
4. If ${\displaystyle x^{(p-1)/2}\not\equiv \pm 1(\mathrm{mod}{f}_z(x))}$ then ${\displaystyle gcd}$ mentioned above provide a non-trivial factorization of ${\displaystyle f_z(x)}$,
5. Otherwise all roots of ${\displaystyle f_z(x)}$ are either residues or non-residues simultaneously and one has to choose another ${\displaystyle z}$.

If ${\displaystyle f(x)}$ is divisible by some non-linear primitive polynomial ${\displaystyle g(x)}$ over ${\displaystyle {\mathbb{Z}}_p}$ then when calculating ${\displaystyle gcd}$ with ${\displaystyle g_0(x)}$ and ${\displaystyle g_1(x)}$ one will obtain a non-trivial factorization of ${\displaystyle f_z(x)/g_z(x)}$, thus algorithm allows to find all roots of arbitrary polynomials over ${\displaystyle {\mathbb{Z}}_p}$.

3.4. Modular Square Root

Consider equation ${\displaystyle x^2\equiv a(\mathrm{mod}p)}$ having elements ${\displaystyle \beta }$ and ${\displaystyle -\beta }$ as its roots. Solution of this equation is equivalent to factorization of polynomial ${\displaystyle f(x)=x^2-a=(x-\beta )(x+\beta )}$ over ${\displaystyle {\mathbb{Z}}_p}$. In this particular case problem it is sufficient to calculate only ${\displaystyle gcd(f_z(x);g_0(x))}$. For this polynomial exactly one of the following properties will hold:

1. GCD is equal to ${\displaystyle 1}$ which means that ${\displaystyle z+\beta }$ and ${\displaystyle z-\beta }$ are both quadratic non-residues,
2. GCD is equal to ${\displaystyle f_z(x)}$which means that both numbers are quadratic residues,
3. GCD is equal to ${\displaystyle (x-t)}$which means that exactly one of these numbers is quadratic residue.

In the third case GCD is equal to either ${\displaystyle (x-z-\beta )}$ or ${\displaystyle (x-z+\beta )}$. It allows to write the solution as ${\displaystyle \beta =(t-z)(\mathrm{mod}p)}$.^[1]

3.5. Example

Assume we need to solve the equation ${\displaystyle x^2\equiv 5(\mathrm{mod}11)}$. For this we need to factorize ${\displaystyle f(x)=x^2-5=(x-\beta )(x+\beta )}$. Consider some possible values of ${\displaystyle z}$:

1. Let ${\displaystyle z=3}$. Then ${\displaystyle f_z(x)=(x-3{)}^2-5=x^2-6x+4}$, thus ${\displaystyle gcd(x^2-6x+4;x^5-1)=1}$. Both numbers ${\displaystyle 3\pm \beta }$ are quadratic non-residues, so we need to take some other ${\displaystyle z}$.

1. Let ${\displaystyle z=2}$. Then ${\displaystyle f_z(x)=(x-2{)}^2-5=x^2-4x-1}$, thus ${\displaystyle gcd(x^2-4x-1;x^5-1)\equiv x-9(\mathrm{mod}11)}$. From this follows ${\displaystyle x-9=x-2-\beta }$, so ${\displaystyle \beta \equiv 7(\mathrm{mod}11)}$ and ${\displaystyle -\beta \equiv -7\equiv 4(\mathrm{mod}11)}$.

A manual check shows that, indeed, ${\displaystyle 7^2\equiv 49\equiv 5(\mathrm{mod}11)}$ and ${\displaystyle 4^2\equiv 16\equiv 5(\mathrm{mod}11)}$.

4. Correctness Proof

The algorithm finds factorization of ${\displaystyle f_z(x)}$ in all cases except for ones when all numbers ${\displaystyle z+{\lambda }_1,z+{\lambda }_2,\dots ,z+{\lambda }_n}$ are quadratic residues or non-residues simultaneously. According to theory of cyclotomy,^[7] the probability of such an event for the case when ${\displaystyle {\lambda }_1,\dots ,{\lambda }_n}$ are all residues or non-residues simultaneously (that is, when ${\displaystyle z=0}$ would fail) may be estimated as ${\displaystyle 2^{-k}}$ where ${\displaystyle k}$ is the number of distinct values in ${\displaystyle {\lambda }_1,\dots ,{\lambda }_n}$.^[1] In this way even for the worst case of ${\displaystyle k=1}$ and ${\displaystyle f(x)=(x-\lambda {)}^n}$, the probability of error may be estimated as ${\displaystyle 1/2}$ and for modular square root case error probability is at most ${\displaystyle 1/4}$.

5. Complexity

Let a polynomial have degree ${\displaystyle n}$. We derive the algorithm's complexity as follows:

1. Due to the binomial theorem $\text{Undefined control sequence binom}$, we may transition from ${\displaystyle f(x)}$ to ${\displaystyle f(x-z)}$ in ${\displaystyle O(n^2)}$ time.
2. Polynomial multiplication and taking remainder of one polynomial modulo another one may be done in ${\displaystyle O(n^2)}$, thus calculation of ${\displaystyle x^{2^k}mod{f}_z(x)}$ is done in ${\displaystyle O(n^2\log p)}$.
3. Binary exponentiation works in ${\displaystyle O(n^2\log p)}$.
4. Taking the ${\displaystyle gcd}$ of two polynomials via Euclidean algorithm works in ${\displaystyle O(n^2)}$.

Thus the whole procedure may be done in ${\displaystyle O(n^2\log p)}$. Using the fast Fourier transform and Half-GCD algorithm,^[8] the algorithm's complexity may be improved to ${\displaystyle O(n\log n\log pn)}$. For the modular square root case, the degree is ${\displaystyle n=2}$, thus the whole complexity of algorithm in such case is bounded by ${\displaystyle O(\log p)}$ per iteration.^[6]