Loading [MathJax]/jax/output/HTML-CSS/jax.js

Monday, 6 August 2012

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=2y, 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

Thursday, 19 July 2012

Irreducibility testing

Fact. A polynomial fFq[x] of degree n1 is irreducible if and only if
1) xqnxmodf and
2) gcd(xqn/tx,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 fFq[x] of degree n. On output test gives an answer "reducible" or "irreducible".
Step 1.  Compute a=xqnremf. If ax return "reducible".
Step 2. Find prime divisors of n.
Step 3. For all prime divisors t of n compute b=xqn/tremf. If gcd(bx,f)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 degf[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.

Saturday, 7 July 2012

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 mN define the mth trace polynomial over F2 by Tm=x2m1+x2m2+...+x2+x.  Let q=2k for some kN, fFq[x] squarefree of degree n, with r2 irreducible factors f1,...,frFq[x],R=Fq[x]/<f>

This modification exploits two facts:

Fact 1. x2m+x=Tm(Tm+1)Tm(α)F2αF2m
Fact 2. Let all irreducible factors of f have degree d. Then Tkd(α)modfiF2iαR

Now one can modify a probabilistic algorithm of equal-degree splitting: instead of computing b=a(qd1)/2remf one have to compute b=Tkd(a)remf.

_____________________

Proof of Fact 1.

1) First note that (x1+...+xk)2n=x2n1+...+x2nkxiF2m. It can be proved by induction on n. For n=0 and n=1 it is obvious. Assuming that this is true for n1 one can get: (x1+...+xk)2n=((x1+...+xk)2n1)2=(x2n11+...+x2n1k)2=x2n1+...+x2nk

2) Then note that aq=aaFq

3) Finally,    Tm(α)(Tm(α)+1)=(Tm(α))2+Tm(α)=
                    (α2m1+α2m1...+α)2+α2m1+...+α2+α=
                    α2m+α2m1...+α2+α2m1+...+α2+α=
                    α2m+α=α+α=0αF2m
It means that Tm(α)=0 or Tm(α)=1, i.e. Tm(α)F2αF2m


Monday, 2 July 2012

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=Fq[x]/<f> is a vector space of dimension n=degf over Fq. Let β:RR be a mapping: β(a)=aqaConsider the kernel of βIf f = f1 ... fr where fi is irreducible i, then for aFq[x] we have
amodfkerβaqamodfaqamodfii
(by Fermat's little theorem)


Thus kerβ=Fq×...×Fq=Frq. amodfkerβ(amodf1,...,amodfr)=(a1modf1,...,armodfr) for some aiFq Denote by Q the matrix representing the Frobenius map S:S(a)=aq with respect to the polynomial basis xn1modf,...,xmodf,1modf of R.


Fact. f is irreducible r=1rank(QI)=n1

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 fFq[x] of degree n
On output we get an irreducible divisor of f or "failure".
Step 1. Compute xqiremf for i=0,...,n1
Step 2. Form matrix Q using the result of Step 1
Step 3. Compute rank of (QI) and basis b1,...,br of kerβ (e.g. using Gaussian elimination for QI). If rank=n1 return f.
Step 4. Choose arbitrary element a=c1b1+...+crbr,ciFq
Step 5. g1=gcd(a,f). If g11 return g1.
Step 6. g2=gcd(a(q1)/2remf1,f). If g21,g2f return g2, 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
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. 


Sunday, 24 June 2012

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 xqixqjFq[x] is divisible by precisely those irreducible polynomials in Fq[x] whose degree divides ij.
The proof is easy: if ij then  xqixqj=(xqijx)qj. From the Theorem from the previous post we know that the factorisation of xqkx 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 fFq[x] of degree n
The output is f[1],...,f[n]Fq[x] such that for 1dn, f[d] is the product of the monic irreducible factors of f of degree d. The algorithm is parameterized by a constant β, with 0β1.

Step 1 (compute baby steps)
Let l=nβ . For 0il, compute hi=xqimodf.
Step 2 (compute giant steps) 
Let m=n2l . For 1jm, compute Hj=xqljmodf
Step 3 (compute interval polynomials) 
For 1jm, compute Ij=0il(Hjhi)modf .
(by Lemma the polynomial Ij is divisible by those irreducible factors of f whose degree divides an integer k with  (j1)l<kjl)
Step 4 (compute coarse DDF) 
In this step, we compute polynomials F1,...,Fm where Fj=f[(j1)l+1]f[(j1)l+2]...f[jl] This is done as follows.
sf;
for j1 to m do
{Fjgcd(s,Ij);ssFj}
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 j1 to m do
{
gFj;
for il1 down to 0 do
    {f[lji]gcd(g,Hjhi);ggf[lji]}
}
if s1 then f[deg(s)]s

Some technical remarks
1.  Lemma. For any positive integer r if we are given h=xqrmodfFq[x], then for any gFq[x], we can compute  gqrmodf as g(h)modf.

Now to perform Step 2 we set H1=hl and then compute each Hj as Hj1(H1)modf.

2. Computing f(g)modh 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.

Monday, 18 June 2012

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 Zq (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 gi (a product of all squarefree irreducible factors of P of degree i), then it decomposes gi into equal degree factors fij (each of degree i) and for all j removes the highest possible degree of fij (say eij) from P simply dividing P by fij. The cycle continues for P=P (P is P without feijij) while degP>1.

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


Theorem. For any d1(xqdx)Fq[x] is the product of all monic irreducible polynomials in Fq[x] whose degree divides d.


Now we can easily construct an algorithm: for d=1,2,... compute h=xqdx and find g=gcd(h,f), continue this procedure for f=fg 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=Pf.

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(q1)/2{1,1}aFq (Fq  is a multiplicative group of  Fq). 


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 aFq[x] with dega<n at random. If dega=0 return "failure". Else compute g1=gcd(a,f). If degg1>0 return g1. Else (if degg1=0) compute b=a(qd1)/2remf (see Lemma) and g2=gcd(b1,f). If g21,g2f return g2, 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).

Tuesday, 5 June 2012

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.