(Part of this.)

We now have our squarefree polynomial $A_i$, which is the product of all the factors of $A$ with exponent $i$. Our goal is to factor it into $A_i = A_{i,1}A_{i,2}\dots A_{i,\ell}$, where $A_{i,d}$ is the product of all the factors of $A_i$ of degree $d$. This is much simpler than the squarefree factorisation algorithm; it is essentially based on the fact that, over a field $\mathbf{F}_p$, the irreducible factors of $X^{p^n}-X$ are precisely the (monic) irreducible polynomials whose degree is a divisor of $n$.

Theory

To ease notation, let $q = p^n$. We first show that all the monic irreducible polynomials of degree $n$ are factors of $X^q-X$. So as to not lengthen our exposition, we will use without proof the fact that there is a field with $q$ elements, say $K$. Now, the multiplicative group $K^\times$ has order $q-1$, meaning that for every element $a \in K^\times$ we have $a^{q-1} = 1$, and so $a^q = a$. Since also $0^q = 0$, we see that all the elements of $K$ are roots of $X^q-X$, and since it cannot have more than $q$ roots, they are all its roots and we have the factorisation
\[ X^q-X = \prod_{a \in K} (X-a). \]

Let now $P \in \mathbf{F}_p[X]$ be monic and irreducible of degree $n$. We know that there exists an extension of $\mathbf{F}_p$ in which $P$ has a root, say $\alpha$. The field $\mathbf{F}_p(\alpha)$ has $q$ elements, and since $\alpha$ is in $\mathbf{F}_p(\alpha)$ we have $\alpha^q = \alpha$. This means that $\alpha$ is a root of $X^q-X$, which in turn means that $X^q-X$ is a multiple of the minimal polynomial of $\alpha$. Since this minimal polynomial is $P$, this proves that $P$ is a factor of $X^q-X$, as desired.

We have shown that all the monic irreducible polynomials of degree $n$ are factors of $X^q-X$, and we now look for the others. Let $P \in \mathbf{F}_p[X]$ be an irreducible factor of $X^q-X$. Since $X^q-X$ has all its roots in $K$, $P$ must have all its roots in $K$ also. Let then $\alpha$ be a root of $P$ in $K$. We have $\mathbf{F}_p \subseteq \mathbf{F}_p(\alpha) \subseteq K$, with $[\mathbf{F}_p(\alpha):\mathbf{F}_p]$ being the degree of $P$. Since we also know that $[K:\mathbf{F}_p] = [K:\mathbf{F}_p(\alpha)]\times [\mathbf{F}_p(\alpha):\mathbf{F}_p] = n$, this shows that the degree of $P$ is a divisor of $n$.

Finally, let $d$ be a divisor of $n$ and $P$ be irreducible of degree $d$. Using again our previous result, we see that $P$ is a factor of $X^{p^d}-X$. Now, since $d$ divides $n$, we have that $p^d-1$ divides $p^n-1$, and so that $X^{p^d-1}-1$ divides $X^{p^n-1}-1$. Finally, $X^{p^d}-X$ divides $X^{p^n}-X$, and so since $P$ is a factor of $X^{p^d}-X$, it is a factor of $X^q-X$.

Algorithm

So, given a polynomial $A$, which we can assume squarefree, we know that the factors of $A$ of degree $d$ are precisely those who are also factors of $X^{p^d}-X$, but not of $X^{p^e}-X$ for any $e < d$. This leads to the following algorithm (again, discarding unit factors):

def ddfact(A):
    p = A.base_ring().characteristic()
    x = A.variables()[0]
    factors = []
    P = x^p
    d = 1
    while A.degree() > 0:
        T = gcd(P-x, A)
        if T != 1:
            factors.append((T, d))
        A = A//T
        d = d+1
        P = P^p
    return factors

which works as expected:

sage: f = (x+1)*(x+2)*(x^2+x+1)*(x^2+x+2)
sage: ddfact(f)
[(x^2 + 3*x + 2, 1), (x^4 + 2*x^3 + 4*x^2 + 3*x + 2, 2)]
sage: (x+1)*(x+2)
x^2 + 3*x + 2
sage: (x^2+x+1)*(x^2+x+2)
x^4 + 2*x^3 + 4*x^2 + 3*x + 2

Speeding it up

This works, but it is also very slow. For example on a polynomial of degree $20$:

sage: f = x^20 + 3*x^19 + 4*x^18 + 4*x^17 + x^16 + 3*x^15 + 2*x^14 + 2*x^13 + 3*x^12 + x^11 + 2*x^10 + 2*x^7 + 4*x^6 + 2*x^5 + 3*x^4 + 3*x^3 + x^2 + x + 2
sage: is_squarefree(f)
True
sage: t = cputime()
sage: ddfact(f)
[(x^2 + 2*x + 3, 2),
 (x^4 + 4*x^2 + 2, 4),
 (x^6 + 3*x^5 + 4*x^4 + 4*x^2 + x + 1, 6),
 (x^8 + 3*x^7 + 2*x^6 + x^5 + x^4 + 2*x^2 + x + 2, 8)]
sage: cputime(t)
1.932749000000058

it takes almost 2 seconds just for the distinct degree factorisation. For comparison, Sage’s factor() function on the same polynomial takes less than a tenth of a second for the complete factorisation:

sage: t = cputime()
sage: factor(f)
(x^2 + 2*x + 3) * (x^4 + 4*x^2 + 2) * (x^6 + 3*x^5 + 4*x^4 + 4*x^2 + x + 1) * (x^8 + 3*x^7 + 2*x^6 + x^5 + x^4 + 2*x^2 + x + 2)
sage: cputime(t)
0.02367099999992206

It is not difficult to identify the bottleneck: in order to compute the product of the irreducible factors of degree $d$, we take the GCD of $A$ and $X^{p^d}-X$. This means that we need to manipulate a polynomial ($X^{p^d}-X$) whose degree grows exponentially in $d$. In the previous example, the largest factor has degree $8$, which means $X^{p^d}-X$ will have degre $5^8 \approx 400,000$. This is clearly unacceptable.

It is also very simple to fix this problem. We are not really interested in $X^{p^d}-X$ itself, only in its common factors with $A$. It is easy to see that if $P$ is a common factor of $X^{p^d}-X$ and $A$, it will also be a factor of $X^{p^d}-X \mod A$ (the remainder of $X^{p^d}-X$ when divided by $A$), and conversely. In other words, we use the well-known result that $\gcd(A,B) = \gcd(A,B\mod A)$. the algorithm becomes

def ddfact(A):
    p = A.base_ring().characteristic()
    x = A.variables()[0]
    factors = []
    P = x^p
    d = 1
    while A.degree() > 0:
        T = gcd(P-x, A)
        if T != 1:
            factors.append((T, d))
        A = A//T
        d = d+1
        P = P^p % A
    return factors

which is much better:

sage: t = cputime()
sage: ddfact(f)
[(x^2 + 2*x + 3, 2),
 (x^4 + 4*x^2 + 2, 4),
 (x^6 + 3*x^5 + 4*x^4 + 4*x^2 + x + 1, 6),
 (x^8 + 3*x^7 + 2*x^6 + x^5 + x^4 + 2*x^2 + x + 2, 8)]
sage: cputime(t)
0.031130999999959386