Berlekamp–Rabin algorithm

In number theory, Berlekamp's root finding algorithm, also called the Berlekamp–Rabin algorithm, is the probabilistic method of finding roots of polynomials over the field ${\displaystyle \mathbb {F} _{p}}$ with ${\displaystyle p}$ elements. The method was discovered by Elwyn Berlekamp in 1970^[1] as an auxiliary to the algorithm for polynomial factorization over finite fields. The algorithm was later modified by Rabin for arbitrary finite fields in 1979.^[2] The method was also independently discovered before Berlekamp by other researchers.^[3]

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^[4] for finding square roots in ${\displaystyle \mathbb {F} _{p}}$.^[5] In 2000 Peralta's method was generalized for cubic equations.^[6]

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

Let ${\textstyle 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 ${\textstyle f_{z}(x)=f(x-z)=(x-\lambda _{1}-z)(x-\lambda _{2}-z)\cdots (x-\lambda _{n}-z)}$ where ${\displaystyle z}$ is some element of ${\displaystyle \mathbb {F} _{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][7]

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 ${\textstyle g_{0}(x)=(x^{(p-1)/2}-1)}$ if ${\displaystyle \lambda }$ is quadratic residue modulo ${\displaystyle p}$,
3. The monomial divides ${\textstyle 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))}$.^[7]

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 ${\textstyle x,x^{2},x^{2^{2}},x^{2^{3}},x^{2^{4}},\ldots ,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 ${\textstyle x^{(p-1)/2}}$ modulo ${\textstyle f_{z}(x)}$,
4. If ${\textstyle x^{(p-1)/2}\not \equiv \pm 1{\pmod {f_{z}(x)}}}$ then ${\displaystyle \gcd }$ mentioned below 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 {F} _{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 {F} _{p}}$.

Consider equation ${\textstyle x^{2}\equiv a{\pmod {p}}}$ having elements ${\displaystyle \beta }$ and ${\displaystyle -\beta }$ as its roots. Solution of this equation is equivalent to factorization of polynomial ${\textstyle f(x)=x^{2}-a=(x-\beta )(x+\beta )}$ over ${\displaystyle \mathbb {F} _{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 ${\textstyle \beta =(t-z){\pmod {p}}}$.^[1]

Assume we need to solve the equation ${\textstyle x^{2}\equiv 5{\pmod {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 ${\textstyle \gcd(x^{2}-4x-1;x^{5}-1)\equiv x-9{\pmod {11}}}$. From this follows ${\textstyle x-9=x-2-\beta }$, so ${\displaystyle \beta \equiv 7{\pmod {11}}}$ and ${\textstyle -\beta \equiv -7\equiv 4{\pmod {11}}}$.

A manual check shows that, indeed, ${\textstyle 7^{2}\equiv 49\equiv 5{\pmod {11}}}$ and ${\textstyle 4^{2}\equiv 16\equiv 5{\pmod {11}}}$.

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},\ldots ,z+\lambda _{n}}$ are quadratic residues or non-residues simultaneously. According to theory of cyclotomy,^[8] the probability of such an event for the case when ${\displaystyle \lambda _{1},\ldots ,\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},\ldots ,\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}$.

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

1. Due to the binomial theorem ${\textstyle (x-z)^{k}=\sum \limits _{i=0}^{k}{\binom {k}{i}}(-z)^{k-i}x^{i}}$, 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 ${\textstyle O(n^{2})}$, thus calculation of ${\textstyle x^{2^{k}}{\bmod {f}}_{z}(x)}$ is done in ${\textstyle 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,^[9] 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.^[7]