(* File:  poliesModP.m *)

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]]];

