First of all, thank you for your answers @vv, @Carl Love and @_Maxim_.

As a general answer I'll say that I tried reading the code of Domains and ¡Holy Jesus! I'm realtively new to Maple and that looked really messy. I did in fact read the implementation of GF (the normal GF, not the one in domains (shall I say the "standard"? is that correct?)) and I found out that GF is kind of a wrapper over modp1. I thought it had a more complex structure but it merely calls modp1 with the apropiated arguments.

I have struggled using the different approaches you advised me and I ended up using a "homemade" method I'll try to explain here.

First of all, we have to set the difference between elements of GF and polynomials over GF. An element in GF is an instance of type/zppoly wereas a polynom should be an instance of type/polynom(zppoly). This is not reachable as I tried the following for creating a simple polynom:

minpoly :=x^4+x+1;
G16 := GF(2, 4, minpoly);
a1 := G16:-ConvertIn(x^5+x^2);
poly := a1*x+2;

which throws Error, invalid terms in product: modp1(ConvertIn(x, x), 2).

Then I tried using the RootOf method. I defined alpha to be root of the minpoly and then I can multiply alpha by x but this gets a bit confusing as the alpha should be a GF element and then, when I try to do something like this:

alias(alpha = RootOf(x^4+x+1));
                                         alpha
alpha*x+2;
                                      alpha x + 2
G16:-ConvertIn(alpha);

which throws Error, (in ConvertIn) only integer polynomials in x can be converted. Also alpha*x+2 was of type/polynom(algnum). Hadn't seen before the type algnum.

That really confused me as the GF and the polynomial were using the same variable x. The inner variable and the one of the polynomials should ceirtainly not be the same as they refer to different things.

I realized that operations like adding, substracting, multiplying and powering with polynomials on integers can't arise with rational coefficients so I could "parse" the polynom in the beggining to make it's coefficients adjust to the field, then convert them out, make the operation and do the same conversion with the resultant polynom. I made a few changes to the definition of G16 and this is what I ended up with:

with(PolynomialTools):
fqpoly := proc(g::symbol)
	option remember;
	module()
		export ope;
		ope:= proc(e)
			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;
	end module;
end proc:

p:=2: k:=4: minpol:=alpha^4+alpha+1:
G16 := GF(p,k,minpol);
a := alpha^2*x^2+2*x*alpha+1;
b := alpha*x+1;

EUCLIDES := proc(a::polynom(rational), b::polynom(rational), g::symbol) #Q[x], Fp[x]
		local r0, r1, r2, lista, i, FX;
		option overload;
		if nargs = 3 then
			FX := fqpoly(g);
			r0 := FX:-ope(a);
			r1 := FX:-ope(b);
		elif nargs = 2 then
			r0, r1 := a, b;
		end if;
		while degree(r1) > 0 do
			r2 := rem(r0, r1, x);
			r0:=r1;
			r1 := r2
		end do;
		if nargs = 2 then
			return r0 / lcoeff(r0)
		end if;

		return r0
	end proc:

mygcd := EUCLIDES(a, b, G16);

fqpoly is a kind of wrapper over the in-and-out conversion with G16. Here I define two polynoms (a and b) but the coefficient of x^1 in the polynom a is not in G16 and has to be converted. This way, before executing the algorithm itself, I make a conversion and I end up with:

a := alpha^2*x^2+1

as it has to be (cause 2*alpha becomes 0 in G16). I thought about making fqpoly a user-defined type but didn't have time to do so. EUCLIDES works with this method (maybe it's just a matter of luck). I think this method is toilsome but effective. It's also interesting because I can use the defined operations for polynoms (such as rem in this example).

After this milestione, I'm trying to work out how to make the division of polynomials inside this method as it sometimes work and sometimes doesn't. This example demonstrates it:

G5 := GF(5,1);
Fx := fqpoly(G5);

a:=4*x^2;
a := Fx:-ope(a);
b:=2*x;
b:= Fx:-ope(b);

division := a/b;
                                       division := 2*x
division := Fx:-ope(division);
                                       division := 2*x

a:=3*x^2;
a := Fx:-ope(a);
b:=2*x;
b:= Fx:-ope(b);

division := a/b;
                                      division := 1.5*x
division := Fx:-ope(division);
Error, (in ConvertIn) only integer polynomials in x can be converted

I think i'll have to work out the division in fqpoly in order to respect the group modular behavior and then call this proc rather than the usual division of polynoms. This gets me to think that I'll also have to define other operations such as rem in order to preserve the integrity.