Singular

http://www.singular.uni-kl.de/Manual/latest/sing_714.htm#SEC766

Edit: see the [b] Note[/b]

http://www.singular.uni-kl.de/Manual/latest/sing_216.htm#SEC256

Singular can (at the moment) not factorize over GF(p^n) for baserings
defined as [i]ring r =(p^n,a),.... [/i]

So you have to choose the fir st form by setting the minpoly yourself.
[code]
> ring ra5 = (25,a),(x,y,z),dp;
> poly f = (x+y)*(z-x);
> factorize(f);
? not implemented
? error occurred in or before STDIN line 9: `factorize(f);`
> basering;
// # ground field : 25
// primitive element : a
// minpoly : 1*a^2+4*a^1+2*a^0
// number of vars : 3
// block 1 : ordering dp
// : names x y z
// block 2 : ordering C
> ring ra25 = (5,a),(x,y,z),lp;
> minpoly = 1*a^2+4*a^1+2*a^0;
> poly f = (x+y)*(z-x);
> factorize(f);
[1]:
_[1]=1
_[2]=x+y
_[3]=-x+z
[2]:
1,1,1
[/code]

[code]
> LIB "general.lib";
> cyclic(3);
_[1]=x+y+z
_[2]=xy+xz+yz
_[3]=xyz-1
> ideal I = _;
> option(redSB);
> facstd(I);
[1]:
_[1]=z+(-2a-1)
_[2]=y+(2a+2)
_[3]=x-1
[2]:
_[1]=z-1
_[2]=y+(-2a-1)
_[3]=x+(2a+2)
[3]:
_[1]=z+(2a+2)
_[2]=y-1
_[3]=x+(-2a-1)
[4]:
_[1]=z-1
_[2]=y+(2a+2)
_[3]=x+(-2a-1)
[5]:
_[1]=z+(2a+2)
_[2]=y+(-2a-1)
_[3]=x-1
[6]:
_[1]=z+(-2a-1)
_[2]=y-1
_[3]=x+(2a+2)
[/code]
Then check also
[i]triang.lib[/i]
http://www.singular.uni-kl.de/Manual/latest/sing_1480.htm#SEC1556

[code]
> triangLfak(std(I));
[1]:
_[1]=z-1
_[2]=y+(-2a-1)
_[3]=x+(2a+2)
[2]:
_[1]=z-1
_[2]=y+(2a+2)
_[3]=x+(-2a-1)
[3]:
_[1]=z+(-2a-1)
_[2]=y-1
_[3]=x+(2a+2)
[4]:
_[1]=z+(-2a-1)
_[2]=y+(2a+2)
_[3]=x-1
[5]:
_[1]=z+(2a+2)
_[2]=y+(-2a-1)
_[3]=x-1
[6]:
_[1]=z+(2a+2)
_[2]=y-1
_[3]=x+(-2a-1)
[/code]

See the discussion in the forum:

[i]Solving systems of polynomials over GF(2)[/i]
http://www.singular.uni-kl.de/forum/viewtopic.php?f=10&t=1908

combined with

[i]...declaration of ring or galois field[/i]
http://www.singular.uni-kl.de/forum/viewtopic.php?f=10&t=1823

As part of my master's thesis I started to write a solver for systems of multivariate polynomial equations over finite fields.

I am curious if there is any work in progress in this field or in finding roots of univariate polynomials with coefficients in GF(p^n)?