Polynomial factorization over finite fields with non-prime order

One can easily factor a polynomial over finite fields of prime order, using Factor command:

Factor[1 + x^2, Modulus -> 2]
===>  (1 + x)^2


Now, is it possible to do this over finite fields of non-prime order? As an example, $x^2 + x + 1$ is reducible over GF(4), and can be decomposed as $(x-a)(x-a^2)$, where $a\not\in \{0,1\}$ is a field member.

Problem is, members of non-prime order finite fields are not denoted by ordinary numbers; rather, they use either polynomial or vector notation.

-
See reference.wolfram.com/mathematica/FiniteFields/tutorial/… It says Factor will not work with the Finite Fields package. –  Michael E2 Jul 18 '13 at 23:54
@MichaelE2: So, Mathematica offers no way to factor over finite fields with non-prime order? –  Sadeq Dousti Jul 18 '13 at 23:56
No built-in functions that I know of. Others might. You could implement an algorithm or you could perhaps use MathLink to call PARI from within Mathematica: pari.math.u-bordeaux.fr/dochtml/mathlink.html See also this related question and the answer. –  Michael E2 Jul 19 '13 at 0:10
You may be able to call a GAP script. See GAP FAQ, Run, etc., this comment, the discussion of this question. –  Michael E2 Jul 19 '13 at 13:15
One way to go about this would be to implement the Cantor-Zassenhaus algorithm over finite extensions to prime fields. –  Daniel Lichtblau Jul 19 '13 at 22:40