The fullfactor PARI/GP package

NOTE: this package requires a recent (development) version of PARI/GP. This is because it uses the new parser (especially t_CLOSUREs) in a fundamental way. It has been tested to work on pari-2.4.2.

Contents

  1. fullfactor()
  2. Factorfuncs
  3. Included factorfuncs
  4. Cyclotomic factorisations
  5. Other functions
  6. Download and Installation

fullfactor()

The main function of this package is fullfactor(x, {factorfunc}). We will get back later to the optional argument factorfunc. When called with one argument x, it does the following: if x is a matrix, it is assumed to be a partial factorisation of something. With a partial factorisation of z, we mean a factorisation of z, but where the factors are not necessarily irreducible. For example the output of factor(z,0) is a partial factorisation. If x is not a matrix, it is something to be factored. In this sense, fullfactor(x) and fullfactor(Mat([x,1])) are equivalent if x is not a matrix. In any case, fullfactor will return the full factorisation. It does this by calling factor(z) on every z in the partial factorisation. For example:

? F = factor(10^50 + 1, 0)
%1 =
[101 1]

[3541 1]

[27961 1]

[60101 1]

[166386582569341608294371141910949901 1]
? fullfactor(F)
%2 =
[101 1]

[3541 1]

[27961 1]

[60101 1]

[7019801 1]

[14103673319201 1]

[1680588011350901 1]

This can also come in handy if you want to factor the value of a reducible polynomial, since the polynomial factorisation automatically gives a factorisation of the value. This was the original motivation for me to write fullfactor. So the factorisation of 1050 + 1 can also be computed as follows:

? F = factor(x^50 + 1)
%3 =
[x^2 + 1 1]

[x^8 - x^6 + x^4 - x^2 + 1 1]

[x^40 - x^30 + x^20 - x^10 + 1 1]

Since x50 + 1 = (x100 - 1)/(x50 - 1), it's easy to check that the above factors should be the cyclotomic polynomials Φ4, Φ20 and Φ100. Now we can get a partial factorisation of 1050 + 1 and then use fullfactor:

? subst(F, x, 10)
%4 =
[101 1]

[99009901 1]

[9999999999000000000099999999990000000001 1]

? fullfactor(%)
%5 =
[101 1]

[3541 1]

[27961 1]

[60101 1]

[7019801 1]

[14103673319201 1]

[1680588011350901 1]

Factorfuncs

This already looks interesting, but what if we don't want to call factor(z) for every partial factor z, but something else? For example, we might want to factor polynomials over a number field using nffactor(nf,x). That's why fullfactor has a second argument factorfunc. It is a t_CLOSURE which says how every partial factor should be factored. Let's say we want to factor x20 + 1 over the field Q(√10). Since we know that x20 + 1 equals Φ40 times Φ8, we can do the following:

? F = [polcyclo(40),1; polcyclo(8),1]
%6 =
[x^16 - x^12 + x^8 - x^4 + 1 1]

[x^4 + 1 1]

? nf = nfinit(t^2 - 10);
? F2 = fullfactor(F, z->nffactor(nf, z))
%8 =
[x^4 + 1 1]

[x^8 + Mod(-t, t^2 - 10)*x^7 + 5*x^6 + Mod(-2*t, t^2 - 10)*x^5 + 7*x^4 + Mod(-2*t, t^2 - 10)*x^3 + 5*x^2 + Mod(-t, t^2 - 10)*x + 1 1]

[x^8 + Mod(t, t^2 - 10)*x^7 + 5*x^6 + Mod(2*t, t^2 - 10)*x^5 + 7*x^4 + Mod(2*t, t^2 - 10)*x^3 + 5*x^2 + Mod(t, t^2 - 10)*x + 1 1]

This factorisation of x20 + 1 actually gives a so-called Aurifeuillian factorisation of 1050 + 1. If we evaluate x20 + 1 at x = 100√10;, we get 1050 + 1. As can be seen from the factors in F2, this substitution applied on F2 will yield integers.

? lift(subst(F2, x, Mod(100*t,nf.pol)))
%9 =
[10000000001 1]

[99004980069800499001 1]

[101005020070200501001 1]

? fullfactor(%)
%10 =
[101 1]

[3541 1]

[27961 1]

[60101 1]

[7019801 1]

[14103673319201 1]

[1680588011350901 1]

To simplify doing several factorisations and substitutions, fullfactor also accepts a vector of factorfuncs which are applied from left to right. So the last example could also be written as follows:

? F = [polcyclo(40),1; polcyclo(8),1]; nf = nfinit(t^2 - 10);
? fullfactor(F, [z->nffactor(nf,z), z->lift(subst(z,x,100*Mod(t,nf.pol))), z->factor(z)])
%12 =
[101 1]

[3541 1]

[27961 1]

[60101 1]

[7019801 1]

[14103673319201 1]

[1680588011350901 1]

Looking at the second component of the factorfuncs vector, you can see that a factorfunc does not need to return a matrix. Returning x is then equivalent to returning Mat([x,1]).

Included factorfuncs

The fullfactor package already includes some useful factorfuncs, we will give an overview of them.

factor_default(x): call factor(x).

factor_coprime(x): this function does nothing and simply returns x. It is useful though, because fullfactor always returns a factorisation in coprimes. For example:

? F = matrix(6, 2, i, j, if(j==1, i+1, 1))
%13 =
[2 1]

[3 1]

[4 1]

[5 1]

[6 1]

[7 1]

? fullfactor(%, factor_coprime)
%14 =
[2 4]

[3 2]

[5 1]

[7 1]

factor_trialdiv(x): call factor(x,0).

factor_power(x): check whether x is a perfect power and return a factorisation of the form an.

factor_squarefree(x): square-free factorisation of the polynomial x. This means that x will be factored as x1 x22 x33xnn, where the xi's are mutually coprime.

factor_above_primelimit(x): this returns x if x is at least primelimit, otherwise it returns the empty factorisation. This way, fullfactor(x, [factor_default, factor_above_primelimit]) can be used to keep only those factors which would not be found by trial division. In other words, they are the factors which would have to be added with addprimes.

factor_check(x): this does not factor x, but checks whether the given factors are irreducible using isprime for integers or polisirreducible for polynomials. If all factors are indeed irreducible, it acts as factor_coprime, otherwise it will give an error. This can also be used to assure that all pseudoprimes in the output of factorint are actually proven primes (however, see also the factor_proven default). Example:

? fullfactor(10^50 + 1, [factor_default, factor_check])
%15 =
[101 1]

[3541 1]

[27961 1]

[60101 1]

[7019801 1]

[14103673319201 1]

[1680588011350901 1]

Now we know for sure that these factors are proven primes. It the check fails, you get the following:

? fullfactor(10^50 + 1, [factor_trialdiv, factor_check])
  ###   user error: The following factor is NOT prime:
166386582569341608294371141910949901

factor_gcd(v): the function factor_gcd is not a factorfunc by itself, but a function which returns a factorfunc (thanks to closures!). The argument v must be a vector or one element which is interpreted as [v]. Then factor_gcd(v)(x) computes gcd's of x with all elements from v, and returns the partial factorisation of x obtained in this manner. See also split_factor_gcd for a similar function which is meant to be used stand-alone (not with fullfactor).

One use of this is to simulate factoring with an addprimes-prime:

? fullfactor(10^50 + 1, [factor_trialdiv, factor_gcd(7019801)])
%16 =
[101 1]

[3541 1]

[27961 1]

[60101 1]

[7019801 1]

[23702464296258769770591950101 1]

A different example:

? fullfactor(2^7*3^8, factor_gcd(6))
%17 =
[2 7]

[3 8]

However,

? fullfactor(2^7*3^7, factor_gcd(6))
%18 =
[6 7]

Finally, remark that all these functions can also be used stand-alone, i.e. without fullfactor(). The main difference will be that the output will not be sorted, because the sorting happens in fullfactor(). For example, one can do

? pol = prod(i=1,5,(x^i-3)^i)*Mod(1,5); factor_squarefree(pol)
%19 =
[Mod(2, 5)*x^2 + Mod(4, 5) 2]

[Mod(1, 5)*x^3 + Mod(2, 5) 3]

[Mod(1, 5)*x^4 + Mod(2, 5) 4]

[Mod(4, 5)*x + Mod(3, 5) 26]

Compare this to

? pol = prod(i=1,5,(x^i-3)^i)*Mod(1,5); fullfactor(pol, factor_squarefree)
%20 =
[Mod(4, 5)*x + Mod(3, 5) 26]

[Mod(2, 5)*x^2 + Mod(4, 5) 2]

[Mod(1, 5)*x^3 + Mod(2, 5) 3]

[Mod(1, 5)*x^4 + Mod(2, 5) 4]

Cyclotomic factorisations

The factor_cyclo(x) factorfunc deserves a section on its own. This function computes cyclotomic factorsations, i.e. it tries to find factors in x which divide be - 1 for “small” b and e. By default, this means b < 100.

The algorithm consists essentially of two steps: first, an algorithm similar to the p - 1 method is used to find common factors between x and be - 1. The main difference with classical p - 1 factoring is that we use several well-chosen bases instead of one (random) base. The second part is then to factor be - 1 using (if possible) Aurifeuillian factorisation and table lookup. Note that the resulting factors are usually not primes, in general we only get a partial factorization.

Aurifeuillian factorisations are a special type of algebraic factorisations. Consider bn - 1 and write b = s2k with k squarefree. If k is 1 (mod 4) and n is k (mod 2k) or k is 2 or 3 (mod 4) and n is 2k (mod 4k), we have an Aurifeuillian factorisation. This comes from non-trivial polynomial factorisations of Φk(k x2) for k ≡ 1 (mod 4) and Φ2k(k x2) for k ≡ 2 or 3 (mod 4). You can compute these polynomial factorisations yourself using polaurif(k,x).

At present, the fullfactor package contains all tables from the Cunningham Project, see also my Cunningham page. These tables contain known primes dividing be ± 1 for b ≤ 12 and e up to a certain bound.

You can set default(debug) to 2, 3 or 4 to get some debug info and see what the algorithm does:

? \g3
   debug = 3
? factor_cyclo(10^50 + 1)
  *** factor: Warning: IFAC: untested integer declared prime.
        166386582569341608294371141910949901
factor_cyclo1: STAGE 1: B1 = 49
    Temp factor 166386582569341608294371141910949901 in 10^5d - 1
factor_cyclo1: Finished stage 1 in 0 ms
factor_cyclo1: STAGE 3: 1 factor to identify
    Processing factor 166386582569341608294371141910949901 in 10^5d - 1
        Sub factor 166386582569341608294371141910949901 in 10^25d - 1
    Processing factor 166386582569341608294371141910949901 in 10^25d - 1
        Sub factor 166386582569341608294371141910949901 in 10^50d - 1
    Processing factor 166386582569341608294371141910949901 in 10^50d - 1
        Sub factor 166386582569341608294371141910949901 in 10^100d - 1
gp_readvec_stream: reaching 16 entries
gp_readvec_stream: reaching 32 entries
gp_readvec_stream: reaching 64 entries
gp_readvec_stream: reaching 128 entries
gp_readvec_stream: found 145 entries
    Aurifeuillian factorisation 99004980069800499001 * 1680588011350901
 ** Final factor 7019801 in 10^100 - 1
 ** Final factor 14103673319201 in 10^100 - 1
 ** Final factor 1680588011350901 in 10^100 - 1
factor_cyclo1: Finished stage 3 in 4 ms
%21 =
[101 1]

[3541 1]

[27961 1]

[60101 1]

[7019801 1]

[14103673319201 1]

[1680588011350901 1]

As mentioned before, by default all bases up to 99 are tried by factor_cyclo. This can be changed however, for example if you are always working with one particular base. Use the function factor_cyclo_setbases(B) to do this. If B is an integer, all the bases 2, …, B will be selected. If B is a vector, the elements of B specify the bases. Note that bases which are a perfect power are always skipped. You can call factor_cyclo_setbases() without argument to retrieve the current set of bases.

Other functions

We describe some other, more specialised, functions in the fullfactor package.

We start with radical(x,{factorfunc}). This computes the radical of x, the product of all the prime factors in x; equivalently, the largest squarefree divisor of x. If x is a matrix, it is assumed to be a factorisation matrix. If x is not a matrix, it will be factored (using factorfunc if specified). If factorfunc is specified, fullfactor(x,factorfunc) will always be called, whether or not x is a matrix.

If you want to factor be - 1, you could do factor_cyclo(b^e - 1). It's more efficient however to use the function partial_factor_pn_1(b,e), which doesn't need to figure out the b and e. As the name suggests, this is only a partial factorisation using the available data. If you want to do a full factorisation, use fullfactor(partial_factor_pn_1(b,e)). The latter should be functionally identical to the libPARI call factor_pn_1(b,e) if b is prime.

The factorfunc factor_cyclo which we mentioned before is actually only a wrapper for factor_cyclo1(x,{flag=0}). Indeed, factor_cyclo(x) is defined as fullfactor(x, [factor_trialdiv, factor_cyclo1]). So it first does trial division. This is necessary because the tables only store factors which are at least 500000 (the default primelimit). If you know for some reason that your number does not contain any small factors, you could use factor_cyclo1 directly instead of factor_cyclo. The function factor_cyclo1(x,{flag=0}) has an optional argument flag. This extra argument only makes sense when factor_cyclo1 is used stand-alone, not as a factorfunc for fullfactor. If flag equals 1, then factor_cyclo1 returns a 4-column factorisation matrix F, where the first two columns form the usual factorisation matrix. The third and fourth columns are such that F[i,1] divides polcyclo(F[i,4],F[i,3]). If no such cyclotomic factor can be found, then F[i,3] and F[i,4] will both be zero.

split_factor_gcd(x,v) is the implementation for the factor_gcd factorfunc. The function split_factor_gcd returns a two-element vector [F,c], where c is the maximal divisor of x which is coprime to all elements of v. Then F is a partial factorisation of x/c, obtained by taking gcd's between x and elements of v.

split_factor_cyclo(x,b,e) is the most low-level function for cyclotomic factoring, called by both factor_cyclo1 and partial_factor_pn_1. Similarly to split_factor_gcd, it returns a vector [F,c]. In this case, c is the maximal divisor of x which is coprime to be - 1. The matrix F is a factorisation matrix of x/c obtained by Aurifeuillian factorisation of Φe(b) and table lookup of the factors of Φe(b). Note that factors of Φd(b) with d dividing e are not separated. Therefore x should be a divisor of Φe(b). Also factors smaller than 500000 are not stored in the tables, so you should do trial division yourself.

Download and Installation

There are two tarballs to be downloaded:

To load the package in your GP session, you just need to type \r fullfactor.gp. When you do this, the location of the data files will be automatically determined. A warning will be given if they can't be found. This location is stored in the variable fullfactor_datadir, which you can set yourself if necessary. If you don't want to use the data, you can set fullfactor_datadir="". This way, only algebraic (cyclotomic and Aurifeuillian) factorisations will be found by factor_cyclo and related functions.

Happy fullfactoring!

Jeroen Demeyer
jdemeyer@cage.ugent.be Last update: 2012-01-26