@Carl Love 

Thank you for your answers and the time you are investing in my question. Thank you very much.

I have few other questions:

• Is there any way to work with CVS and Maple at the same time? Maple worksheets are internally coded in sth that looks like xml and Git doesn't conform to it.
• I'm mainly self-taught in Maple and I feel it's just HUGE. Do you know any way of getting used to Maple that isn't the help pages? I mean, the help pages are really precise and if you are looking for something concrete it's the way to go but I don't feel like reading all the help pages one by one is the most efficient way of learning this.

I think that's all. I won't have time to complete the program I was working in until next week but I had some time for testing things and I think it will work with mod.

@Carl Love 

Wow, I think I was too obsessed with modp1 and GF. I think I tried using mod but as I am relatively new to Maple, I might have introduced the parameters in wrong order or something and, as it failed (my fault indeed), I discarded the method. Now I have tried what you are proposing and read ?mod and I have a few questions.

Ceirtainly, what you exposed work and is far (really far) more simple than my solution. It looks exactly what I need. I'm gonna leave aside my previous work with wrappers and complex methods in favour of this one. But I have 2 problems:

• There are some functions that ceirtainly work with mod and aren't listed in ?mod. For example there exist Irreduc which says on ?Irreduc that is intended to work with mod but it's not listed in mod "allowed functions". How can I see all the options available for mod?
• My other question is that when I define two polynomials like the ones below, I supposed a*b mod p would also "reduce" the q's but instead it gives q^4. Shall I call Expand every time after an arithmetic calculation or am I doing something wrong?

restart; 
p := 2; 
n := 4;

irr := `mod`(Randprime(n, z), p);
                               4        
                       irr := z  + z + 1
alias(q = RootOf(irr));

a := q^3*x;
                                 3  
                           a := q  x
b := q*x;
                            b := q x
a*b mod p;
                              4  2
                             q  x 
Expand(%) mod p;
                              2    2
                           q x  + x 

• Finally, asking about your previous answer (in which we can continue the discussion from now on if you prefer in order that when someone arrives at this question in the future they see the related topics in the same answer) I don't quite understand what does this line do (see below). I get the text between quotes is the executing instruction but it's the first time I see the symbol $, don't know what 2: refers to and thought % was used to refer the previous result, never seen <%> before.

p||(1..2):= 'Randpoly(4,x,q) mod p' $ 2: <%>;

Thank you for your time

@Jasagredo 

It wasn't working correctly as when any x had a coefficient where alpha was in the denominator, the process would break.

This is the fix I worked out:

quot := proc(a::polynom, b::polynom)
			local aux, lista, pol, i, an, j, afq, bfq;
			afq := crear(a);
			bfq := crear(b);
			aux := quo(afq,bfq,x);
			lista := CoefficientList(aux, x);
			lista := map(r -> if not type(r,polynom) then numer(r)*g:-ConvertOut(g:-inverse(g:-ConvertIn(denom(r)))) else r end if, lista);
			lista := map(r -> if type(r, polynom) and not type(r, rational) then CoefficientList(r, alpha) else r end if, lista);
			pol := 0;
...
...

and the method goes on unchanged. The first map ensures that alpha will always be on the numerator of the coefficient by multiplying the numerator by the inverse (calculated inside the GF) of the denominator and that makes the trick, as shown in this example:

G16 := GF(2, 4, alpha^4+alpha+1);
                       G16 := 𝔽[16]
FX := fqpoly(G16);
                       FX := module() ... end module;
a := FX:-crear((alpha^3+1)*x^3);
                           /     3    \  3
                      a := \alpha  + 1/ x 
b := FX:-crear(alpha*x);
                          b := alpha x
FX:-crear(quo(a, b, x));
Error, invalid input: crear expects its 1st argument, e, to be of type polynom, but received x^2*(alpha^3+1)/alpha
quott := FX:-quot(a, b); remm := FX:-rem(a, b);
                         2 /     3        2    \
               quott := x  \alpha  + alpha  + 1/
                           remm := 0
a = FX:-crear(expand(b*quott+remm));
               /     3    \  3   /     3    \  3
               \alpha  + 1/ x  = \alpha  + 1/ x 

that works as expected. What do you think about this method?

@Carl Love 

I have tried adapting my code to modp2 representation in order to be able to do the division. It didn't work. I tried this:

p := 7; 
a := 3*x^6+2*x^2+x+5; 
b := 7*x^4+x^3+2*x+4;

tryquo := proc (a, b, p) 
    local afq, bfq, quott; 
    afq := modp2(ConvertIn(a, x, alpha), p); 
    bfq := modp2(ConvertIn(b, x, alpha), p); 
    quott := modp2(Pquo(afq, bfq, x, alpha), p); 
    return modp2(ConvertOut(quott, x, alpha), p) 
end proc;

tryquo(a,b,p);

and then it start throwing errors and goes crazy. It assigns 1mod7 to the x variable and does a lot of strange things. But I adopted your idea and think in the polynomials over Fq as polynomials in 2 variables. That way I worked out a method that breaks the polynom in a list of coefficients or coefficients list (in case the coeff is a polynomial in alpha or not) and then transform the result of quo.

It's easier to follow in the example:

with(PolynomialTools):
fqpoly := proc(g::symbol)
	option remember;
	local p;
	p := NumberTheory[PrimeFactors](g:-size())[1];
	module()
		export crear, suma, resta, mult, quot, rem;
		crear:= proc(e::polynom)
			local lista, r0, i;
			lista := CoefficientList(e, x);
			r0 := 0;
			for i to degree(e,x) + 1 do
				r0 := r0 + x^(i-1) * g:-ConvertOut(g:-ConvertIn(lista[i]))
			end do;
			return r0;
		end proc;
		suma := (a::polynom, b::polynom) -> crear(a+b);
		resta := (a::polynom, b::polynom) -> crear(a-b);
		mult := (a::polynom, b::polynom) -> crear(a*b);
		quot := proc(a::polynom, b::polynom)
			local aux, lista, pol, i, an, j, afq, bfq;
			afq := crear(a);
			bfq := crear(b);
			aux := quo(afq,bfq,x);
			lista := CoefficientList(aux, x);
			lista := map(r -> if type(r, polynom) and not type(r, rational) then CoefficientList(r, alpha) else r end if, lista);
			pol := 0;
			for i to numelems(lista) do
				an := 0;
				if type(lista[i],indexable) then
					for j to numelems(lista[i]) do
						if not type(lista[i][j],integer) then
							lista[i][j] := numer(lista[i][j])/denom(lista[i][j]) mod p;
						end if;
						an := an + lista[i][j]*alpha^(j-1);
					end do;
				else  
					if not type(lista[i],integer) then
						an:= numer(lista[i])/denom(lista[i]) mod p;
					else
						an := lista[i];
					end if;
				end if;
				pol := pol + an*(x^(i-1));
			end do;
			return crear(pol);
		end proc;
		rem := proc(a::polynom, b::polynom)
                        local quott, mu, res;
			quott := quot(a,b);
			mu := mult(b,quott);
			res := resta(a, mu);
			return res;
		end proc;
	end module;
end proc:

that is my re-definition of fqpoly. I shall warn that the name of the functions are in spanish and they mean: crear = create, suma = sum, resta = substract, mult = multiply, quot = quotient, rem = remainder.

What I do is obtain the quotient considering them as polynomials over the x variable and then I break the result in a coefficients list (for example alpha*1/2*x + 1 -> [1,[0,1/2]]) considering the first level as multiplied by x and the second level multiplied by alpha. then I run through all the coefficients and if they are fractional I do the division they imply but in modulus. I do that in the line reading:

lista[i][j] := numer(lista[i][j])/denom(lista[i][j]) mod p;

After that I build the polynom again and make the usual parse (with the crear method) of it as now its coefficients will be integer polynomials in alpha (and integer polynomials in alpha are elements of GF).

I know it's quite tricky but i think it works cause I could work out this two examples:

G7 := GF(7, 1, x+5);
                        G7 := 𝔽[7]
FX := fqpoly(G7);
                       FX := module() ...  end module;
a := FX:-crear(3*x^6+2*x^2+x+5);
                            6      2        
                    a := 3 x  + 2 x  + x + 5
b := FX:-crear(7*x^4+x^3+2*x+4);
                             3          
                       b := x  + 2 x + 4
quott := FX:-quot(a, b);
                                 3        
                     quott := 3 x  + x + 2
remm := FX:-rem(a, b);
                           remm := 4
a = FX:-crear(expand(b*quott+remm));
              6      2              6      2        
           3 x  + 2 x  + x + 5 = 3 x  + 2 x  + x + 5

No problem in there as the division is integer but this example shows a non-integer division:

a := FX:-crear(3*x);
                            a := 3 x
b := FX:-crear(2*x);
                            b := 2 x
FX:-crear(quo(a, b, x));
Error, (in ConvertIn) only integer polynomials in x can be converted
quott := FX:-quot(a, b);
                           quott := 5
remm := FX:-rem(a, b);
                           remm := 0
a = FX:-crear(expand(b*quott+remm));
                           3 x = 3 x

I think this can do the trick. What do you think?

@vv First of all, thanks for answering.

If I use your solution, I face the problem of the interpretation of the elements as (for example) for DUP elements and Gauss(Z) elements they both evaluate their type to function and so it's quite difficult to manage the overload of the proc (I'd have to use the Domain I'm working with inside the proc). I managed to get the euclides algorithm "working" with all the fields I intend to use without using the Domains package. Do you think this is a good way to proceed? (The interesting parts are the proc headings and the examples which are below EUCLIDES definition)

restart;
with(GaussInt, GIrem, GInormal):
EUCLIDES := overload([
	proc (a::integer, b::integer, p::integer) #Z, Zp
		local r0, r1, r2; 
		option overload;
		if nargs > 2 then 
			r0, r1 := abs(a) mod p, abs(b) mod p;
			while r1 <> 0 do 
				r2 := irem(r0, r1) mod p; 
				r0 := r1; 
				r1 := r2;
				end do; 
			return r0 mod p 
		end if;
		r0, r1 := abs(a), abs(b);
		while r1 <> 0 do 
			r2 := irem(r0, r1);
			r0 := r1;
			r1 := r2
			end do;
		return r0
		end proc, 
		
	proc (a::(complex(integer)), b::(complex(integer))) #GaussIntegers
		local r0, r1, r2; 
		option overload;
		
		r0, r1 := a, b; 
		while r1 <> 0 do 
			r2 := GIrem(r0, r1); 
			r0 := r1; 
			r1 := r2 
			end do;  
		return GInormal(r0)
	end proc,
			
	proc(a::polynom(rational), b::polynom(rational)) #Q[x], Fp[x]
		local r0, r1, r2, IDT;
		option overload;
		if nops(indets(a) union indets(b)) <> 1 then
			error("Polinomios en distintas variables.");
		end if;
		IDT := (indets(a) union indets(b))[1];
		r0, r1 := a, b;
		while degree(r1) > 0 do
			r2 := rem(r0, r1, IDT);
			r0:=r1;
			r1 := r2
			end do;
		return r0 / lcoeff(r0)
	end proc,

	proc(a::zppoly, b::zppoly, g::symbol) #Fp, Fq
		option overload;
		local r0, r1, r2;
		r0, r1 := a, b;
		while r1 <> g:-zero do
			r2 := modp1(Rem(r0,r1),NumberTheory[PrimeFactors](g:-size())[1]);
			r0 := r1;
			r1 := r2
		end do;
		return r0
	end proc,

	proc(a::polynom(algnum), b::polynom(algnum)) #Fq[x]
		option overload;
		local r0, r1, r2, IDT;
		if nops(indets(a) union indets(b)) <> 1 then
			error("Polinomios en distintas variables.");
		end if;
		IDT := (indets(a) union indets(b))[1];
		r0, r1 := a, b;
		while degree(r1) > 0 do
			r2 := rem(a,b,IDT);
			r0 := r1;
			r1 := r2
		end do;
		return r0
	end proc
	]):
	

#Ejemplo con Z
print(`Ejemplo con Z`);
a := 2;
b := 140;
mcd := EUCLIDES(a,b);

#Ejemplo con Zp
print(`Ejemplo con Zp`);
p := 3;
a := 8;
b := 140;
mcd := EUCLIDES(a,b,p);

#Ejemplo con GaussInt
print(`Ejemplo con GaussInt`);
p1 := Complex(2,-3):
p2 := Complex(5,2):
p3 := Complex(5,1):
a := p1*p2;
b := p1*p3;
mcd := EUCLIDES(a,b);

#Ejemplo con Q[x]
print(`Ejemplo con Q[x]`);
a := (x^2 + 2*x - 15)/7;
b := (x^2 + 13*x + 40)/5;
mcd := EUCLIDES(a,b);

#Euclides con Fp
print(`Ejemplo con Fp`);
p:=7:
G7:=GF(p,1,x+1);
a:=G7:-ConvertIn(6);
b := G7:-ConvertIn(3*x^6+2*x^2+x+5);
mcd:=EUCLIDES(a,b, G7);


#Euclides con Fp[x]
print(`Ejemplo con Fp[x]`);
a:= modp1(ConvertIn(6^2*x^2 + 2*6*x + 1,x),7);
b := modp1(ConvertIn(6+x,x),7);
mcd := EUCLIDES(a,b, G7);

#Euclides con Fq
print(`Ejemplo con Fq`);
p:=2: k:=4: minpol:=x^4+x+1:
G16 := GF(p,k,minpol);
a := G16:-ConvertIn(x^3+x^2);
b := G16:-`*`(a,a);
mcd := EUCLIDES(a, b, G16);

#Euclides con Fq[x]
print(`Ejemplo con Fq[x]`);
alias(alpha=RootOf(minpol)):
a := alpha^2*x^2+2*x*alpha+1;
b := alpha*x+1;
mcd := EUCLIDES(a,b);

If I use Domains I end up having to pass by parameter the Domain I'm working with and I don't really like it.

Thanks for your advice and please take a look at that code.