
(** FiniteFields.m  **)

(* Program by Ilan Vardi, December, 1989.

   Copyright Wolfram Research, Inc.
   P.O. Box 6059, Champaign, IL 61820, 217-398-0700   *)

(* Adjustment to Mathematica version 2 by
   Vedran Krcadinac, November, 1995.         *)

(*
   A finite field of p^n elements can be constructed as polynomials
   f[y] with coefficients (mod p), and the polynomials themselves
   are taken (mod Generator[y]), which is an irreducible polynomial
   of degree n (mod p). Multiplication and addition are just as usual
   but taken (mod p, Generator[y]).

   The only hard part is taking inverses (i.e., 1/f[y]), which ends up
   being equivalent to computing the extended GCD (mod p). In other words,
   given f[y], Generator[y], find a[y], b[y], such that

  a[y] f[y] + b[y] Generator[y] = GCD[Generator[y], f[y], Modulus-> p]

   If GCD = 1, then a[y] is the required inverse. This is the only
   nontrivial part. The rest of the code is just a routine set up of
   finite field rules.

   I set up a data structure FiniteField[f, y, {p, n}] which represents finite
   field elements. Once the rules of addition, multiplication,
   and division have been generalized to this data structure, other
   functions will operate on them like field elements. Thus
   PolynomialQuotient[f, g, x] will give the correct answer when
   applied to polynomials f, g, in x with coefficients FiniteField[f, y, {p, n}].

*)

(* This is the irreducible polynomial (mod p) that generates
   the finite field of p^n elements  *)

Generator[y_, p_, n_] := Generator[y, p, n] =
   First[Select[First /@
   FactorList[Cyclotomic[p^n-1,y], Modulus -> p],
   Exponent[#, y] == n&]]

(* An element in the finite field of p^n elements is a polynomial
   (mod Generator) and (mod p)                                       *)

FiniteField[f_, y_, {p_, n_}] :=
          FiniteField[PolynomialMod[PolynomialRemainder[PolynomialMod[f, p],
                   Generator[y, p, n], y], p],y, {p, n}] /;
     Exponent[f, y] >= n || Select[CoefficientList[f, y], # >= p&] != {}

(* The field laws for elements in a finite field. The only hard
   part is defining inverses, which is equivalent to finding the
   ExtendedGCD for polynomials (mod p)                           *)

FiniteField /: FiniteField[f_, y_, {p_, n_}] + FiniteField[g_, y_, {p_, n_}] :=
   FiniteField[PolynomialMod[f + g, p], y, {p, n}]

FiniteField /: FiniteField[f_, y_, {p_, n_}] + m_Integer :=
   FiniteField[PolynomialMod[f + m,p], y, {p, n}]

FiniteField /: FiniteField[f_, y_, {p_, n_}] * FiniteField[g_, y_, {p_, n_}] :=
   FiniteField[PolynomialMod[f * g, p] , y, {p, n}]

FiniteField /: FiniteField[f_, y_, {p_, n_}] * m_Integer :=
   FiniteField[PolynomialMod[f * m, p], y, {p, n}]

FiniteField /: FiniteField[0, y_, {p_, n_}] * x_ :=
   FiniteField[0, y, {p, n}]

(*    Taking 1/x, the only hard part    *)

FiniteField /: Power[FiniteField[f_, y_, {p_, n_}], -1] :=
   FiniteField[Last[FFExtendedGCD[Generator[y, p, n], f, y, p]], y, {p, n}]  /;
   PolynomialGCD[Generator[y, p, n], f, Modulus -> p] == 1

FiniteField /: FiniteField[f_, y_, {p_, n_}]^k_ :=
   FiniteField[PolynomialMod[f^k, p], y, {p, n}]    /;  k >= 0


(*   **********************************************************

     Computing ExtendedGCD (mod p). This is a little awkward
     since PolynomialRemainder[f, g, x], and PolynomialQuotient[f, g, x]
     don't have the option Modulus -> p, so if g is not monic
     you have to compute inverses (mod p) to make them monic.
     For example,

PolynomialRemainder[233 x^5 + 641 x^3 + 5 x + 17, 101 x^3 + 54 x^2 +67, x]

   =>  -(466108906/1030301) + (893999 x)/10201 - (549034337 x^2 )/1030301

    ************************************************************    *)


FFExtendedGCD[f_,g_, y_, p_] :=
    Reverse[FFExtendedGCD[g, f, y, p]] /;
            Exponent[PolynomialMod[f, p], y] < Exponent[PolynomialMod[g, p], y]

FFExtendedGCD[f_, g_, y_, p_Integer] :=
    {1, 0} /;  Exponent[PolynomialMod[g, p], y] == -Infinity

FFExtendedGCD[f_, g_, y_, p_] :=
     Block[{fp = PolynomialMod[f, p], gp = PolynomialMod[g, p],
            extf = {1, 0}, extg = {0, 1},
            q, r = y, leading},
       If[Exponent[gp, y] == 0, {0, PowerMod[gp, -1, p]},
           While[Exponent[r, y] > 0,
                 leading = Coefficient[gp, y, Exponent[gp, y]];
                  monic = PowerMod[leading, -1, p];
            q = PolynomialQuotient[PolynomialMod[monic fp, p], PolynomialMod[monic gp, p], y];
            r = leading *
                   PolynomialRemainder[PolynomialMod[monic fp, p], PolynomialMod[monic gp, p], y];
           {fp, gp, extf, extg}  = PolynomialMod[#, p] & /@
                                   {gp, r, extg, extf - q extg};
                 ];
             If[r == 0, Null, PowerMod[r, -1, p] extg]]]



(************************************************)
(*        Additions to FiniteFields.m           *)
(************************************************)

(* Vedran Krcadinac, March, 1996. *)

(* All elements of the field of order p^n *)

Poly[x_, a_] := a-1
Poly[x_, b__,a_] := a-1 + x Poly[x, b]

AllElements[x_,{p_,n_}] := AllElements[x,{p,n}] =
  Array[ FiniteField[Expand[Poly[x,##]],x,{p,n}] &, 
	 Table[p,{n}] ]  // Flatten

(* Define addition and multiplication in the field with
   p^n elements. For large fields this may take a while *)

Operations[{p_,n_}] :=
  Module[ {ff,i,j,x},
	  Clear[FFPlus];
	  Attributes[FFPlus]={Flat,Orderless,Listable};
	  Clear[FFTimes];
	  Attributes[FFTimes]={Listable};
	  FFOrder = p^n;
	  ff = AllElements[x,{p,n}];
	  For[i=1, i<=p^n, ++i,
	      For[j=1, j<=p^n, ++j,
		  FFPlus[i-1,j-1] = Position[ff,ff[[i]]+ff[[j]]] [[1,1]]-1;
		  FFTimes[i-1,j-1] = Position[ff,ff[[i]]*ff[[j]]] [[1,1]]-1
                 ]
             ]
        ]


(* Save the operations to a specified file *)

SaveFF[{p_,n_},name_] :=
 ( Operations[{p,n}]; Put[Definition[FFOrder,FFPlus,FFTimes],name] )


