performance tuning - How to do the polynomial stuff over finite fields extensions fast?

This question is raised from the problem of package FiniteFields being very slow (please, see the corresponding question): I have had an evidence that Mathematica takes the exponential time from count of multiplications/additions to compute, say, just the value of polynomial at specified point. Please, see the following example:

<< FiniteFields`;

Table[
  AbsoluteTiming@
    Sum[
      GF[2, 12][RandomInteger[{1, 2}, 12]] x^i,
      {i, 0, deg}
    ] //. {x -> GF[2, 12][{0, 1}]},
  {deg, 12, 22} 

]

For Intel(R) Core(TM) i5-2540M CPU @ 2.60 GHz I have received the following results of computation time (in seconds):

{0.0440, 0.0780, 0.1550, 0.3350, 0.7230, 1.5470, 3.1891, 6.7934, 14.3468, 30.1807, 64.2497}

I have to do common polynomial stuff over finite fields extensions. For instance, factorization (looking for roots) and multiplication. Does any fast solution exist for this?

Thank you!

Answer

I do not have functionality to cover anywhere near e.g. factorization (except over prime fields). I will provide code that I use for some of the simpler tasks though. I will assume that the reader (probably taking flight by the next line or so) is familiar with basics of representing elements in extrension fields, e.g. as polynomials of degree less than that of a given defining polynomial.

First we'll want a way to define an extension field. For that we need, at minimum, an irreducible of the correct degree.

randpoly[x_, p_, deg_] := 
 Table[Random[Integer, {0, p - 1}], {deg}].Table[
    x^j, {j, 0, deg - 1}] + x^deg

getIrreducible[x_, p_, deg_] := Catch[Module[{poly},
   While[True, poly = randpoly[x, p, deg];
    If[Length[FactorList[poly, Modulus -> p]] == 2, Throw[poly]];]]]

Example:

p = 293;

deg = 15;
SeedRandom[1111]
irred = getIrreducible[x, p, deg]

(* 38 + 117 x + 244 x^2 + 234 x^3 + 212 x^4 + 142 x^5 + 103 x^6 + 60 x^7 + 
 203 x^8 + 124 x^9 + 183 x^10 + 96 x^11 + 225 x^12 + 123 x^13 + 
 251 x^14 + x^15 *)

I use an explicit polynomial representation above but for practical purposes I'd work instead with the list of coefficients.

If you want to use the Zech logs representation then a primitive polynomial is required.

isPrimitive[lin_, poly_, p_, deg_, facs_] := 
 Catch[(For[j = 1, j <= Length[facs], j++, 
    If[PolynomialMod[lin^facs[[j]], {poly, p}] === 1, 
      Throw[False]];];
   True)]

getPrimitivePolynomial[x_, p_, deg_] := Catch[Module[
   {j, irred, facs},
   facs = (p^deg - 1)/Map[First, FactorInteger[p^deg - 1]];
   While[True,

    irred = getIrreducible[x, p, deg];
    If[isPrimitive[x, irred, p, deg, facs], 
     Throw[CoefficientList[irred, x]]]
    ];
   ]]

Example:

defpoly = getPrimitivePolynomial[x, p, deg]

(* {183, 163, 260, 257, 142, 15, 36, 65, 110, 80, 138, 143, 19, 188, 221,

  1} *)

Here I will show how to find a multiplicative inverse when we use defpoly to represent our extension field (I will use the standard representation, not Zech logs). We'll take a random element in the extension field.

r1 = randomListPoly[1, deg, p][[1]]

(* {142, 155, 24, 238, 167, 267, 184, 133, 264, 96, 171, 85, 280, 242, \
266} *)

Now I invert it.

r1inv = Algebra`PolynomialExtendedGCDModList[r1, defpoly, p][[2, 1]]


(* {34, 171, 230, 158, 48, 127, 202, 151, 212, 56, 257, 119, 168, 58, \
290} *)

Quick check using ordinary polynomial representation and PolynomialMod:

rpoly = r1.x^Range[0, Length[r1] - 1];
invpoly = r1inv.x^Range[0, Length[r1inv] - 1];
defpoly2 = defpoly.x^Range[0, Length[defpoly] - 1];
PolynomialMod[rpoly*invpoly, {defpoly2, p}]


(* 1 *)

Next we move on to univariate polynomials over this field GF(p^deg). Treat them as a coefficient list where each coeff is a "polynomial" representing an element in GF(p^deg). We use lists for these as well. The code below will multiply a pair of these polynomials modulo the defining equation (an irreducible) for the field and the characteristic.

polynomialMultiplyMod[l1_, l2_, irred_, p_] := 
 Module[{ndim1, ndim2, ntot, flatl1, flatl2, prod, res, 
   deg = Length[irred]},
  ndim1 = Dimensions[l1][[2]];
  ndim2 = Dimensions[l2][[2]];
  ntot = ndim1 + ndim2 - 1;
  flatl1 = Flatten[Map[PadRight[#, ntot] &, l1, {1}]];

  flatl2 = Flatten[Map[PadRight[#, ntot] &, l2, {1}]];
  prod = Algebra`PolynomialTimesModList[flatl1, flatl2, p];
  res = Partition[prod, ntot];
  Map[PadRight[
     Algebra`PolynomialRemainderModList[#, irred, p], deg - 1] &, 
   res]
  ]

polynomialAddMod[l1_, l2_, irred_, p_] := 
 Module[{len1 = Length[l1], len2 = Length[l2], minl, llong, 

   deg = Length[irred]},
  minl = Min[len1, len2];
  If[len1 == minl, llong = l2, llong = l1];
  Join[
   Map[PadRight[#, deg - 1] &, 
    Table[Algebra`PolynomialPlusModList[l1[[j]], l2[[j]], p], {j, 
      minl}]],
   Table[llong[[j]], {j, minl + 1, Max[len1, len2]}]]
  ]

Here is an example.

randomListPoly[m_, n_, p_] := RandomInteger[p - 1, {m, n}]
p = 2;
deg = 16;
r1 = randomListPoly[p^3, deg, p];
r2 = randomListPoly[p^3, deg, p];
prod = polynomialMultiplyMod[r1, r2, primpoly, p]

(* {{1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1}, {1, 0, 1, 1, 1, 0, 
  1, 0, 0, 1, 0, 1, 1, 0, 1, 1}, {1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 

  1, 1, 0, 1}, {1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1}, {0, 
  0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1}, {0, 0, 0, 0, 0, 1, 1, 
  0, 0, 1, 1, 0, 1, 0, 0, 1}, {1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 
  0, 0, 1}, {1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1}, {1, 0, 
  1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0}, {1, 1, 0, 0, 1, 0, 0, 1, 
  0, 0, 0, 1, 1, 0, 1, 1}, {1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 
  0, 0}, {1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0}, {0, 1, 0, 
  0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1}, {1, 0, 1, 0, 1, 1, 0, 1, 0, 
  0, 1, 0, 1, 0, 1, 0}} *)

I'll show this as a polynomial in x, with coefficients as polynomials in y.

Collect[
 FromDigits[Reverse[Expand[Map[FromDigits[Reverse[#], x] &, prod]]], 
  y], x]

(* 1 + y + y^2 + y^3 + y^6 + y^7 + y^8 + y^9 + y^10 + y^11 + \
y^13 + x^9 (y + y^2 + y^4 + y^5 + y^6 + y^7 + y^8 + y^10) + 
 x^8 (y^2 + y^3 + y^6 + y^7 + y^8 + y^11) + 
 x^3 (y + y^2 + y^3 + y^8 + y^10 + y^11) + 
 x^15 (1 + y + y^2 + y^3 + y^4 + y^5 + y^6 + y^7 + y^9 + y^12) + 

 x (y^2 + y^6 + y^7 + y^9 + y^10 + y^12) + 
 x^6 (y + y^2 + y^4 + y^5 + y^10 + y^11 + y^12) + 
 x^13 (1 + y^2 + y^3 + y^7 + y^10 + y^11 + y^12) + 
 x^11 (y + y^3 + y^7 + y^8 + y^9 + y^10 + y^11 + y^12) + 
 x^10 (y^2 + y^3 + y^4 + y^5 + y^6 + y^7 + y^8 + y^11 + y^13) + 
 x^14 (y + y^4 + y^9 + y^11 + y^13) + 
 x^2 (y + y^2 + y^8 + y^10 + y^11 + y^13) + 
 x^12 (y + y^2 + y^3 + y^4 + y^5 + y^6 + y^7 + y^8 + y^9 + y^10 + 
    y^11 + y^13) + 
 x^5 (1 + y^4 + y^5 + y^7 + y^8 + y^10 + y^12 + y^13) + 

 x^7 (1 + y^2 + y^3 + y^4 + y^6 + y^9 + y^11 + y^12 + y^13) + 
 x^4 (y + y^2 + y^4 + y^8 + y^9 + y^10 + y^11 + y^12 + y^13) *)

So this is a start. I have not written code to do polynomial division. A good way to go about this can be found in

J. von zur Gathen and J. Gerhard. Modern Computer Algebra. Ch. 9.

An alternative is to convert to explicit polynomials and use a certain Groebner basis tactic to invert the divisor modulo a power of the polynomial variable, then multiply that result by the dividend. This is, suffice it to say, much slower. It is useful if you will be working many times with the same divisor.

A third tactic is to do a pedestrian polynomial one-term-at-a-time iteration. Somewhat slower, but basically the same idea, is to convert to explicit polynomials and use PolynomialReduce, giving the extension defining polynomial as a second divisor and Modulus->2.

Likewise one may want a gcd or extended gcd. This can most readily be coded by iterated quotient/remainder sequence.

I will finally mention that full blown factorization is not needed if you only seek roots in the base field (that is, your GF(p^n) rather than a proper extension thereof). In that case one only needs a first step of distinct-degree factorization, to gather all roots of degree 1 over the base field.