Attach("RelAlgKTheory.m");
Attach("ao.m");
ZX<x> := PolynomialRing(Integers());

load "A4Mipos.m";
load "D5Mipos.m";
load "D6Mipos.m";
load "D7Mipos.m";

/* Now you have the following list of polynomials:
   A4_pols, D5_pols, etc
*/

/* Compute the associated order and a generating element for the j-th polynomial
   in A4_pols
*/
j := 111;
X := A4_nb(j, A4_pols);
L := X[1]; theta := X[2]; Atheta := X[3]; AssOrd := X[4]; nb := X[5]; h := X[6]; QG := X[7];

/* Compute the associated order and a generating element for the j-th polynomial
   in D5_pols
*/
j := 11;
X := Dn_nb(j, D5_pols);
L := X[1]; theta := X[2]; Atheta := X[3]; AssOrd := X[4]; nb := X[5]; h := X[6]; QG := X[7];

/* Compute the associated order and a generating element for an abelian extension.
*/
f  := CyclotomicPolynomial(25);
X := Abelian_nb(f);



f := x^4 + 3*x^2 - 2*x + 6;   
L := SplittingField(f);

/* The following sequence of statements checks whether OL is locally free. It suffices to
   check local freeness for primes p dividing the group order.
*/
G, Aut, h := AutomorphismGroup(L);
h := map<Domain(h)->Codomain(h) | g:->h(g^-1)>;
OL := MaximalOrder(L);
theta := NormalBasisElement(OL, h);
Ath := ComputeAtheta(OL, h, theta);
QG := GroupAlgebra(Rationals(), G);
AssOrd := ModuleConductor(QG, Ath, Ath);
rho := RegularRep(QG);
M := ZGModuleInit(Ath`hnf, rho);
for p in PrimeDivisors(Order(G)) do
    isfree, w := IsLocallyFree(QG, AssOrd, M, p);
    if isfree then
        print "OL is free over its associated order at p = ", p;
    else
        print "OL is NOT free over its associated order at p = ", p;
    end if;
end for;

/* The following poynomials have Galois group S_4. Note that in Example 7 and 8
   OL is not free over the associated order.
*/

/* unramified */
f := x^4 +x +1;     /* Example 1 */
f := x^4 +3*x +1;   /* Example 2 */

/* 2 wildly ramified */
f := x^4 + 2*x + 2; /* Example 3 */
f := x^4 + 6*x + 6; /* Example 4 */

/* 3 wildly ramified */
f := x^4 - x^3 -6*x^2 + 3; /* Example 5 */

/* 2 and 3 wildly ramified */
f := x^4 - 2*x^3 + 3*x^2 +2*x - 1;  /* Example 7 */
f := x^4 + 3*x^2 - 2*x + 6;         /* Example 8 */


/* Compute a normal basis element in the galois closure of Q[x]/(f(x)) */
f := x^4 + x + 1;
X := S4_nb(f);
L := X[1]; theta := X[2]; Atheta := X[3]; AssOrd := X[4]; nb := X[5]; h := X[6]; QG := X[7];

/* The following sequence of statements checks whether the computed results are correct.
   It can be used with any of the above examples.
*/
lambda := QG ! ElementToSequence(nb);
delta := QGAction(lambda, theta, h);    /* the integral normal basis element as an element of L */
IsIntegral(delta);
OL := MaximalOrder(L);
delta := OL ! delta;
conj := [ElementToSequence( OL ! QGAction(QG!e, OL ! delta, h)) : e in RowSequence(AssOrd`hnf) ];
Z := Matrix(conj);
Determinant(Z);           /* this determinant must be +-1 */