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

## Contents

- fullfactor()
- Factorfuncs
- Included factorfuncs
- Cyclotomic factorisations
- Other functions
- 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 10^{50} + 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 *x*^{50} + 1 = (*x*^{100} - 1)/(*x*^{50} - 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 10^{50} + 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 *x*^{20} + 1 over the field **Q**(√10).
Since we know that *x*^{20} + 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 *x*^{20} + 1 actually gives a so-called
Aurifeuillian factorisation of 10^{50} + 1.
If we evaluate *x*^{20} + 1 at *x* = 100√10;,
we get 10^{50} + 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 *a*^{n}.

factor_squarefree(*x*):
square-free factorisation of the polynomial *x*.
This means that *x* will be factored as
*x*_{1} *x*_{2}^{2} *x*_{3}^{3}… *x*_{n}^{n},
where the *x*_{i}'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 *b*^{e} - 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 *b*^{e} - 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 *b*^{e} - 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 *b*^{n} - 1
and write *b* = *s*^{2}*k* with *k* squarefree.
If *k* is 1 (mod 4) and *n* is *k* (mod 2*k*)
or *k* is 2 or 3 (mod 4) and *n* is *2k* (mod 4*k*),
we have an Aurifeuillian factorisation.
This comes from non-trivial polynomial factorisations of
Φ_{k}(*k x*^{2}) for *k* ≡ 1 (mod 4)
and Φ_{2k}(*k x*^{2}) 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 *b*^{e} ± 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 *b*^{e} - 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 *b*^{e} - 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:

- fullfactor.tgz (2012-01-26): the relevant GP scripts.
Put these files somewhere in GP's
`default(path)`. - fullfactor-data.tgz (always updated):
data files containing factors of
`polcyclo(e,b)`. Extract these files either in GP's`default(datadir)`or`default(path)`(be sure to keep the`data/fullfactor`directory structure).

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!