Is there a fast way to factor polynomials with finite field $\operatorname{GF}(p^k)$ coefficients?

Question

For example the 33-term polynomial $x^{32} + c_{31} x^{31} + \ldots + c_1 x + c_0$, where the coefficients are 16-bit finite field numbers. For the cases I'm considering, the 33-term polynomial will have 16 quadratic factors of the form $x^2 + a x + b$, again where the coefficients are 16 bit finite field numbers. Currently I'm doing a somewhat optimized brute force over the nearly $2^{32}$ (4 billion) possible combinations of $a$ and $b$ ($a$ is $\geq 1$).

This can be used for finding a "compatible" (isomorphic) composite mapping of $\operatorname{GF}(2^{32})$ with primitive element $α(x) = x (\operatorname{hex} 2)$ to a quadratic of $\operatorname{GF}((2^{16})^2)$ with primitive element $β(x) = x (\operatorname{hex} 10000).$

Let $f(x)$ be the defining primitive polynomial for $\operatorname{GF}(2^{32}).$ If the 1 bit coefficients of $f(x)$ are treated as elements in $\operatorname{GF}(2^{16}),$ there will be 16 primitive quadratic factors of $f(x).$ Any of those 16 quadratic factors can be used for the mapping meeting the isomorphic requirements that $\operatorname{map}(a + b) = \operatorname{map}(a) + \operatorname{map}(b)$ and $\operatorname{map}(a b) = \operatorname{map}(a) \operatorname{map}(b).$

For $\operatorname{GF}(2^{32})$ to $\operatorname{GF}((2^{16})^2),$ even though a brute force search is nearly 4 billion cases, it can be done in a few minutes on a typical PC. I'm wondering if there is some faster way so that a case like mapping from $\operatorname{GF}(2^{64})$ to $\operatorname{GF}((2^{32})^2)$ could be searched for in a similar way.

With carryless multiply instructions like pclmulqdq on X86 or parallel table look up instructions like pshufb on X86, it may be a pointless mapping (primarily used to calculate $1/x$).

---

Example for $\operatorname{GF}(2^{16}),$ primitive polynomial $x^{16} + x^{12} + x^3 + x + 1.$ The coefficients of $f(x)$ are $0$ and $1,$ but are elements of $\operatorname{GF}(2^{16})$

f(x) = x^32 + x^22 + x^2 + x + 1

There are 16 quadratic factors, values shown in hex:

x^2 + 04c4 x + 118d
x^2 + 09ad x + 1cec
x^2 + 0e38 x + 1cb7
x^2 + 16b6 x + 1dbc
x^2 + 173b x + 0cf9
x^2 + 1c89 x + 1cf0
x^2 + 40ab x + 4be1
x^2 + 524f x + 0a76
x^2 + 5e62 x + 0716
x^2 + 5eec x + 0a37
x^2 + b67f x + 4188
x^2 + be6f x + fbf0
x^2 + e079 x + 17d4
x^2 + effe x + ed71
x^2 + f62b x + 07d5
x^2 + fd83 x + 17dd

---

An alternative for mapping from $\operatorname{GF}(2^{32})$ to $\operatorname{GF}((2^{16})^2)$ is to pre-select a primitive polynomial for $\operatorname{GF}((2^{16})^2)$, such as $x^2 + x + 2000_{16}$, then do a brute force search for any primitive element α(x) of $\operatorname{GF}(2^{32})$ based on polynomial $f(x) = x^{32} + x^{22} + x^2 + x + 1$, that will result in an isomorphic mapping. Turns out for that there are 32 instance of α(x) (out of 2^31 possible α(x)), but the search is much slower.

Answer 1

Once you have found one quadratic factor you can find the others (assuming they are all Galois conjugates) by applying the Frobenius Automorphism:

from sage.rings.finite_rings.hom_finite_field import FiniteFieldHomomorphism_generic

K.<T> = GF(2^16)

R.<x> = PolynomialRing(K)

p = x^32 + x^22 + x^2 + x + 1

# factor(p)

f = x^2 + (T^12 + T^10 + T^9 + T^6 + T^4 + 1)*x + T^15 + T^12 + T^10 + T^5 + T^4 + T + 1

Frob = K.frobenius_endomorphism(); Frob

Frobenius endomorphism T |--> T^2 on Finite Field in T of size 2^16

Frob(T^12 + T^10 + T^9 + T^6 + T^4 + 1)

T^13 + T^12 + T^11 + T^10 + T^9 + T^6 + T^5 + T^2 + 1

Frob(T^15 + T^12 + T^10 + T^5 + T^4 + T + 1)

T^14 + T^13 + T^11 + T^9 + T^8 + T^7 + T^6 + T^3 + T

---

$T^{12} + T^{10} + T^9 + T^6 + T^4 + 1$

Represents $1011001010001 = 0x1651$