(* File:  poliesModP4.m, Erhard Aichinger, 3.11.2006, 
   nur zur Benutzung fuer die Lehrveranstaltung MAT3
   der FH OOe, CMS. *)

Deg [f_, var_] := Length [ CoefficientList [ f, var ] ] - 1 ;

Lcf [f_, var_] := CoefficientList [f, var] [[Deg[f,var] + 1]];

PolynomialQuotientP [f_,g_,var_,p_] :=
Module[{lf, lg, f1, g1,q},
   f1 = PolynomialMod[f,p];
   g1 = PolynomialMod[g,p];
   If [Deg[g1, var] > Deg[f1, var],
         0,
      lf = Lcf [f1, var];
      lg = Lcf [g1, var];
      q = lf * PowerMod[lg, -1, p] * var^(Deg[f1,var]-Deg[g1,var]);
      (* compute return value *)
      q + PolynomialQuotientP [
             PolynomialMod [f1 - q * g1, p],
             g1,
             var,
             p]
   ]];
       

PolynomialRemainderP [f_,g_,var_,p_] :=
   PolynomialMod [f - g * PolynomialQuotientP [f,g,var,p], p]; 


GCDPoliesP [f_, g_, var_, p_] :=
Module [{matrix, linesofmatrix, q, gcdline},
   matrix = PolynomialMod [{{f,1,0},{g,0,1}},p];
   linesofmatrix = 2;
   While [Deg [ matrix [[linesofmatrix]] [[1]], var] >= 0,
      q = PolynomialQuotientP [matrix [[linesofmatrix - 1]] [[1]],
                               matrix [[linesofmatrix]] [[1]],
                               var, p];
      matrix = Append [matrix, PolynomialMod [Expand[Together[matrix [[linesofmatrix - 1]] -
                               q * matrix [[linesofmatrix]]]],p]];
      linesofmatrix++;
   ];
   gcdline = matrix [[linesofmatrix - 1]];
   gcdline = gcdline *  
             PowerMod [ 
               CoefficientList [gcdline [[1]], var] 
                 [[Deg [gcdline [[1]], var] + 1]], 
               -1,
               p
             ];
   i = 1;
   gcdline = PolynomialMod [gcdline, p];
   {matrix, gcdline}
]

   
PolynomialExtendedGCDP [f_, g_, var_, p_] :=
   Module [{lastline},
      lastline = GCDPoliesP[f,g,var,p][[2]];
      {lastline[[1]], {lastline[[2]], lastline[[3]]}}];


PolynomialInverseP [f_, g_, var_, p_] :=
(* Computes h such that f*h is congruent to 1 modulo g *)
  If [PolynomialExtendedGCDP[f,g,var,p] [[1]] =!= 1,
        "no inverse",
     PolynomialExtendedGCDP[f,g,var,p] [[2]] [[1]]];


PolyPowerModP [basis_, exp_, module_, var_, p_] :=
(* basis, exp, module are polynomials in Z_p.
   Computes basis^exp mod module by square and multiply method. *)
   If [exp == 0, 
         1,
      If [EvenQ [exp],
            basisSquare = PolynomialRemainderP [basis * basis, module, var, p];
            PolyPowerModP [basisSquare, (exp/2), module, var, p],
         PolynomialRemainderP [basis * PolyPowerModP [basis, exp-1, module,
                                                     var, p],
                               module, var , p]
      ]
  ]


IsIrreduciblePoly [f_, var_, p_] :=
Module [{n, i, f1,  canStillBeIrreducible, ggt, varToPToI},
   f1 = PolynomialMod [f, p];
   n = Deg [f1, var];
   canStillBeIrreducible = True;
   i = 1;
   While [canStillBeIrreducible && i <= (n/2),
      varToPToI = PolyPowerModP [var, p^i, f1, var, p];
      ggt = PolynomialExtendedGCDP [f1, varToPToI - var, var, p] [[1]];
      canStillBeIrreducible = (Deg [ggt, var] == 0);
      i = i+1;
   ];
   canStillBeIrreducible
]; 

RandomPoly [var_, maxdegree_, p_] := Module[{i},
                   Sum [Random[Integer, {0,p-1}] * var^i,
                       {i,0,maxdegree}]];

IsPrimitiveElement [poly_, module_, var_, p_] :=
Module[{primeFactors, n, canStillBePrimitive, j},
   n = p^Deg[module,var] - 1;
   primeFactors := Map [Function [x, x[[1]]], FactorInteger[n] ];
   canStillBePrimitive = True;
   j = 1;
   While [j <= Length [primeFactors] && canStillBePrimitive,
      canStillBePrimitive = 
         (PolyPowerModP [poly, n/primeFactors[[j]], module, var, p]
           =!= 1);
          
      j = j + 1;
   ];
   canStillBePrimitive
]


SeedRandom[Floor[AbsoluteTime[]]];




