Profiling

Sometimes it is very helpful to have a function which makes an automatic choice between available factorisation algorithms depending on an input polynomial. For this purpose I wrote some profiling code in nmod_poly and fmpz_mod_poly_factor. The results are as follows (along Oy axis: $y, n=2^y$, where $n$ is a degree of input polynomial, along Ox: $x=bits(q)$ is number of bits in modulo, the letters denote the boarder where one algorithm becomes faster than another: B - Berlekamp, CZ - Cantor-Zassenhaus, KS - Kaltofen-Shoup).

nmod_poly:

fmpz_mod_poly:

 I tried to make a simple formula which fits these data. Finally, I got:
nmod_poly:
if bits < 5 and n > 128, then use B,
else if bits >= 5 and 2 * bits + n > 74, then use KS
else use CZ

fmpz_mod_poly:
if 5 * bits + n > 75, then use KS, else use CZ

Irreducibility testing

Fact. A polynomial $f\in F_q[x]$ of degree $n\geq 1$ is irreducible if and only if
1) $x^{q^n}\equiv x \mod f$ and
2) $\gcd (x^{q^{n/t}}-x, f)=1$ for all prime divisors $t$ of $n$.

From this fact one can get the Rabin's irreducibility test:

Irreducibility test (Rabin)
Input and Output  On input we have $f\in F_q[x]$ of degree $n$. On output test gives an answer "reducible" or "irreducible".
Step 1.  Compute $a=x^{q^n} \mathrm{rem} f$. If $a\neq x$ return "reducible".
Step 2. Find prime divisors of $n$.
Step 3. For all prime divisors $t$ of $n$ compute $b=x^{q^{n/t}} \mathrm{rem} f$. If $\gcd(b-x, f)\neq 1$ return "reducible".
Step 4. Return "irreducible".

Another possibility is to use distinct-degree factorisation to test irreducibility.

Irreducubility test (ddf)
Step 1.  Check if $f$ is squarefree (by computing it's derivative). If not, return "reducible".
Step 2. Build ddf $f^{[1]},...,f^{[n]}$ for $f$.
Step 3. If $\deg f^{[n]}=n$ return "irreducible", else return "reducible".

Since fmpz_mod_poly doesn't have fast gcd computation and fast modular composition implemented, the ddf-based algorithm turns out to be faster for it.

Making factor_equal_deg_prob to work for q=2

In my previous posts I described the variant of equal-degree splitting for odd prime powers. This algorithm requires some modification for fields with characteristic 2.

For $m\in \mathbb{N}$ define the $m$th trace polynomial over $F_2$ by $T_m=x^{2^{m-1}}+x^{2^{m-2}}+...+x^2+x.$  Let $q = 2^k$ for some $k\in \mathbb{N}$, $f\in F_q[x]$ squarefree of degree $n$, with $r\geq 2$ irreducible factors $f_1,...,f_r\in F_q[x], R=F_q[x]/<f>$

This modification exploits two facts:

Fact 1. $x^{2^m}+x = T_m(T_m+1) \Rightarrow T_m(\alpha)\in F_2 \,\,\forall \alpha\in F_{2^m}$
Fact 2. Let all irreducible factors of $f$ have degree $d$. Then $T_{kd}(\alpha) \mod f_i\in F_2 \,\,\forall i \,\,\forall \alpha\in R$

Now one can modify a probabilistic algorithm of equal-degree splitting: instead of computing $b=a^{(q^d-1)/2}\mathrm{rem} f$ one have to compute $b = T_{kd}(a)\mathrm{rem} f$.

_____________________

Proof of Fact 1.

1) First note that $(x_1+...+x_k)^{2^n} = x_1^{2^n}+...+x_k^{2^n} \,\, \forall x_i\in F_{2^m}$. It can be proved by induction on $n$. For $n=0$ and $n=1$ it is obvious. Assuming that this is true for $n-1$ one can get: $(x_1+...+x_k)^{2^n}=\left((x_1+...+x_k)^{2^{n-1}}\right)^2 = (x_1^{2^{n-1}}+...+x_k^{2^{n-1}})^2=x_1^{2^n}+...+x_k^{2^n}$

2) Then note that $a^q = a \,\, \forall a\in F_q$

3) Finally,    $T_m(\alpha)(T_m(\alpha)+1)=(T_m(\alpha))^2+T_m(\alpha) = $
                    $(\alpha^{2^{m-1}}+\alpha^{2^{m-1}}...+\alpha)^2+\alpha^{2^{m-1}}+...+\alpha^2+\alpha=$
                    $\alpha^{2^m}+\alpha^{2^{m-1}}...+\alpha^2+\alpha^{2^{m-1}}+...+\alpha^2+\alpha=$
                    $\alpha^{2^m}+\alpha=\alpha+\alpha=0 \,\, \forall \alpha\in F_{2^m}$
It means that $T_m(\alpha) = 0$ or $T_m(\alpha) = 1$, i.e. $T_m(\alpha)\in F_2\,\,\forall \alpha\in F_{2^m}$

Berlekamp

Yesterday I finished porting Berlekamp's algorithm from nmod_poly module to fmpz_mod_poly_factor module.

Unlike Cantor-Zassenhaus and Kaltofen-Shoup algorithms this algorithm uses linear algebra -- not number theory. Berlekamp algorithm consists of two parts. The first part is standard: squarefree factorization. The second part is more specific.

Assume for simplicity that q is an odd prime number. $R = F_q[x]/<f>$ is a vector space of dimension $n = \deg f$ over $F_q$. Let $\beta : R \to R$ be a mapping: $\beta(a)  = a^q-a$. Consider the kernel of $\beta$. If $f$ = $f_1$ ... $f_r$ where $f_i$ is irreducible $\forall i$, then for $a \in F_q[x]$ we have
$$a \mod f \in \ker \beta \Leftrightarrow a^q \equiv a \mod f \Leftrightarrow a^q \equiv a \mod f_i \,\,\, \forall i$$
(by Fermat's little theorem)


Thus $\ker \beta = F_q \times ... \times F_q = F_q^r$. $a \mod f \in \ker \beta \Leftrightarrow (a \mod f_1, ... , a \mod f_r) = (a_1 \mod f_1, ... , a_r \mod f_r)$ for some $a_i \in F_q$ Denote by $Q$ the matrix representing the Frobenius map $S: S(a) = a^q$ with respect to the polynomial basis $x^{n-1} \mod f, ... , x \mod f, 1 \mod f$ of $R$.


Fact. $f$ is irreducible $\Leftrightarrow r = 1 \Leftrightarrow \mathrm{rank}(Q-I) = n-1$

Now we can formulate the second part of the algorithm. 

Berlekamp's algorithm (second part).
Input and Output.
On input we have a monic squarefree polynomial $f\in F_q[x]$ of degree $n$. 
On output we get an irreducible divisor of $f$ or "failure".
Step 1. Compute $x^{qi} \mathrm{rem} f$ for $i=0,...,n-1$
Step 2. Form matrix $Q$ using the result of Step 1
Step 3. Compute rank of $(Q-I)$ and basis $b_1, ... , b_r$ of $\ker \beta$ (e.g. using Gaussian elimination for $Q-I$). If $\mathrm{rank} = n-1$ return $f$.
Step 4. Choose arbitrary element $a = c_1 b_1 + ... + c_r b_r, c_i \in F_q$
Step 5. $g_1 = \gcd(a,f)$. If $g_1 \neq 1$ return $g_1$.
Step 6. $g_2 = \gcd(a^{(q-1)/2} \mathrm{rem} f - 1, f)$. If $g_2 \neq 1, g_2 \neq f$ return $g_2$, else return "failure"


To get the full variant of factorization algorithm one has to call the second part many times until f is irreducible.

I also made Kaltofen-Shoup algorithm to use Brent-Kung algorithm for modular composition. I measured working time of Cantor-Zassenhaus, Kaltofen-Shoup and Berlekamp algorithms in the simplest case when an input polynomial has 1..6 factors of degree 1..8 and multiplicity 1..31. On each iteration (one iteration is a complete factorization of one polynomial) all parameters of input polynomial are randomly chosen (with uniform distribution). There is a table describing results on my computer (NI $-$ number of iterations, CZ $-$ Cantor-Zassenhaus, B $-$ Berlekamp, KS $-$ Kaltofen-Shoup):

NI                  :    100      200     300     400     500     600     700     800     900    1000

CZ time (sec) :   10.71   18.98   27.20  35.53  46.02  55.36   64.10  72.11  82.54  91.79

B time (sec)   :    1.65     2.98    4.34    6.65    8.36    9.79    11.44  12.92  14.53  15.91

KS time (sec) :    1.55     2.92    4.09    5.80    7.23    8.56    10.06  11.56  12.91  14.13

Of course it's necessary to make more tests because algorithm effectiveness can depend on some special properties of input polynomials or on field characteristic. 

Baby/giant step strategy

I wrote in my previous post that Cantor and Zassenhaus' algorithm of polynomial factorisation over finite fields can be divided into three stages:

• squarefree factorisation
• distinct-degree factorisation
• equal-degree factorisation

(in the standard Cantor-Zassenhaus algorithm the squarefree factorisation step is merged with the algorithm itself)

According to Modern Computer Algebra, "the second stage consumes the bulk of the computing time". In 1998 Erich Kaltofen and Victor Shoup designed a new approach to distinct-degree factorisation which allows to decrease it's cost. This strategy is called "baby/giant step" and exploits the fact from number theory:

Lemma. For nonnegative integers $i$ and $j$, the polynomial $x^{q^i} - x^{q^j} \in F_q [x]$ is divisible by precisely those irreducible polynomials in $F_q [x]$ whose degree divides $i−j$.

The proof is easy: if $i \geq j$ then  $x^{q^i} - x^{q^j} = (x^{q^{i-j}} - x)^{q^j}$. From the Theorem from the previous post we know that the factorisation of $x^{q^k} - x$ consists of all irreducible factors whose degree is a divisor of $k$. 

Now it's possible to formulate the algorithm:

Fast distinct-degree factorisation (DDF)

Input and output

This algorithm takes as input a square-free polynomial $f \in F_q [x]$ of degree $n$. 

The output is $f ^{[1]} , . . . , f^{[n]} \in F_q [x]$ such that for $1 \leq d \leq n$, $f ^{[d]}$ is the product of the monic irreducible factors of f of degree $d$. The algorithm is parameterized by a constant $\beta$, with $0 \leq \beta \leq 1$.

Step 1 (compute baby steps)

Let $l = \lceil n\beta \rceil$ . For $0\leq i\leq l$, compute $h_i = x^{q^i} \mod f$.
Step 2 (compute giant steps) 

Let $m = \lceil \frac{n}{2l} \rceil$ . For $1\leq j\leq m$, compute $H_j = x^{q^{lj}} \mod f$. 

Step 3 (compute interval polynomials) 

For $1\leq j\leq m$, compute $I_j = \prod\limits_{0\leq i\leq l}(H_j - h_i) \mod f$ .

(by Lemma the polynomial $I_j$ is divisible by those irreducible factors of $f$ whose degree divides an integer $k$ with  $(j-1)l < k \leq jl$)

Step 4 (compute coarse DDF) 

In this step, we compute polynomials $F_1, ... , F_m $ where $F_j = f^{[(j-1)l+1]} f^{[(j-1)l+2]} ... f^{[jl]}$.  This is done as follows.
$s\leftarrow f;$
for $j \leftarrow 1$ to $m$ do
$\left\{F_j \leftarrow \gcd(s, I_j); s\leftarrow \frac{s}{F_j}\right\}$
Step 5 (compute fine DDF) 

In this step, we compute the output polynomials $f^{[1]}, ... , f ^{[n]}$. 

First, initialize $f^{[1]}, ... , f^{[n]}$ to 1. Then do the following.

for $j \leftarrow1$ to $m$ do

{

$g\leftarrow F_j ;$

for $i\leftarrow l - 1$ down to 0 do
    $\left\{ f^{[lj-i]} \leftarrow \gcd(g, H_j - h_i); g\leftarrow \frac{g}{f^{[lj-i]}} \right\}$

}

if $s \neq 1$ then $f^{[\deg(s)]}\leftarrow s$

Some technical remarks

1.  Lemma. For any positive integer $r$ if we are given $h = x^{q^r} \mod f \in F_q [x]$, then for any $g \in F_q [x]$, we can compute  $g^{q^r} \mod f$ as $g(h) \mod f$.

Now to perform Step 2 we set $H_1 = h_l$ and then compute each $H_j$ as $H_{j−1} (H_1) \mod f$.

2. Computing $f(g) \mod h$ is a so-called "modular composition" problem. For solving this problem one can use e.g. Horner scheme or Brent and Kung's algorithm (this also requires fast matrix multiplication).

By now I implemented the simplest version of the fast CZ algorithm: it uses Horner scheme in Step 2.

The next step might be to implement matrix multiplication and Brent and Kung's algorithm.

Cantor-Zassenhaus algorithm

Last Saturday I finished the first big part of my project: I ported Cantor-Zassenhaus algorithm with all helper functions from nmod_poly module to fmpz_mod_poly module.

I think it was an important part of my job because it helped me to familiarise myself with the code: I found some important differences between nmod and fmpz modules (I spent a lot of time debugging my stupid errors relates to these differences).

During the last month I got a good impression of the code of FLINT. I like it's clearness: every function has it's own file and detailed description, the function names are easy-to-understand (I mean that one can guess what the function should do according to it's name), pieces of code doing one distinct thing are organised in separate functions (so the complex algorithms don't look ugly).


Now I'm going to describe Cantor-Zassenhaus algorithm (the complete description can be found in Modern Computer Algebra book).

Input and output
The algorithm gets on input a nonconstant polynomial $P$ over a field $Z_q$ ($q$ is an odd prime number). On output it gives all monic irreducible factors of P with their multipliсities.

Short description
In brief the algorithm can be described as follows: in a cycle for $i=1,2,...$ it computes distinct-degree polynomial $g_i$ (a product of all squarefree irreducible factors of $P$ of degree $i$), then it decomposes $g_i$ into equal degree factors $f_{ij}$ (each of degree $i$) and for all $j$ removes the highest possible degree of $f_{ij}$ (say $e_{ij}$) from $P$ simply dividing $P$ by $f_{ij}$. The cycle continues for $P=P'$ ($P'$ is $P$ without $f_{ij}^{e_{ij}}$) while $\deg P' > 1$.

Distinct-degree factorisation
This algorithm is based on Fermat's little theorem, more precisely on the consequence of the following


Theorem. For any $d\geq 1 (x^{q^d}-x) \in F_q[x]$ is the product of all monic irreducible polynomials in $F_q[x]$ whose degree divides $d$.


Now we can easily construct an algorithm: for $d = 1, 2, ...$ compute $h = x^{q^d} - x$ and find $g=\gcd(h,f)$, continue this procedure for $f=\frac{f}{g}$ until $f=1$.


Equal-degree factorisation
It is a recursive algorithm: it calls a probabilistic procedure of finding a factor $g$ of a given degree of a polynomial $P$, then it removes the found factor $g$ from $P$ and calls itself for $P=\frac{P}{f}$.

Equal-degree splitting
It is a probabilistic procedure of finding a factor $g$ of a given degree of a polynomial $P$. It is based on the following 


Lemma. Let $q$ be an odd prime number. Then $a^{(q-1)/2}\in\{1,-1\} \,\,\, \forall a\in F_q^{*}$ ($F_q^{*}$  is a multiplicative group of  $F_q$). 


Let the polynomial $P$ have the degree $n=rd$ and every it's single factor have the degree $d$. Then the probabilistic algorithm is as follows: choose $a\in F_q[x]$ with $\deg a < n$ at random. If $\deg a = 0$ return "failure". Else compute $g_1=\gcd(a,f)$. If $\deg g_1 > 0$ return $g_1$. Else (if $\deg g_1 = 0$) compute $b=a^{(q^d-1)/2}\mathrm{rem} f$ (see Lemma) and $g_2=\gcd(b-1, f)$. If $g_2\neq 1, g_2\neq f$ return $g_2$, else return "failure".


Some technical remarks
To compute powers of q the binary exponentiation method can be applied.
To compute the gcd  euclidean method can be used.


According to my project proposal I'm out of schedule, so I would like to change the order of the next two algorithms and to continue now with "baby/giant step" algorithm (and to implement Berlekamp after it if there is time left).

Adding helper functions, part 1

Now I'm working on the second step of porting code for Cantor-Zassenhaus factorisation. It is necessary to move functions for polynomial arithmetic such as mulmod, remove, pow_trunc, pow_trunc_binexp, powmod_ui_binexp. For now I moved the first two functions.