p := 3;
m := 2;
n := 3;

// When working with Z/p for p prime it is best to explicitly 
// construct it as a field, rather than using Integers(n)
//  Magma will then take greater advantage of the field properties.

// Both of these give finite fields.  You can also use the argument
// 25 instead of 5,2

G := GF(5, 2);

GG := FiniteField(5,2);


// Let's construct extensions of the field with p elements.
Fp := GF(p);

f := PrimitivePolynomial(Fp, m);

//  Now extend Fp using a root of f.

Fp_m<al> := ext<Fp | f>;

// Here is a function to get all the primitive polynomials;
// AllIrreduciblePolynomials() is a Magma function, but
// They don't seem to have AllPrimitivePolynomials
// So I wrote it.

AllPrimitivePolynomials := function(F, m)
   S := AllIrreduciblePolynomials(F, m);
   T := {@ @};
   for s in S do         
     if IsPrimitive(s) then
        Include(~T,s);        
     end if;
   end for;
return(T);
end function;

// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// Here we construct the field of p^m elements and then the field
// of p^mn.  In the line ff<x>, the <x> is not necessary.  
// It is just used to make ff a polynomial in x rather than 
// a less readable variable that Magma produces.

T := AllPrimitivePolynomials(Fp, m);
f := T[1];
Fp_m<al> := ext<Fp | f>;

TT := AllPrimitivePolynomials(Fp_m,n);
ff<x> := TT[1];
Fp_m_n<beta> := ext< Fp_m | ff>;

Fp_m; Fp_m_n;

// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//  We could also construct Fp_mn  directly.

S := AllPrimitivePolynomials(Fp,m*n);

// Let's  randomly select an element of S instead of using S[1].
h := Random(S);  
Fp_mn<gam> := ext< Fp | h>;
Fp_mn;


// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// Or we could construct Fp_mn by extending first to p^n then to the p^nm.


U := AllPrimitivePolynomials(Fp, n);
g := U[1];
Fp_n<alo> := ext<Fp | g>;

UU := AllPrimitivePolynomials(Fp_n,m);
gg<y> := UU[1];
Fp_n_m<beto> := ext< Fp_n | gg>;

Fp_n; Fp_n_m;

// Notice that f, g, h have $.1 as the variable, since I didn't specify
// the variable myself.
print f, ff, g, gg, h;

// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// You might play around with this function for the different fields.

DefiningPolynomial(Fp_n);

//%%%%%%%%%%%%%%%%%%%%%%%
// Here is a for loop to list all the primitive elements.

for a in Fp_m do
   if IsPrimitive(a) then
     print a;
   end if;
end for;

// Here is a function to list to list all of the roots of a 
// given polynomial f in a given finite field F.

FindRoots := function(F, f)
S :={};            // The set to contain the roots.
for a in Fp_m_n do 
   if Evaluate(h, a) eq 0 then
      Include(~S,a);
   end if;
end for;
return S;
end function;

//  Experiment with this function on the polynomials above, with the 
//  different fields constructed.

FindRoots(Fp_m_n, h);

//  These don't work because the fields are unrelated as 
//  far as Magma is concerned.

print Fp_mn ! beta;
print Fp_m_n ! gam;
      
//  We just tell Magma to relate the two fields.

Embed(Fp_mn, Fp_m_n);

print Fp_mn ! beta;
print Fp_m_n ! gam;

// Notice the result of the final one is of the roots of h
//  that we found using FindRoots.


//Experiment with these functions:
//Order
//FactoredOrder
//Primitive
//IsPrimitive