Английская Википедия: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 <math>\mathbb F_p</math> with <math>p</math> 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]
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[4] for finding square roots in <math>\mathbb F_p</math>.[5] In 2000 Peralta's method was generalized for cubic equations.[6]
Statement of problem
Let <math>p</math> be an odd prime number. Consider the polynomial <math display="inline">f(x) = a_0 + a_1 x + \cdots + a_n x^n</math> over the field <math>\mathbb F_p\simeq \mathbb Z/p\mathbb Z</math> of remainders modulo <math>p</math>. The algorithm should find all <math>\lambda</math> in <math>\mathbb F_p</math> such that <math display="inline">f(\lambda)= 0</math> in <math>\mathbb F_p</math>.[2][7]
Algorithm
Randomization
Let <math display="inline">f(x) = (x-\lambda_1)(x-\lambda_2)\cdots(x-\lambda_n)</math>. 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 <math display="inline">f_z(x)=f(x-z) = (x-\lambda_1 - z)(x-\lambda_2 - z) \cdots (x-\lambda_n-z)</math> where <math>z</math> is some any element of <math>\mathbb F_p</math>. If one can represent this polynomial as the product <math>f_z(x)=p_0(x)p_1(x)</math> then in terms of the initial polynomial it means that <math>f(x) =p_0(x+z)p_1(x+z)</math>, which provides needed factorization of <math>f(x)</math>.[1][7]
Classification of <math>\mathbb F_p</math> elements
Due to Euler's criterion, for every monomial <math>(x-\lambda)</math> exactly one of following properties holds:[1]
- The monomial is equal to <math>x</math> if <math>\lambda = 0</math>,
- The monomial divides <math display="inline">g_0(x)=(x^{(p-1)/2}-1)</math> if <math>\lambda</math> is quadratic residue modulo <math>p</math>,
- The monomial divides <math display="inline">g_1(x)=(x^{(p-1)/2}+1)</math> if <math>\lambda</math> is quadratic non-residual modulo <math>p</math>.
Thus if <math>f_z(x)</math> is not divisible by <math>x</math>, which may be checked separately, then <math>f_z(x)</math> is equal to the product of greatest common divisors <math>\gcd(f_z(x);g_0(x))</math> and <math>\gcd(f_z(x);g_1(x))</math>.[7]
Berlekamp's method
The property above leads to the following algorithm:[1]
- Explicitly calculate coefficients of <math>f_z(x) = f(x-z)</math>,
- Calculate remainders of <math display="inline">x,x^2, x^{2^2},x^{2^3}, x^{2^4}, \ldots, x^{2^{\lfloor \log_2 p \rfloor}}</math> modulo <math>f_z(x)</math> by squaring the current polynomial and taking remainder modulo <math>f_z(x)</math>,
- Using exponentiation by squaring and polynomials calculated on the previous steps calculate the remainder of <math display="inline">x^{(p-1)/2}</math> modulo <math display="inline">f_z(x)</math>,
- If <math display="inline">x^{(p-1)/2} \not \equiv \pm 1 \pmod{f_z(x)}</math> then <math>\gcd</math> mentioned above provide a non-trivial factorization of <math>f_z(x)</math>,
- Otherwise all roots of <math>f_z(x)</math> are either residues or non-residues simultaneously and one has to choose another <math>z</math>.
If <math>f(x)</math> is divisible by some non-linear primitive polynomial <math>g(x)</math> over <math>\mathbb F_p</math> then when calculating <math>\gcd</math> with <math>g_0(x)</math> and <math>g_1(x)</math> one will obtain a non-trivial factorization of <math>f_z(x)/g_z(x)</math>, thus algorithm allows to find all roots of arbitrary polynomials over <math>\mathbb F_p</math>.
Modular square root
Consider equation <math display="inline">x^2 \equiv a \pmod{p}</math> having elements <math>\beta</math> and <math>-\beta</math> as its roots. Solution of this equation is equivalent to factorization of polynomial <math display="inline">f(x) = x^2-a=(x-\beta)(x+\beta)</math> over <math>\mathbb F_p</math>. In this particular case problem it is sufficient to calculate only <math>\gcd(f_z(x); g_0(x))</math>. For this polynomial exactly one of the following properties will hold:
- GCD is equal to <math>1</math> which means that <math>z+\beta</math> and <math>z-\beta</math> are both quadratic non-residues,
- GCD is equal to <math>f_z(x)</math>which means that both numbers are quadratic residues,
- GCD is equal to <math>(x-t)</math>which means that exactly one of these numbers is quadratic residue.
In the third case GCD is equal to either <math>(x-z-\beta)</math> or <math>(x-z+\beta)</math>. It allows to write the solution as <math display="inline">\beta = (t - z) \pmod{p}</math>.[1]
Example
Assume we need to solve the equation <math display="inline">x^2 \equiv 5\pmod{11}</math>. For this we need to factorize <math>f(x)=x^2-5=(x-\beta)(x+\beta)</math>. Consider some possible values of <math>z</math>:
- Let <math>z=3</math>. Then <math>f_z(x) = (x-3)^2 - 5 = x^2 - 6x + 4</math>, thus <math>\gcd(x^2 - 6x + 4 ; x^5 - 1) = 1</math>. Both numbers <math>3 \pm \beta</math> are quadratic non-residues, so we need to take some other <math>z</math>.
- Let <math>z=2</math>. Then <math>f_z(x) = (x-2)^2 - 5 = x^2 - 4x - 1</math>, thus <math display="inline">\gcd( x^2 - 4x - 1 ; x^5 - 1)\equiv x - 9 \pmod{11}</math>. From this follows <math display="inline">x - 9 = x - 2 - \beta</math>, so <math>\beta \equiv 7 \pmod{11}</math> and <math display="inline">-\beta \equiv -7 \equiv 4 \pmod{11}</math>.
A manual check shows that, indeed, <math display="inline">7^2 \equiv 49 \equiv 5\pmod{11}</math> and <math display="inline">4^2\equiv 16 \equiv 5\pmod{11}</math>.
Correctness proof
The algorithm finds factorization of <math>f_z(x)</math> in all cases except for ones when all numbers <math>z+\lambda_1, z+\lambda_2, \ldots, z+\lambda_n</math> are quadratic residues or non-residues simultaneously. According to theory of cyclotomy,[8] the probability of such an event for the case when <math>\lambda_1, \ldots, \lambda_n</math> are all residues or non-residues simultaneously (that is, when <math>z=0</math> would fail) may be estimated as <math>2^{-k}</math> where <math>k</math> is the number of distinct values in <math>\lambda_1, \ldots, \lambda_n</math>.[1] In this way even for the worst case of <math>k=1</math> and <math>f(x)=(x-\lambda)^n</math>, the probability of error may be estimated as <math>1/2</math> and for modular square root case error probability is at most <math>1/4</math>.
Complexity
Let a polynomial have degree <math>n</math>. We derive the algorithm's complexity as follows:
- Due to the binomial theorem <math display="inline">(x-z)^k = \sum\limits_{i=0}^k \binom{k}{i} (-z)^{k-i}x^i</math>, we may transition from <math>f(x)</math> to <math>f(x-z)</math> in <math>O(n^2)</math> time.
- Polynomial multiplication and taking remainder of one polynomial modulo another one may be done in <math display="inline">O(n^2)</math>, thus calculation of <math display="inline">x^{2^k} \bmod f_z(x)</math> is done in <math display="inline">O(n^2 \log p)</math>.
- Binary exponentiation works in <math>O(n^2 \log p)</math>.
- Taking the <math>\gcd</math> of two polynomials via Euclidean algorithm works in <math>O(n^2)</math>.
Thus the whole procedure may be done in <math>O(n^2 \log p)</math>. Using the fast Fourier transform and Half-GCD algorithm,[9] the algorithm's complexity may be improved to <math>O(n \log n \log pn)</math>. For the modular square root case, the degree is <math>n = 2</math>, thus the whole complexity of algorithm in such case is bounded by <math>O(\log p)</math> per iteration.[7]
References
Шаблон:Number-theoretic algorithms
