#
# AlgebraicExtension(D,m)
#
# Creates a simple algebraic extension defined by m in D, the minimal
# polynomial, an irreducible univariate polynomial over a field F
#
# Author: MBM 1989
#
AlgebraicExtension := proc() local d,k,A,F,Rep,e,env;

	if nargs <> 2 then ERROR(`wrong number of arguments`) fi;
	Rep := args[1];
	if not hasCategory(Rep,UnivariatePolynomial) then 
		ERROR(`1st argument must be a univariate polynomial`) fi;

	F := Rep[CoefficientRing];
	if not hasCategory(F,Field) then
	ERROR(`1st argument must be a univariate polynomial over a field`) fi;

	e := args[1][Normal](args[2]);
	d := Rep[Degree](e);
	if d = 0 then ERROR(`extension cannot be a constant`) fi;
	env := ['P' = Rep,'D' = 'A'];

	k := d*F[AbsoluteDegree];
	if hasCategory(F,Finite) then
		A := FiniteField(F[Characteristic],k);
	else	A := Field();
	fi;
	A[AbsoluteDegree] := k;
	A[DomainName] := AlgebraicExtension;
	defOperation( BaseField, [] &-> Field, A );
	defOperation( Extension, Rep, A );
	defOperation( Degree, Integer, A );
	defOperation( Representation, [] &-> UnivariatePolynomial, A );

	A[Extension] := e; # Should check that e is irreducible over F
	A[BaseField] := Rep[CoefficientRing];
	A[Representation] := Rep;
	A[Characteristic] := Rep[Characteristic];
	A[Degree] := d;
	A[Random] := subs(env, proc() P[Random](D[Degree]-1) end);
	A[Type] := eval(Rep[Type]);
	A[Output] := eval(Rep[Output]);
	A[Input] := subs(env, proc(x) local i;
		i := P[Input](x);
		if i = FAIL then i else P[Rem](i,D[Extension]) fi
	    end);

	A[0] := Rep[0];
	A[1] := Rep[1];
	A[`=`] := eval(Rep[`=`]);
	A[`+`] := eval(Rep[`+`]);
	A[`-`] := eval(Rep[`-`]);
	A[`*`] := subs(env, proc(x) local r,y;
		if nargs = 0 then RETURN(P[1]) fi;
		if type(x,integer) then r := P[Coerce](x) else r := x fi;
		for y in [args[2..nargs]] do
		    r := P[Rem](P[`*`](y,r),D[Extension])
		od;
		r
	    end);
	A[Inv] := subs(env, proc(x) local s;
		if x = P[1] then RETURN( x ) fi;
		if P[Gcdex](x,D[Extension],s) = P[1] then s else FAIL fi
	    end);

#A[ModularHomomorphism] := SAEModularMapping1(Rep,A);
#A[ModularMapping] := SAEModularMapping2(Rep,A);

	op(A)
end:

SAEModularMapping1 := proc(P,F) local e,x,i,m,b,g;
    if not isDomain(F[BaseField],Rational) then RETURN(proc() FAIL end) fi;
    e := F[Extension];
    e := convert( [ seq(e[i] * x^(i-1), i=1..nops(e)) ], `+` );
    if not type(e,polynom,integer) then RETURN(proc() FAIL end) fi;
    b := 10001;
    do
        m := abs( subs(x=b,e) );
#	Product of all primes in the range [0,100]
	g := igcd(2305567963945518424753102147331756070,m);
	while m > 1 and g > 1 do
		m := iquo(m,g);
		g := igcd(2305567963945518424753102147331756070,m);
	od;
#	Product of all primes in the range [100,1000]
g := igcd( 8496969489233418110532339909187349965926062586648932736611545426342\
203893270769390909069477309509137509786917118668028861499333825097682386722983\
737962963066757674131126736578936440788157186969893730633113066478620448624949\
257324022627395437363639038752608166758661255956834630697220447512298848222228\
550062683786342519960225996301315945644470064720696621750477244528915927867113,
m);
	while m > 1 and g > 1 do
		m := iquo(m,g);
g := igcd( 8496969489233418110532339909187349965926062586648932736611545426342\
203893270769390909069477309509137509786917118668028861499333825097682386722983\
737962963066757674131126736578936440788157186969893730633113066478620448624949\
257324022627395437363639038752608166758661255956834630697220447512298848222228\
550062683786342519960225996301315945644470064720696621750477244528915927867113,
m);
	od;
	if m < b then b := 2*b+1; next fi;

lprint(`modulus chosen`,b);
RETURN( subs( _P = P, _M = m, _B = b, 
   proc()
	proc(a) local d,r,k;
		d := _P[Degree](a);
		r := _P[Coeff](a,d);
		for k from d-1 by -1 to 0 do
		    r := r*_B+_P[Coeff](a,k) mod _M
		od;
		r
	    end,
	_M
   end ) )

    od;
end:

SAEModularMapping2 := proc(P,F) local e,f,G,x,i,m;
    if not isDomain(F[BaseField],Rational) then RETURN(proc() FAIL end) fi;
    e := F[Extension];
    e := convert( [ seq(e[i] * x^(i-1), i=1..nops(e)) ], `+` );
    m := 1;
    do
        m := modp1( Prime(m) );
        while lcoeff(e) mod m = 0 do m := modp1( Prime(m) ) od;
	f := Factors(e) mod m;
	# Choose an interesting factor for now
	f := f[2];
	f := f[1][1];
lprint(`factor chosen`,f);
	G := GaloisField(m,degree(f),f);
RETURN( subs( _M=m, _P=P, _f=modp1(ConvertIn(f,x),m), _G=G, proc()
	proc(a) local r;
		r := _P[ListCoeffs](a) mod _M;
		r := modp1( ConvertIn(r), _M );
		[modp1( Rem(r,_f), _M  )]
	    end, _G end ) )
    od;
end:

save `AE.m`:
quit
