#
# GaloisField(p,k);
#
# Creates a Galois Field with p^k elements where the representation
# chosen is to use integers mod p if k = 1, otherwise to use Maples
# modp1 polynomial representation.  An irreducible polynomial is
# computed at random.
#
# Input/Output: This is done using integers, i.e. the elements
# of the field are encoded in [0,p^k-1]
#
# Author: MBM 1990
#

macro( `mod` = modp1 );
GaloisField := proc(p,k) local ext,F,env; option remember;

    if not type(p,integer) or not isprime(p) then
	ERROR(`1st argument must be a prime integer`) fi;

    if nargs = 1 then
	F := Zmod(p);
	F[DomainName] := GaloisField;
	RETURN( op(F) )
    fi;

    F := FiniteField(p,k);
    F[DomainName] := GaloisField;

    # Rep := Maple Modp1 over GF(p)
    addProperty( F, CannonicallyUniquelyRepresented );
    defOperation( Generator, [] &-> F, F );
    defOperation( Convert, [F] &-> Expression, F );
    defOperation( order, F &-> Integer, F );
    defOperation( isGenerator, F &-> Boolean, F );

    env := [modulus = p, degree = k, 'R' = F];

    F[Index] := subs(env, proc(n)
    	    [ConvertIn( convert(n,'base',modulus) ) mod modulus]
        end);
    F[Input] := eval(F[Index]);
    F[Lookup] := subs(env, proc(a) local x,t;
	    t := ConvertOut(a[1],x) mod modulus;
	    subs( x=modulus, taylor( t,x,infinity ) )
        end);
    F[Output] := eval(F[Lookup]);
    F[Random] := subs(env, proc() [Randpoly(degree-1) mod modulus] end);
    F[Convert] := subs(env, proc(a) local x;
	    ConvertOut(a[1],x) mod modulus
        end);

    ext := Randpoly(k) mod p;
    while not Irreduc(ext) mod p do ext := Randpoly(k) mod p od;
    ext := UNormal(ext) mod p;
    F[Extension] := [ext];

    F[`=`] := <evalb(x=y)>;
    F[0] := [Zero() mod p];
    F[1] := [One() mod p];
    F[`+`] := subs(env, proc() [Add(op(map(op,[args]))) mod modulus] end);
    F[`-`] := subs(env, proc(x,y)
	    if nargs = 1 then R[`-`](R[0],x)
	    else [Subtract(x[1],y[1]) mod modulus]
	    fi
        end);
    F[`*`] := subs(env, proc(x) local s,y;
	    if nargs = 0 then RETURN(R[0]) fi;
	    if type(x,integer) then s := modp(x,modulus) else s := x[1] fi;
	    for y in [args[2..nargs]] do
	        s := Rem(Multiply(s,y[1]),R[Extension][1]) mod modulus
	    od;
	    [s]
        end);
    F[Inv] := subs(env, proc(x) local s;
	    # if type(x,integer) then RETURN(R[Inv](R[Coerce](x))) fi;
	    if x = R[0] then ERROR(`division by zero`) fi;
	    if Gcdex(x[1],R[Extension][1],s) mod modulus = One() mod modulus
		then [s] else FAIL fi
        end);
    op(F)

end:

save `GF.m`;
quit
