# FiniteFields package is very slow. Any fast substitute for Mathematica?

I want to compute the inverse of matrix, say with dimensions $100 \times 100$, defined over a large finite field extension such as $GF(2^{120})$. I am using the package FiniteFields, but Mathematica's computation time is exponential with respect to matrix dimensions. The following code illustrates the problem:

<< FiniteFields;
Table[
With[{ext = 12},
First@AbsoluteTiming@
Inverse[
Array[GF[2, ext][RandomInteger[{0, 1}, ext]] &, {n, n}]
]
],
{n, 1, 11}
]


I am using an Intel Xeon X5680 @ 3.33GHz (64-bit OS) and Mathematica v8.0.4.0. I have received the following timing results:

{0.0030, 0.0080, 0.0210, 0.0630, 0.1860, 0.5110, 1.3350, 3.3840, 8.9340, 23.0090, 57.4660}


I believe the source of the problem is that the FiniteFields package defines many UpValues, DownValues and SubValues of Times and Plus for head GF and, consequently, the pattern matching of arguments is degraded.

Does anyone know if a patch for the FiniteFields package or a faster substitute providing a similar interface?

Many thanks!

-

--- edit ---

More recent versions of the work containing the code mentioned below are found at these links: SymbolicFAQ and PCwGB

--- end edit ---

This will not scale to dimension 100 but will be an improvement on what you now have. It is cribbed from the section "Linear Algebra over Galois Fields here as well as the section "Groebner bases over modules and related computations" in this notebook.

deg = 12;
flen = 3;
j = 0;
While[flen > 2 && j++ < 100,
defpoly =
x^deg + 1 + RandomInteger[{0, 1}, deg - 1].x^Range[1, deg - 1];
flen = Length[FactorList[defpoly, Modulus -> 2]]
]

dim = 20;
mat = Table[
RandomInteger[{0, 1}, deg].x^Range[0, deg - 1], {dim}, {dim}];

newvars = Array[z, 2*dim];
augmat = Transpose[
ArrayFlatten[{mat, IdentityMatrix[Length[mat]]}, 1]];
auxpolys = Union[Flatten[Outer[Times, newvars, newvars]]];
linpolys = Join[augmat.newvars, {defpoly}, auxpolys];
allvars = Append[newvars, x];


Here is the bulk of the computation.

Timing[gb = GroebnerBasis[linpolys, allvars, Modulus -> 2];]

(* {168.240000, Null} *)


Now extract the inverse matrix.

modulegb = Complement[gb, Join[auxpolys, {defpoly}]];
redmat = Reverse[Sort[Outer[D, modulegb, newvars]]];
invmat = Transpose[redmat[[All, dim + 1 ;; 2*dim]]];


We'll check the result.

PolynomialMod[invmat.mat, {defpoly, 2}] === IdentityMatrix[dim]

(* True *)


A serious bottleneck is the need, using this approach, to have $O(\tt{dim}^2)$ auxiliary polynomials.

One can improve on this by writing a direct row reduction using a list representation form the field elements. This would require working with undocumented functionality in the Algebra context, such as AlgebraPolynomialTimesModList. This TMJ article may give slightly more information should you opt to go that route.

--- edit ---

Here is an approach that will sometimes work. Treat the elements as algebraics over Q instead of GF[2]. Now you can use some built in functionality. If you are lucky you get substantially the same row reduction and can recover the mod-2 result. I'll show what I mean with the example above.

We first set up a matrix of AlgebraicNumber opjects.

alg = algnum[root[defpoly], CoefficientList[#1, x]] &;
mat2 = Map[alg, mat, {2}] /. x -> # /. root[a_] :> Root[a &, 1] /.
algnum -> AlgebraicNumber;
augmat = Transpose[
ArrayFlatten[{mat2, IdentityMatrix[Length[mat2]]}, 1]];


Now row reduce it and hope we do not get any denominators divisible by 2.

Timing[
rred = RowReduce[augmat, Method -> "OneStepRowReduction"];]

(* {36.53499999999994, Null} *)


Extract the inverse part.

invmat = Transpose[
rred[[All, dim + 1 ;; 2*dim]] /.
AlgebraicNumber[aa_, bb_] :>
PolynomialMod[bb, 2].x^Range[0, deg - 1]];


Check:

PolynomialMod[invmat.mat, {defpoly, 2}] === IdentityMatrix[dim]

(* True *)


--- end edit ---

-
Thank you for great references! It seems that I has to reimplement finite fields stuff in Mathematica to improve performance. This is bad news for me :( – Piotr Semenov Nov 25 '12 at 11:28
I think the key is AlgebraPolynomialExtendedGCDModList because it allows you to find field inverses. You invoke it using polynomials represented as their coefficient lists (constant term on the left and ascending from there). If defpoly is the list for your field defining polynomial, then to invert the field element r1 you would do r1inv=AlgebraPolynomialExtendedGCDModList[r1,defpoly,p][[2,1]]. In this notation p is the prime charatceristic, so it would be 2 for the cases of apparent interest. – Daniel Lichtblau Nov 26 '12 at 19:04
@Daniel Lichtblau Thank you for the point! I see, the easiest ways are to implement finite field extensions myself or do the patch for FiniteFields. If I get some progress, I will share it in the community. – Piotr Semenov Nov 27 '12 at 10:11
Yes, please do share with community. Here and/or library.wolfram.com would be good places. Is your work only with linear algebra on field elements, or does it involve e.g. polynomials over extension fields? if the latter, there may be some existing functionality of interest I could locate. – Daniel Lichtblau Nov 27 '12 at 15:35
@DanielLichtblau My work includes both, the linear algebra and polynomials over extension fields. For example, I have to compute the determinant of big matrices over field extensions and to search roots of polynomials by their factorization. Sure, it will be great! Thanks. – Piotr Semenov Dec 3 '12 at 12:32

The textbook algorithm, that is, row reducing {A, I} to {I, A^(-1)} works.

n = 1000;
m = Array[RandomInteger[{0, 1}]&, {n, n}];
m // MatrixPlot

MatrixRank[m, Modulus -> 2] (* see if m is full rank *)

inverse = With[{n = Length@#},
RowReduce[MapThread[Join, {#, IdentityMatrix[n]}], Modulus -> 2][[;;,n+1;;]]]&;
inverseM = inverse[m]; // AbsoluteTiming

Mod[inverseM. m, 2] // MatrixPlot


On my laptop, computing inverse takes 15.262331 for a full rank matrix.

Often, you use the inverse matrix to solve A x == b, and Mathematica has fast LinearSolve in the FiniteFields package.

n = 1000;
LinearSolve[Array[RandomInteger[{0, 1}]&, {n, n}],
Array[RandomInteger[{0, 1}]&,n], Modulus -> 2] ;// AbsoluteTiming


On my laptop this last gives {5.110337, Null}`

-
Works well over prime fields but not extensions thereof. – Daniel Lichtblau May 27 '13 at 19:46