(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
Great article! But, what about working in $F_q$, where q is a power of a prime number? Is the same theory valid? Thank you!
Yes, it’s essentially the same thing, it mostly just makes the implementation more messy.
Please correct me if I am wrong: if there are no common divisors up to the degree of the polynomial the while loop will continue incrementing d above the degree of the polynomial. Is it a correct behavior of the algorithm?
Or should it stop at d == A.degree() (with A the original polynomial, not the reduced one), otherwise it never exits (it may exit if $X^{p^d}−X = KA$ for some $d$, giving $A$ as GCD, but that’s not what we are looking for… amirite?)
Side note: every time T == 1 you can save time by skipping A = A//T (A should not change), right?
BTW, great article and even great series of articles!
PS: I am not sure about the syntax to use for the formulas, and there is no preview: crossing fingers…
“if there are no common divisors up to the degree of the polynomial”
That is not possible. Every factor of $A$ must have degree $d \le \deg(A)$, so it will be a common factor of $A$ and $X^{p^d}-X$.
Ok, sorry.
I am implementing the whole thing in c++ for polynomials of degree up to 512 in [latex]GF(2)[/latex], and the algorithm was failing only for polynomials of degree above 31…
It has been hard to find but the problem was in the computation of P = P^2 % A: I was using the same data size of the polynomial in input and the exponentiation was going in overflow. Solved by doubling the size of P (2 * A.size()).
Thank you.