AmodpModule := recformat<M : ModTupFld,
                         phi : Map,
                         lift : Map
                    >;

/******************************** NormalBasisElement ***************************************
 *
 * OL is the ring of algebraic integers in the absolute number field L, h is a map as
 * provided by  
 * 		G, Aut, h := AutomorphismGroup(L);
 *              h := map<Domain(h)->Codomain(h) | g:->h(g^-1)>;
 * 
 * The function "computes" a normal basis element by randomly testing some integral elements.
 * This is very dirty, but works well.
 *
 *******************************************************************************************/

intrinsic NormalBasisElement(OL :: RngOrd, h :: Map) -> RngOrdElt
{}
    local G, b, D, found;

    G := Domain(h);

    for b in Basis(OL) do
        D := Matrix([ ElementToSequence( h(g)(b) ) : g in G ]);
	if Determinant(D) ne 0 then
	    return OL!b;
	end if;
    end for;

    found := false;
    while  not found  do
        b := OL ! [ Random(3) : i in [1..#G] ];
	D := Matrix([ ElementToSequence( h(g)(b) ) : g in G ]);
        if Determinant(D) ne 0 then
	   return OL!b;
	end if;
    end while;

    if not found then
        print "ERROR: No normal basis element found!!!";
    end if;
end intrinsic;


/******************************** ComputeAtheta ***************************************
 *
 * OL is the ring of algebraic integers in the absolute number field L, h is a map as
 * provided by
 *              G, Aut, h := AutomorphismGroup(L);
 *              h := map<Domain(h)->Codomain(h) | g:->h(g^-1)>;
 * theta is a normal basis element.
 *
 * Computes the Z-module Atheta := { lambda in QG | lambda(theta) in OL}.
 *
 *******************************************************************************************/

intrinsic ComputeAtheta(OL :: RngOrd, h :: Map, theta :: RngOrdElt) -> Rec
{}
    local n, A, b, i, x, Atheta;;

    G := Domain(h);
    n := #G;
    A := Matrix( [ElementToSequence( OL ! h(g)(theta)  ) : g in G] );
    b := ZeroMatrix(Rationals(), 1, n);
    Atheta := [];
    for i:=1 to n do
        b[1,i] := Rationals() ! 1;
        x := Solution(A, b);
	Append(~Atheta, x);
	b[1,i] := Rationals() ! 0;
    end for;
    Atheta := ModuleInit(Matrix(Atheta));
    return Atheta;
end intrinsic;



/******************************** QGAction ***************************************
 *
 * lambda is an element in QG, a an element in OL, and h a map 
 * provided by
 *              G, Aut, h := AutomorphismGroup(L);
 *              h := map<Domain(h)->Codomain(h) | g:->h(g^-1)>;
 *
 * Computes lambda(a), which is an element of L.
 *
 *******************************************************************************************/

intrinsic QGAction(lambda :: AlgGrpElt, a :: FldOrdElt, h :: Map) -> FldOrdElt
{}
    
    local g;

    return &+[ Coefficient(lambda, g) * h(g)(a) : g in Domain(h) ];

end intrinsic;
    

intrinsic QGAction(lambda :: AlgGrpElt, a :: RngOrdElt, h :: Map) -> RngOrdElt
{}

    local g;

    return &+[ Coefficient(lambda, g) * h(g)(a) : g in Domain(h) ];

end intrinsic;


/******************************** TestAtheta ***************************************
 *
 * Checks if the result of ComputeAtheta is correct. Only for testing.
 *
 *******************************************************************************************/
    
intrinsic TestAtheta(QG :: AlgGrp, Atheta :: Rec, OL :: RngOrd, theta :: RngOrdElt, h ::Map) -> FldOrdElt
{}
    local w, b;

    w := [QGAction(QG ! b, theta, h) : b in RowSequence(Atheta`hnf)];
    return Determinant( Matrix( [ElementToSequence(OL ! x) : x in w] ) );
end intrinsic;


/******************************** AlgAssAction ***************************************
 *
 * S is a Rec AmodpModule, lambda an element in Amodp, a an element in S.
 *
 * Computes lambda(a), which is an element in S.
 *
 *******************************************************************************************/

intrinsic AlgAssAction(S :: Rec, lambda :: AlgAssElt, a :: ModTupFldElt) -> ModTupFldElt
{}
    local B, basis, coeff, i;

    B := Parent(lambda);
    basis := Basis(B);
    coeff := ElementToSequence(lambda);
    return &+[ coeff[i] * a * Transpose(S`phi(basis[i])) : i in [1..#basis] ];
end intrinsic;

/******************************** AlgAssSubModule ***************************************
 *
 * S is a Rec AmodpModule, Amodp the acting AlgAss, J an ideal in Amodp (e.g. the 
 * Jacobson radical).
 *
 * Computes the submodule J*S.
 *
 *******************************************************************************************/

intrinsic AlgAssSubModule(S :: Rec, Amodp :: AlgAss, J :: AlgAss) -> Rec
{}
    local gens, lambda, a,  Q;


    gens := [ AlgAssAction(S, Amodp ! lambda, a) : lambda in Basis(J), a in Basis(S`M) ];

    Q := rec<AmodpModule | >;
    Q`M := sub<S`M | gens >;
    Q`phi := S`phi;
    return Q;
end intrinsic;

/******************************** AlgAssSubModule ***************************************
 *
 * S is a Rec AmodpModule, R the acting AlgAss, gens a sequence of elements of S.
 *
 * Computes the submodule generated by gens over R.
 *
 *******************************************************************************************/

intrinsic AlgAssSubModule(S :: Rec, R :: AlgAss, gens :: SeqEnum) -> Rec
{}
    local L, x, b, Q;

    L := [];
    for x in gens do
        for b in Basis(R) do
            Append(~L, AlgAssAction(S, b , x));
        end for;
    end for;
    Q := rec<AmodpModule | >;
    Q`M := sub<S`M | L >;
    Q`phi := S`phi;
    return Q;
end intrinsic;

/******************************** FSubModule2 ***************************************
 *
 * S is a Rec AmodpModule, gens a sequence of elements of S.
 *
 * Computes the F-submodule generated by gens.
 * Attention: This is not necessarily a Amodp-module, although we return a Rec AmodpModule.
 *
 *******************************************************************************************/

intrinsic FSubModule2(S :: Rec, gens :: SeqEnum) -> Rec
{}
    local L, Q;

    L := [];
    Q := rec<AmodpModule | >;
    Q`M := sub<S`M | gens >;
    Q`phi := S`phi;
    return Q;
end intrinsic;


/******************************** AtoAmodpModule ***************************************
 *
 * Amodp is an AlgAss which comes from an order A by reduction mod p, M is an A-module,
 * given as a Z-module.
 *
 * Computes the Amodp-module M / pM as a Rec AmodpModule.
 *
 *******************************************************************************************/

intrinsic AtoAmodpModule(Amodp :: AlgAss, M :: Rec, p :: RngIntElt) -> Rec
{}
    local n, D, T, BasisM, Bs, b, lambda, B;


    n := Dimension(Codomain(M`phi));

    D := rec<AmodpModule | >;
    D`M := RSpace(GF(p), n);
    T := MatrixRing(GF(p), n);
    BasisM := [ Vector(Rationals(), v) : v in RowSequence(M`hnf) ];
    Bs := [];
    for b in Basis(Amodp) do
        lambda := Amodp`lift(b);
        B := Matrix([ElementToSequence( GroupRingAction(M, lambda, c)*M`inv ) : c in BasisM ]);
        Append(~Bs, b->MakeMatrixIntegral(Transpose(B)));
    end for;

    // D`phi := map<Basis(D`M)->T | Bs >;
    D`phi := map<Basis(Amodp)->T | Bs >;
    D`lift := map<Basis(D`M)->BasisM | [ Basis(D`M)[i]->BasisM[i] : i in [1..#Basis(D`M)]]>;

    return D;
end intrinsic;


/******************************** LiftToM ***************************************
 *
 * Mp is a Rec AmodpModule which comes from an A-module M by reduction mod p, 
 * v is an element in Mp`M.
 * given as a Z-module.
 *
 * Computes the lift of v as an element in M.
 *
 *******************************************************************************************/

intrinsic LiftToM(Mp :: Rec, v :: ModTupFldElt) -> ModTupRngElt
{}
    local i;

    return &+[Integers()!v[i] * Mp`lift(Basis(Mp`M)[i]) : i in [1..#Basis(Mp`M)]];

end intrinsic;


/******************************** AlgAssQuotientModule ***************************************
 *
 * S is a Rec AmodpModule, U is a Rec AmodpModule of the form J*S, 
 * psi is a map Amodp -> R:= Amodp / J.
 *
 * Computes the R-Module S/U and the map S -> S/U.
 *
 *******************************************************************************************/

intrinsic AlgAssQuotientModule(S :: Rec, U :: Rec, psi :: Map) -> Rec, Map
{}
    local F, Q, T, Amodp, R, BasisQ, Bs, b, lambda, B, c,tau;

    F := Parent( S`phi(Basis(S`M)[1])[1,1] );

    Q := rec<AmodpModule | >;
    Q`M, tau := quo<S`M | U`M >;
    T := MatrixRing(F, Dimension(Q`M));

    Amodp := Domain(psi);
    R := Codomain(psi);

    BasisQ := [ (tau^-1)(b) : b in Basis(Q`M) ];
    Bs := [];
    for b in Basis(R) do
        lambda := Amodp ! (psi^-1)(b);
        B := Matrix([ElementToSequence( tau( AlgAssAction(S, lambda, c) ) ) : c in BasisQ ]);
        Append(~Bs, b->Transpose(B));
    end for;

    Q`phi := map<Basis(R)->T | Bs >;

    return Q, tau;
end intrinsic;


/******************************** FindFBasis ***************************************
 *
 * V is a Rec AmodpModule, actually a module over the matrix algebra M = M_n(F). V is given
 * by a basis over GF(p). F is a finite extension of GF(p). Q is a list needed to compute
 *		beta : M --> component of R
 * R is a semisimple algebra over GF(p) with a component isomorphic to M.
 *
 * Computes a F-basis of V.
 *
 *******************************************************************************************/

intrinsic FindFBasis(V :: Rec, M :: AlgMat, Q :: List, R :: AlgAss) -> SeqEnum
{}

    local F, n, Fbasis, Vbasis, BasisOverF, W, b, d, i, WN, N;

    F := CoefficientRing(M);
    n := Integers() ! (Dimension(V`M) / Dimension(F));
    Fbasis := Basis(F);
    Vbasis := Basis(V`M);
    BasisOverF := [ Vbasis[1] ];
    W := Matrix( [ AlgAssAction(V, (R!beta(M!b, Q)), Vbasis[1]) : b in Fbasis ] );

    d := 1;
    i := 2;
    while d lt n do
        N := Matrix( [ AlgAssAction(V, (R!beta(M!b, Q)), Vbasis[i]) : b in Fbasis ] );
        WN := VerticalJoin(W, N);
        if Dimension(Kernel(WN)) eq 0 then
            W := WN;
            Append(~BasisOverF, Vbasis[i]);
            d +:= 1;
        end if;
        i +:= 1;
    end while;
    return BasisOverF;

end intrinsic;

/******************************** MoritaBasis ***************************************
 *
 * M is a matrix algebra M_n(F), which is assumed to be isomorphic to a component of the
 * semisimple GF(p)-algebra R. Here F is a finite extension of GF(p).
 * P is a M-module, and Q is a list needed to compute
 *              beta : M --> component of R
 *
 * If P is free over M, then the function computes a M-basis of P and returns true and the basis. 
 * Otherwise returns false, [].
 *
 *******************************************************************************************/

intrinsic MoritaBasis(M :: AlgMat, Q :: List, P :: Rec,  R :: AlgAss) -> BoolElt, List
{}

    local n, E, gens, V, b, F, FBasisV, i, l, w, rk;


    n := Degree(M);
    F := CoefficientRing(M);
    rk := Dimension(P`M) / (Dimension(F) * n^2);
    if not IsIntegral(rk) then
        return false, [];
    end if;
    rk := Integers() ! rk;

    E := Matrix(n, n, [(R ! beta(b, Q)) : b in Basis(M)]);
    gens := [ AlgAssAction(P, E[1,1], b) : b in Basis(P`M) ];
    V := FSubModule2(P, gens);
    FBasisV := FindFBasis(V, M, Q, R);

    if #FBasisV ne rk*n then
        return false, [];
    end if;

    i := 0; l := #FBasisV;
    w := [];
    while i lt #FBasisV do
        Append(~w, &+[ AlgAssAction(P, E[j,1], FBasisV[i+j])  : j in [1..Min(n, l)] ]);
        l := l-n; i := i+n;
    end while;

    return true, w;
end intrinsic;


/******************************** IsLocallyFree ***************************************
 *
 * QG is a groupring over the rationals, A is a Z-order in QG, M is a A-module and p a 
 * rational prime.
 *
 * If M is locally free over A at p, then the function computes a A_p-basis of M_p and 
 * returns true and the basis. Otherwise returns false, [].
 *
 *******************************************************************************************/

intrinsic IsLocallyFree(QG :: AlgGrp, A :: Rec, M :: Rec, p :: RngIntElt) -> BoolElt, List
{}

    local F, Abasis, w, T, b1, b2, Amodp, i, j, J, Mp, R, psi, JMp, MpmodJ, tau, idem, comp, C,
          e, gens, B, MatM, Q, v, isfree; 

    F := GF(p);
    Abasis := [ QG ! w : w in RowSequence(A`hnf) ];
    T := [ ElementToSequence(b2*b1) : b1 in Abasis, b2 in Abasis ];
    T := Matrix(T) * A`inv;
    Amodp := AssociativeAlgebra< F, Nrows(A`hnf) | ElementToSequence(T) : Check := false >;
    AddAttribute(AlgAss, "lift");
    Amodp`lift := map< Basis(Amodp) -> Abasis | [Basis(Amodp)[i]->Abasis[i] : i in [1..#Abasis] ] >;
    J := JacobsonRadical(Amodp);

    Mp := AtoAmodpModule(Amodp, M, p);
    R, psi := quo<Amodp | J>;
    // print "Compute JMp";
    JMp := AlgAssSubModule(Mp, Amodp, J);
    // print "Compute DpmodJ";
    MpmodJ, tau := AlgAssQuotientModule(Mp, JMp, psi);

    if Dimension(R) ne Dimension(MpmodJ`M) then
        return false, [];
    end if;

    idem, comp := CentralIdempotents(R);
    C := [];
    for i:=1 to #comp do
        // A := comp[i];
        // print "i = ", i, "\nA = ", A;
        // e := (psi^-1)(R ! idem[i]);
        e := idem[i];
        gens := [ AlgAssAction(MpmodJ, e, b) : b in Basis(MpmodJ`M) ];
        // print "#gens = ", #gens;
        Append(~C, AlgAssSubModule(MpmodJ, R, gens));
    end for;

    B := [* *];
    for i:=1 to #comp do
        /* Berechne die Basis in der i-ten Komponente */
        MatM, Q := CentralSimpleAlgRep(comp[i]);
        isfree, v := MoritaBasis(MatM, Q, C[i], R);
        if not isfree then
            return false, [];
        end if;
        Append(~B, v);
    end for;

    /* In jeder Komponente C[i] ist nun eine comp[i]-Basis berechnet. Diese Basen sind nun
       zusammenzusetzen.
    */

    w := [ (tau^-1)( &+[B[i][j] : i in [1..#B]] )  : j in [1..#B[1]] ];
    // print "w = ", w;
    return true, [ LiftToM(Mp, v) : v in w ];
end intrinsic;


/******************************** FindOFBasis ***************************************
 *
 * QG is a groupring over the rationals. We assume that the components of QG are known and
 * stored in QG`W. ind specifies a component of QG, namely QG`W[ind] = M_n(F), where F is a
 * finite extension of Q. V is a free M_n(OF)-module.
 *
 * Computes an OF-basis of V.
 * Attention: We always assume that F has class number 1.
 *
 *******************************************************************************************/

intrinsic FindOFBasis(QG :: AlgGrp, ind :: RngIntElt, V :: Rec) -> SeqEnum
{}

    local F, n, Fbasis, Vbasis, BasisOverF, J, FbasisInQG, lambda, W, b, d, WN, N, i,
          OF, A, s, I, c, BasisOverOF, BasisS, M, S;

    F := CoefficientRing(QG`W[ind]);
    assert ClassNumber(F) eq 1;

    n := Integers() ! (Nrows(V`hnf) / Degree(F));
    Fbasis := Basis(F);
    Vbasis := [ Vector(Rationals(), c) : c in RowSequence(V`hnf) ];
    BasisOverF := [ Vbasis[1] ];
    
    J := [* Zero(w) : w in QG`W *];
    FbasisInQG := [];
    for b in Fbasis do
            J[ind] := ScalarMatrix(QG`W[ind], b);
            lambda := WedIsoInv(QG, J);
            Append(~FbasisInQG, lambda);
    end for;

    W := Matrix( [ GroupRingAction(V, lambda, Vbasis[1]) : lambda in FbasisInQG ] );

    d := 1;
    i := 2;
    while d lt n do
        N := Matrix( [ GroupRingAction(V, lambda, Vbasis[i]) : lambda in FbasisInQG ] );
        WN := VerticalJoin(W, N);
        if Dimension(Kernel(WN)) eq 0 then
            W := WN;
            Append(~BasisOverF, Vbasis[i]);
            d +:= 1;
        end if;
        i +:= 1;
    end while;

    OF := MaximalOrder(F);
    V_Qbasis := [];
    for v in BasisOverF do
        w := [ GroupRingAction(V, lambda, v) : lambda in FbasisInQG ];
        V_Qbasis cat:= w ;
    end for;

    s := #Fbasis; r := #BasisOverF;
    A := Matrix(V_Qbasis);
    b := Solution(A, Vbasis);
    c := [ [ F ! [x[k] : k in [(j-1)*s+1..j*s] ] : j in [1..r] ] : x in b ];
    d := Denominator( ideal<OF | &cat[x : x in c]> );
    J := [* Zero(w) : w in QG`W *];
    d := Denominator( ideal<OF | &cat[x : x in c]> );
    c := [ [OF!(d*y) : y in x] : x in c ]; 
    M := Module(OF, r);
    gens := [M!x : x in c];
    S := sub< M |gens >;
    BasisS := PseudoBasis(S);

    BasisOverOF := [];
    for i:= 1 to #BasisS do
        v := [];
        for j:=1 to #BasisOverF do
	    J[ind] := ScalarMatrix(QG`W[ind], BasisS[i][2][j]);
	    lambda :=  WedIsoInv(QG, J);
	    Append(~v, GroupRingAction(V, lambda, BasisOverF[j]) );
	end for;
	gen := &+v;

	I := ideal<OF | [OF!x : x in Generators(BasisS[i][1])]>;
	isprincipal, alpha := IsPrincipal(I);
	J[ind] := ScalarMatrix(QG`W[ind], alpha);
	lambda := WedIsoInv(QG, J);
	Append(~BasisOverOF, GroupRingAction(V, 1/d * lambda, gen));
    end for;
    
    return BasisOverOF;

end intrinsic;


/******************************** VeryNiceMoritaBasis ***************************************
 *
 * QG is a groupring over the rationals. We assume that the components of QG are known and
 * stored in QG`W. ind specifies a component of QG, namely QG`W[ind] = M_n(F), where F is a
 * finite extension of Q. P is a M_n(OF)-module.
 *
 * Computes an M_n(OF)-basis of V, if one exists.
 * Attention: We always assume that F has class number 1.
 *
 *******************************************************************************************/

intrinsic VeryNiceMoritaBasis(QG :: AlgGrp, ind :: RngIntElt, P :: Rec) -> BoolElt, List
{}

    local n, E, gens, V, b, F, FBasisV, i, l, w, rk, J, BasisW, k, BasisP;


    n := Degree(QG`W[ind]);
    F := CoefficientRing(QG`W[ind]);
    rk := Nrows(P`hnf) / (Degree(F) * n^2);
    if not IsIntegral(rk) then
        return false, [];
    end if;
    rk := Integers() ! rk;

    J := [* Zero(w) : w in QG`W *]; 
    E := ZeroMatrix(QG, n, n);
    BasisW := [ w : w in Basis(QG`W[ind])];
    for k := 1 to n do
        for l:=1 to n do
	    J[ind] := BasisW[(k-1)*n+l];
	    E[k, l] := WedIsoInv(QG, J);
	end for;
    end for;
	    
    BasisP := [Vector(Rationals(),c) : c in RowSequence(P`hnf)];	    
    gens := [ GroupRingAction(P, E[1,1], b) : b in BasisP ];
    V := ZGSubModule(P, gens);
    OFBasisV := FindOFBasis(QG, ind, V);

    if #OFBasisV ne rk*n then
        return false, [];
    end if;

    i := 0; l := #OFBasisV;
    w := [];
    while i lt #OFBasisV do
        Append(~w, &+[ GroupRingAction(P, E[j,1], OFBasisV[i+j])  : j in [1..Min(n, l)] ]);
        l := l-n; i := i+n;
    end while;

    return true, w;
end intrinsic;


/******************************** VeryNiceBasis ***************************************
 *
 * QG is a groupring over the rationals. We assume that the components of QG are known and
 * stored in QG`W. M is a module over the very nice maximal order 
 *                 MaxOrd = \bigoplus M_{n_i}(OF_i).
 *
 * Computes MaxOrd-basis of M, if one exists.
 * Attention: We always that all fields F_i have class number 1.
 * This is not necessary, the implementation just needs additional effort!!  
 *
 *******************************************************************************************/

intrinsic VeryNiceBasis(QG :: AlgGrp, M :: Rec) -> BoolElt, List
{}

    local F,  w, T, i, j, J, R, psi,  C,
          e, gens, B, MatM, Q, v, isfree;

    // E := QIdempotents(QG, QG`X, QG`H);
    E := QIdempotents(QG);
    BasisM := [Vector(Rationals(), c) : c in RowSequence(M`hnf)];
    C := [];
    for i:=1 to #E do
        gens := [ GroupRingAction(M, E[i], b) : b in BasisM ];
        Append(~C, ZGSubModule(M, gens));
    end for;

    B := [* *];
    for i:=1 to #QG`W do
        print i, "-te Komponente";
        /* Berechne die Basis in der i-ten Komponente */
        isfree, v := VeryNiceMoritaBasis(QG, i, C[i]);
        if not isfree then
            return false, [];
        end if;
        Append(~B, v);
    end for;

    /* In jeder Komponente C[i] ist nun eine comp[i]-Basis berechnet. Diese Basen sind nun
       zusammenzusetzen.
    */
    w := [ &+[B[i][j] : i in [1..#B]]   : j in [1..#B[1]] ];
    return true, w;
end intrinsic;


/******************************** TestLocalBasis ***************************************
 *
 * Checks whether the basis computed by IsLocallyFree is actually a AssOrd basis at p.
 * Only used for testing.
 *
 *******************************************************************************************/

intrinsic TestLocalBasis(QG :: AlgGrp, AssOrd :: Rec, P :: Rec, basis :: SeqEnum, p :: RngIntElt) -> BoolElt
{}

    local AssOrdBasis, gens, W, v, b, C, ind;

    AssOrdBasis := [ QG ! b : b in RowSequence(AssOrd`hnf) ];
    gens := [ ElementToSequence(GroupRingAction(P, b, v)) : v in basis, b in AssOrdBasis ];
    W := Matrix(Rationals(), gens);
    C := ZGModuleInit(W, P`phi);
    ind := ZGModuleIndex(P, C);
    print "ind = ", ind;
    return not (IsDivisibleBy(Numerator(ind), p) or IsDivisibleBy(Denominator(ind), p));
end intrinsic;
 

/******************************** TestBasisOverMaxOrd ***************************************
 *
 * Checks whether the basis computed by VeryNiceBasis is actually a MaxOrd basis.
 * Only used for testing.
 *
 *******************************************************************************************/

intrinsic TestBasisOverMaxOrd(QG :: AlgGrp, MaxOrd :: Rec, P :: Rec, basis :: SeqEnum) -> RngIntElt
{}

    local MaxOrdBasis, gens, W, v, b, C, ind;

    MaxOrdBasis := [ QG ! b : b in RowSequence(MaxOrd`hnf) ];
    gens := [ ElementToSequence(GroupRingAction(P, b, v)) : v in basis, b in MaxOrdBasis ];
    W := Matrix(Rationals(), gens);
    C := ZGModuleInit(W, P`phi);
    ind := ZGModuleIndex(P, C);
    print "ZGModuleCompare:  ", ZGModuleCompare(P, C);
    return ind;
end intrinsic;


/******************************** S3_rep ***************************************
 *
 * No longer needed. Use Dn_rep instead.
 *
 *******************************************************************************************/

intrinsic S3_rep(QG :: AlgGrp) -> List, SeqEnum
{}
    local Gl, G, chi, iota, sigma, tau, Sigma, Tau, M, N;

    Gl := [ GL( Integers()!Degree(QG`X[i][1]), QG`H[i][1] ) : i in [1..#QG`X]];
    W := [ MatrixRing( QG`H[i][1], Integers()!Degree(QG`X[i][1]) ) : i in [1..#QG`X]];

    G := Group(QG);

    chi := [QG`X[i][1] : i in [1..#QG`X]];
    iota := [QG`H[i][3] : i in [1..#QG`X]];
    sigma := G.1;
    tau := G.2;
    Sigma := [* Gl[1]![iota[1](chi[1](sigma))], Gl[2]![iota[2](chi[2](sigma))] *];
    Tau := [* Gl[1]![iota[1](chi[1](tau))], Gl[2]![iota[2](chi[2](tau))] *];
    M := Matrix([ [iota[3](0), iota[3](1)], [iota[3](-1), iota[3](-1)] ]);
    N := Matrix([ [iota[3](0), iota[3](1)], [iota[3](1), iota[3](0)] ]);
    Append(~Sigma, Gl[3] ! M);
    Append(~Tau, Gl[3] ! N);

    return [* hom<G -> Gl[i] | [sigma->Sigma[i], tau->Tau[i]]> : i in [1..#QG`X] *], W;
end intrinsic;


/******************************** A4_rep ***************************************
 *
 * Computes the representations of G, which is assumed to be isomorphic to A4.
 *
 *******************************************************************************************/

intrinsic A4_rep(QG :: AlgGrp) -> List, SeqEnum
{}
    local Gl, G, chi, iota, sigma, nu, Sigma, Nu, M, N, A4, e, n, T, Omega, X, H,
          b, iso;

    G := Group(QG);

    A4 := AlternatingGroup(4);
    e := Exponent(A4); n := #A4;
    T := CharacterTable(A4);
    Omega := CompOmega(A4);      /* Omega = Galois group of Q(zeta_e) / Q */
    X := OmegaRepresentatives(T, Omega);
    H := ComputeKi(G, X);
    Gl := [ GL( Integers()!Degree(X[i][1]), H[i][1] ) : i in [1..#X]];
    W := [ MatrixRing( H[i][1], Integers()!Degree(X[i][1]) ) : i in [1..#X]];
    chi := [X[i][1] : i in [1..#X]];
    iota := [H[i][3] : i in [1..#X]];
    sigma := A4.1;
    nu := A4.2;
    Sigma := [* Gl[1]![iota[1](chi[1](sigma))], Gl[2]![iota[2](chi[2](sigma))] *];
    Nu := [* Gl[1]![iota[1](chi[1](nu))], Gl[2]![iota[2](chi[2](nu))] *];
    M := Matrix([ [iota[3](1), iota[3](0), iota[3](0)], 
                  [iota[3](0), iota[3](-1), iota[3](0)],
                  [iota[3](0), iota[3](0), iota[3](-1)] 
                ]);
    N := Matrix([ [iota[3](0), iota[3](1), iota[3](0)],
                  [iota[3](0), iota[3](0), iota[3](1)],
                  [iota[3](1), iota[3](0), iota[3](0)]
                ]);
    Append(~Sigma, Gl[3] ! M);
    Append(~Nu, Gl[3] ! N);

    b, iso := IsIsomorphic(G, A4);
    return [* iso * hom<A4 -> Gl[i] | [sigma->Sigma[i], nu->Nu[i]]> : i in [1..#X] *], W;
end intrinsic;

/******************************** S4_rep ***************************************
 *
 * Computes the representations of G, which is assumed to be isomorphic to S4.
 *
 *******************************************************************************************/

intrinsic S4_rep(G :: GrpPerm) -> List, SeqEnum
{}
    local i, S4, x,y,z,s,t,S,T,X,Y,V4,Q,pi,U,d,dt,S3,sigma,tau,e,n,Omega,H,Gl,W,
          iota, Sigma,Tau,M,N,b,iso,Reps;

    S4 := SymmetricGroup(4);
    x := S4 ! [2,1,4,3];
    y := S4 ! [3,4,1,2];
    z := S4 ! [4,3,2,1];
    s := S4 ! [2,3,1,4];
    t := S4 ! [2,1,3,4];

    V4 := sub< S4 | [x,y,z] >;
    Q, pi := quo< S4 | V4 >;

    S3 := SymmetricGroup(3);
    sigma := S3.1; tau := S3.2;
    
    e := Exponent(S4); n := #S4;
    T := CharacterTable(S4);
    Omega := CompOmega(S4);      /* Omega = Galois group of Q(zeta_e) / Q */
    X := OmegaRepresentatives(T, Omega);
    H := ComputeKi(G, X);
    Gl := [ GL( Integers()!Degree(X[i][1]), H[i][1] ) : i in [1..#X]];
    W := [ MatrixRing( H[i][1], Integers()!Degree(X[i][1]) ) : i in [1..#X]];

    iota := [H[i][3] : i in [1..#X]];
    Sigma := [* Gl[1]![iota[1](1)], Gl[2]![iota[2](1)] *];
    Tau := [* Gl[1]![iota[1](1)], Gl[2]![iota[2](-1)] *];
    M := Matrix([ [iota[3](0), iota[3](-1)],
                  [iota[3](1), iota[3](-1)]
                ]);
    N := Matrix([ [iota[3](1), iota[3](-1)],
                  [iota[3](0), iota[3](-1)]
                ]);
    Append(~Sigma, Gl[3] ! M);
    Append(~Tau, Gl[3] ! N);

    b, iso := IsIsomorphic(Q, S3);
    Reps := [* pi*iso*hom<S3 -> Gl[i] | [sigma->Sigma[i], tau->Tau[i]]> : i in [1,2,3] *];

    CharTab := CharacterTable(V4);
    chi := [CharTab[i] : i in [2,3,4]];
    S := Gl[4] ! Matrix([ [0,1,0], [0,0,1],[1,0,0] ]);
    T := Gl[4] ! Matrix([ [0,0,1], [0,1,0],[1,0,0] ]);
    X := Gl[4] ! Matrix([ [chi[1](x), 0, 0], [0, chi[2](x), 0], [0,0,chi[3](x)] ]);
    Y := Gl[4] ! Matrix([ [chi[1](y), 0, 0], [0, chi[2](y), 0], [0,0,chi[3](y)] ]);
    U := sub<S4 | [s,t,x,y]>;
    b, iso := IsIsomorphic(S4, U);
    d := iso*hom<U->Gl[4] | [s->S, t->T, x->X, y->Y]>;

    S := ScalarMatrix(3, Reps[2](s)[1,1] ) * S;
    T := ScalarMatrix(3, Reps[2](t)[1,1] ) * T;
    X := ScalarMatrix(3, Reps[2](x)[1,1] ) * X;
    Y := ScalarMatrix(3, Reps[2](y)[1,1] ) * Y;
    dt := iso*hom<U->Gl[4] | [s->S, t->T, x->X, y->Y]>;
    
    Reps cat:= [* d, dt *];
    b, iso := IsIsomorphic(G, S4);
    Reps := [* hom<G->Codomain(d) | [x->d(iso(x)) : x in Generators(G)]> : d in Reps *];

    return Reps, W;
end intrinsic;

/******************************** Dn_rep ***************************************
 *
 * Computes the representations of G, which is assumed to be isomorphic to Dn, the
 * dihedral group of order 2n.
 *
 *******************************************************************************************/

intrinsic Dn_rep(QG :: AlgGrp) -> List, SeqEnum
{}
    local Gl, G, chi, iota, sigma, tau, Sigma, Tau, M, N, Dn, e, n, T, Omega, X, H,
          b, iso;

    G := Group(QG);

    Dn := DihedralGroup(Integers() ! (#G / 2));
    e := Exponent(Dn); n := #Dn;
    T := CharacterTable(Dn);
    Omega := CompOmega(Dn);      /* Omega = Galois group of Q(zeta_e) / Q */
    X := OmegaRepresentatives(T, Omega);
    H := ComputeKi(Dn, X);
    Gl := [ GL( Integers()!Degree(X[i][1]), H[i][1] ) : i in [1..#X]];
    W := [ MatrixRing( H[i][1], Integers()!Degree(X[i][1]) ) : i in [1..#X]];
    chi := [X[i][1] : i in [1..#X]];
    iota := [H[i][3] : i in [1..#X]];
    sigma := Dn.1;
    tau := Dn.2;
    Sigma := [* *];  Tau := [* *];
    for i:=1 to #chi do
        if Degree(chi[i]) eq 1 then
	    Append( ~Sigma, Gl[i] ! [iota[i](chi[i](sigma))] );
	    Append( ~Tau, Gl[i] ! [iota[i](chi[i](tau))] );
	end if;
	if Degree(chi[i]) eq 2 then
	    M := Matrix([ [iota[i](0), iota[i](-1)],
	                  [iota[i](1), iota[i](chi[i](sigma))]
			]);
	    N := Matrix([ [iota[i](1), iota[i](chi[i](sigma))],
                          [iota[i](0), iota[i](-1)]
                        ]);
	    Append( ~Sigma, Gl[i] ! M );
            Append( ~Tau, Gl[i] ! N );
        end if;
    end for;

    b, iso := IsIsomorphic(G, Dn);
    return [* iso * hom<Dn -> Gl[i] | [sigma->Sigma[i], tau->Tau[i]]> : i in [1..#X] *], W;
end intrinsic;


/******************************** Abelian_rep ***************************************
 *
 * Computes the representations of G, which is assumed to be abelian.
 *
 *******************************************************************************************/

intrinsic Abelian_rep(G :: GrpPerm) -> List, SeqEnum
{}
    local Gl, chi, iota, T, Omega, X, H;

    T := CharacterTable(G);
    Omega := CompOmega(G);      /* Omega = Galois group of Q(zeta_e) / Q */
    X := OmegaRepresentatives(T, Omega);
    H := ComputeKi(G, X);
    Gl := [ GL( Integers()!Degree(X[i][1]), H[i][1] ) : i in [1..#X]];
    W := [ MatrixRing( H[i][1], Integers()!Degree(X[i][1]) ) : i in [1..#X]];
    chi := [X[i][1] : i in [1..#X]];
    iota := [H[i][3] : i in [1..#X]];
    homs := [* hom<G -> Gl[i] | [g->Gl[i]![iota[i](chi[i](g))] : g in Generators(G)]> : i in [1..#X] *];
    return homs, W;
end intrinsic;



/******************************** ComputePhiMatrix ***************************************
 *
 * QG is a groupring over the rationals, rep the list of irreducible QG-representations and
 * W the list of corresponding matrix rings of the Wedderburn decomposition. We always assume
 * that all Schur indices (over Q) are 1.
 *
 * Computes a matrix which describes the isomprphism
 *		Phi : QG --> Wedderburn decomposition
 * with respect to naturally chosen basis.
 *
 *******************************************************************************************/

intrinsic ComputePhiMatrix(QG :: AlgGrp, rep :: List, W :: SeqEnum) -> Mtrx
{}
    local Z, g, c, i, D, k, l, Phi;

    Z := [];

    for g in Group(QG) do
        c := [];
	for i := 1 to #W do
	   D := rep[i](g);
	   for k:=1 to Nrows(D) do
	       for l:=1 to Ncols(D) do
		   c cat:= ElementToSequence(D[k,l]);
	       end for;
	   end for;
	end for;
	Append(~Z, c);
	// Z cat:= c;
    end for;
    // Phi := Transpose( Matrix(Z) );
    Phi := Matrix(Z) ;

    return Phi;
end intrinsic;


/******************************** WedIso ***************************************
 *
 * lambda is an element in QG. We assume that QG`W, QG`Phi and QG`PhiInv are initialized.
 *
 * Computed Phi(lambda), where 
 *              Phi : QG --> Wedderburn decomposition.
 *
 *******************************************************************************************/

intrinsic WedIso(lambda :: AlgGrpElt) -> List
{}
    local QG, c, T, first, i, n, K, length, D, k, l;

    QG := Parent(lambda);
    // c := Vector( ElementToSequence(lambda) ) * Transpose(QG`Phi);
    c := Vector( ElementToSequence(lambda) ) * QG`Phi;

    T := [* *];
    first := 1;
    for i:=1 to #QG`W do
        n := Degree(QG`W[i]);
	K := CoefficientRing(QG`W[i]);
	length := Degree(K);
	D := ZeroMatrix(K, n, n);
	for k:=1 to n do
	    for l:=1 to n do
	       D[k,l] := K ! [c[first+t] : t in [0..length-1]];
	       first := first + length;
	    end for;
	end for;
	Append(~T, D);
    end for;
    return T;
end intrinsic;
	
/******************************** WedIsoInv ***************************************
 *
 * We assume that QG`W, QG`Phi and QG`PhiInv are initialized. T is list of matrices in
 * \bigoplus QG`W[i].
 *
 * Computed Phi^-1(T), where
 *              Phi : QG --> Wedderburn decomposition.
 *
 *******************************************************************************************/
    
intrinsic WedIsoInv(QG :: AlgGrp, T :: List) -> AlgGrpElt
{}
    local c, i, k, l, lambda;

    c := [];
    for i:=1 to #QG`W do
        for k:=1 to Nrows(T[i]) do
	    for l:=1 to Ncols(T[i]) do
	       c cat:= ElementToSequence(T[i][k, l]);
	    end for;
	end for;
    end for;
    lambda := QG ! ElementToSequence( Vector(c)*QG`PhiInv );
    return lambda;
end intrinsic;


/******************************** xtoi ***************************************
 *
 * Only used in A4-Enumeration.
 *
 *******************************************************************************************/

intrinsic xtoi(x :: ModMatRngElt, m :: RngIntElt) -> RngIntElt
{}
    local i;

    // return Integers()!x[2,1]*m + Integers()!x[1,1] + 1;
    return &+[ Integers()!x[i,1]*m^(i-1) : i in [1..Nrows(x)] ] + 1;

end intrinsic;

/********************************  ElementaryMatrices ***************************************
 *
 * Only used in A4-Enumeration. Should be repleced by GlGenerators.
 *
 *******************************************************************************************/

intrinsic ElementaryMatrices(n :: RngIntElt) -> SeqEnum
{}
    local E, i, j;

    E := [];
    for i:=1 to n do
        for j:=1 to n do
            if i ne j then
                A := ScalarMatrix(n, 1);
                A[i,j] := 1;
                Append(~E, A);
            end if;
        end for;
    end for;
    for i:=1 to n do
        A := ScalarMatrix(n, 1);
        A[i,i] := -1;
        Append(~E, A);
    end for;
    Append(~E, ScalarMatrix(n, -1));

    return E;
end intrinsic;

/********************************  GlGenerators ***************************************
 *
 * M is a matrix ring M_n(F), F a finite extension of Q.  Suppose frg is an integral ideal of OF.
 *
 * Then this function computes preimages of generators of the image of Gl_n(OF) in Gl_n(OF/frg).
 *
 *******************************************************************************************/

intrinsic GlGenerators(M :: AlgMat) -> SeqEnum
{}
    local E, n, F, OF, U, f, A, b, i, j;

    E := [];
    n := Degree(M);
    F := CoefficientRing(M);
    OF := MaximalOrder(F);
    U, f := UnitGroup(OF);
    for u in Generators(U) do
        A := ScalarMatrix(n, OF!1);
        A[1,1] := f(u);
        Append(~E, A);
    end for;

    for i:=1 to n do
        for j:=1 to n do
            if i ne j then
	        for b in Basis(OF) do
                    A := ScalarMatrix(n, OF!1);
                    A[i,j] := b;
                    Append(~E, A);
		end for;
            end if;
        end for;
    end for;

    return E;
end intrinsic;

/********************************  Enumeration ***************************************
 *
 * QG is a group ring with QG`W etc.initialized. 
 * Ideals is a list of ideals (the central conductors), wnb is a basis over the maximal order
 * and Ath is a free AssOrd-module. Tries to find a generator over AssOrd by Enumeration.
 *
 *******************************************************************************************/

intrinsic Enumeration(QG :: AlgGrp, Ideals :: List, wnb :: ModTupFldElt, Ath :: Rec) -> List
{}
    local n, m, a, b,R, M, X, S, E, Q, L, W, rho, phi, F, BasisF, psi, K, U, pi,
          U1, U2, u1, u2, v, lambda, beta, T;

    U := [* *];
    for i:=1 to #Ideals do
        n := Degree(QG`W[i]);
	OF := MaximalOrder( CoefficientRing(QG`W[i]) );
	R, f := quo< OF | Ideals[i] >;
	B := BasisMatrix(Ideals[i]);
	T := CartesianProduct( [ [1..Abs(B[i,i])] : i in [1..Nrows(B)] ] );
	EltsOfR := [ f( OF ! [x : x in t] ) : t in T ];
        M := MatrixRing(R, n);
	X := CartesianPower(EltsOfR, n);
        X := [ Matrix(n, 1, [y: y in x]) : x in X ];
        print "i = ", i, "    #X = ", #X;
	S := Sym(#X);
	E := GlGenerators(QG`W[i]);
        Q := [ S ! [Index(X, M!g*x) : x in X] : g in E ];
        L := sub< S | Q >;
        W := WordGroup(L);
        rho :=InverseWordMap(L);
        W := Codomain(rho);
        phi := hom< W->L | Q >;
        F := FreeGroup(#Generators(W));
        BasisF := [F.i : i in [1..#Generators(F)]];
        psi := hom< W->F | BasisF >;
	V := [E[1]^0];
	print "#L = ", #L;
	counter := 1;
	for y in L do
	    if counter mod 1000 eq 0 then print "  ", counter; end if;
	    v := ElementToSequence(psi(rho(y)));
	    if #v gt 0 then
	        w := [];
	        for i in v do
		    if i ge 0 then
		        Append(~w, E[i]);
		    else
		        Append(~w, E[-i]^-1);
		    end if;
		end for;
		T := &*w;
	        // T := &*[E[i] : i in v];
		Append(~V, T);
	    end if;
	    counter := counter+1;
	end for;
	Append(~U, V);
    end for;

    print "Start to enumerate: Size of set ----> ", &*[#V : V in U]; 
    counter := 1;
    EnumSet := CartesianProduct([ [1..#U[i]] : i in [1..#U] ]);
    for y in EnumSet do
        lambda := WedIsoInv(QG, [* U[i][y[i]] : i in [1..#U] *]);
        if counter mod 100 eq 0 then printf "%o  ", counter; end if;
        beta := GroupRingAction(Ath, lambda, wnb);
        if IsZGModuleElt(ElementToSequence(beta), Ath) then
             print "HURRA HURRA***************************************************";
             return beta;
        end if;
        counter := counter+1;
    end for;
    return beta;
end intrinsic;


/********************************  S4_Enumeration ***************************************
 *
 * QG is a group ring with QG`W etc.initialized.
 * Ideals is a list of ideals (the central conductors), wnb is a basis over the maximal order
 * and Ath is a free AssOrd-module. Tries to find a generator over AssOrd by Enumeration.
 *
 * The set of units in MaxOrd is much tooooo big! We enumerate only a subset!
 * This does not work, only used for experiments.
 *
 *******************************************************************************************/

intrinsic S4_Enumeration(QG :: AlgGrp, Ideals :: List, wnb :: ModTupFldElt, Ath :: Rec) -> List
{}
    local n, m, a, b,R, M, X, S, E, Q, L, W, rho, phi, F, BasisF, psi, K, U, pi,
          U1, U2, u1, u2, v, lambda, beta, T;

    U := [* *];
    for i:=1 to #Ideals do
        n := Degree(QG`W[i]);
        OF := MaximalOrder( CoefficientRing(QG`W[i]) );
        R, f := quo< OF | Ideals[i] >;
        B := BasisMatrix(Ideals[i]);
        T := CartesianProduct( [ [1..Abs(B[i,i])] : i in [1..Nrows(B)] ] );
        EltsOfR := [ f( OF ! [x : x in t] ) : t in T ];
        M := MatrixRing(R, n);
        X := CartesianPower(EltsOfR, n);
        X := [ Matrix(n, 1, [y: y in x]) : x in X ];
        print "i = ", i, "    #X = ", #X;
        S := Sym(#X);
        E := GlGenerators(QG`W[i]);
        Q := [ S ! [Index(X, M!g*x) : x in X] : g in E ];
        L := sub< S | Q >;
        W := WordGroup(L);
        rho :=InverseWordMap(L);
        W := Codomain(rho);
        phi := hom< W->L | Q >;
        F := FreeGroup(#Generators(W));
        BasisF := [F.i : i in [1..#Generators(F)]];
        psi := hom< W->F | BasisF >;
        V := [E[1]^0];
        print "#L = ", #L;
        counter := 1;
        for y in L do
            if counter mod 100 eq 0 then print "  ", counter; end if;
            v := ElementToSequence(psi(rho(y)));
            if #v gt 0 then
                w := [];
                for i in v do
                    if i ge 0 then
                        Append(~w, E[i]);
                    else
                        Append(~w, E[-i]^-1);
                    end if;
                end for;
                T := &*w;
                // T := &*[E[i] : i in v];
                // Append(~V, T);
            end if;
            counter := counter+1;
            /* if counter gt 200 then
                break;
            end if; 
            */
        end for;
        Append(~U, V);
    end for;
    return U;

    print "Start to enumerate: Size of set ----> ", &*[#V : V in U];
    counter := 1;
    EnumSet := CartesianProduct([ [1..#U[i]] : i in [1..#U] ]);
    for y in EnumSet do
        lambda := WedIsoInv(QG, [* U[i][y[i]] : i in [1..#U] *]);
        if counter mod 100 eq 0 then printf "%o  ", counter; end if;
        beta := GroupRingAction(Ath, lambda, wnb);
        if IsZGModuleElt(ElementToSequence(beta), Ath) then
             print "HURRA HURRA***************************************************";
             return beta;
        end if;
        counter := counter+1;
    end for;
    return beta;
end intrinsic;




/********************************  A4_Enumeration ***************************************
 *
 * QG is a group ring with QG`W etc.initialized.
 * Ideals is a list of ideals (the central conductors), wnb is a basis over the maximal order
 * and Ath is a free AssOrd-module. Tries to find a generator over AssOrd by Enumeration.
 * This is much faster than the funtion Enumeration.
 *
 *******************************************************************************************/

intrinsic A4_Enumeration(QG :: AlgGrp, Ideals :: List, wnb :: ModTupFldElt, Ath :: Rec) -> List
{}
    local n, m, a, b,R, M, X, S, E, Q, L, W, rho, phi, F, BasisF, psi, K, U, pi,
          U1, U2, u1, u2, v, lambda, beta, T;

    n := 3;
    a,b := IsPrincipal(Ideals[3]);
    m := Abs( Integers() ! b);
    R := quo< Integers() | m*Integers() >;
    M := MatrixRing(R, n);

    X := [ Matrix(n, 1, [r,s,t]) : r in R, s in R, t in R ];
    S := Sym(#X);
    E := ElementaryMatrices(n);
    Q := [ S ! [xtoi(M!g*x, m) : x in X] : g in E ];
    L := sub< S | Q >;
    // W := WordGroup(L);
    rho :=InverseWordMap(L);
    W := Codomain(rho);
    phi := hom< W->L | Q >;
    F := FreeGroup(#Generators(W));
    BasisF := [F.i : i in [1..#Generators(F)]];
    psi := hom< W->F | BasisF >;

    K := CoefficientRing(QG`W[1]);
    U, pi := UnitGroup(K);
    U1 := [Matrix([[pi(u)]]) : u in U];

    K := CoefficientRing(QG`W[2]);
    U, pi := UnitGroup(K);
    U2 := [Matrix([[pi(u)]]) : u in U];

    print "Start to enumerate: Size of set ----> ",  #U1*#U2*#L;
    counter := 1;
    for y in L do
        v := psi(rho(y));
	v := ElementToSequence(v);
	if #v gt 0 then
	    T := &*[E[i] : i in v];
	    for u1 in U1 do
	        for u2 in U2 do
		    if counter mod 100 eq 0 then printf "%o  ", counter; end if;
	            lambda := WedIsoInv(QG, [* u1, u2, T *]);
		    beta := GroupRingAction(Ath, lambda, wnb);
		    if IsZGModuleElt(ElementToSequence(beta), Ath) then
		        print "HURRA HURRA***************************************************";
			return beta;
		    end if;
		    counter := counter+1;
	        end for;
	    end for;
	end if;
    end for;
    return beta;
end intrinsic;

/********************************  A4_nb ***************************************
 *
 * Finds a generator over AssOrd for L := NumberField(pols[j]). Can be used with
 * the lists of A4-polynomials extracted from data base of Klueners/Malle in the file A4Mipos.m.
 *
 *******************************************************************************************/

intrinsic A4_nb(j :: RngIntElt, pols :: SeqEnum) -> List
{}
    local L, G, Aut, h, OL, theta, Ath, QG, AssOrd, ZG, rho, M, C, rep, W, MaxOrd, MaxAth,
          wnb, Zbasis, Z, F, CentralF, Ideals, InfinitePlaces, nb;

    L := NumberField(pols[j][1]);
    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);
    ZG := ModuleInit( IdentityMatrix(Rationals(), #G) );  /* the integral group ring */

    ModuleCompare(QG, ZG, AssOrd);
    IsOrder(QG, AssOrd);

    rho := RegularRep(QG);
    M := ZGModuleInit(Ath`hnf, rho);

    C := LocallyFreeClassgroup(G, true);
    print "cl(ZG) = ", C`cl;  /* We never use this information */

    QG := C`QG;
    rep, W := A4_rep(QG);
    AddAttribute(AlgGrp, "W");
    AddAttribute(AlgGrp, "Phi");
    AddAttribute(AlgGrp, "PhiInv");
    QG`W := W;
    QG`Phi := ComputePhiMatrix(QG, rep, W);
    QG`PhiInv := (QG`Phi)^-1;

    MaxOrd := MaximalOrder(QG, AssOrd);
    MaxAth := ZGModuleInit( ModuleProd(QG, MaxOrd, Ath)`hnf, rho );
    Ath := ZGModuleInit( Ath`hnf, rho );
    b, w := VeryNiceBasis(QG, MaxAth);
    TestBasisOverMaxOrd(QG, MaxOrd, MaxAth, w)  ;
    wnb := w[1];

    Zbasis := QGCentre(QG, QG`X, QG`H);
    Z := ModuleInit (Matrix( [ ElementToSequence(Zbasis[i]) : i in [1..#Zbasis] ] ) );
    F := ModuleConductor(QG, QG`MaxOrd, AssOrd);
    CentralF := ModuleIntersection(QG, F, Z);
    Ideals, InfinitePlaces := ComputeGi(QG, CentralF, QG`OC);
    Ideals;

    nb := A4_Enumeration(QG, Ideals, wnb, Ath);
    // nb := Enumeration(QG, Ideals, wnb, Ath); /* Works also, but is slower, since more general */
    TestBasisOverMaxOrd(QG, AssOrd, Ath, [nb]);
    nb*Ath`inv;


    return [* L, theta, Ath, AssOrd, nb, h, QG *];
end intrinsic;


/********************************  Dn_nb ***************************************
 *
 * Finds a generator over AssOrd for L := NumberField(pols[j]). Can be used with
 * the lists of Dn-polynomials extracted from data base of Klueners/Malle in the files 
 * D4Mipos.m, D5Mipos.m, D6Mipos.m or D7Mipos.m.
 *
 *******************************************************************************************/

intrinsic Dn_nb(j :: RngIntElt, pols :: SeqEnum) -> List
{}
    local L, G, Aut, h, OL, theta, Ath, QG, AssOrd, ZG, rho, M, C, rep, W, MaxOrd, MaxAth,
          wnb, Zbasis, Z, F, CentralF, Ideals, InfinitePlaces, nb;

    L := NumberField(pols[j][1]);
    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);
    ZG := ModuleInit( IdentityMatrix(Rationals(), #G) );  /* the integral group ring */

    ModuleCompare(QG, ZG, AssOrd);
    IsOrder(QG, AssOrd);

    rho := RegularRep(QG);
    M := ZGModuleInit(Ath`hnf, rho);

    C := LocallyFreeClassgroup(G, true); 
    print "cl(ZG) = ", C`cl;  /* We never use this information */
    QG := C`QG;
    rep, W := Dn_rep(QG);
    AddAttribute(AlgGrp, "W");
    AddAttribute(AlgGrp, "Phi");
    AddAttribute(AlgGrp, "PhiInv");
    QG`W := W;
    QG`Phi := ComputePhiMatrix(QG, rep, W);
    QG`PhiInv := (QG`Phi)^-1;
    
    MaxOrd := MaximalOrder(QG, AssOrd);
    while not IsVeryNiceOrder(QG, MaxOrd) do
        print "Search a very nice maximal order"; 
        MaxOrd := MaximalOrder(QG, AssOrd);
    end while;
    MaxAth := ZGModuleInit( ModuleProd(QG, MaxOrd, Ath)`hnf, rho );
    Ath := ZGModuleInit( Ath`hnf, rho );
    b, w := VeryNiceBasis(QG, MaxAth);
    TestBasisOverMaxOrd(QG, MaxOrd, MaxAth, w)  ;
    wnb := w[1];

    Zbasis := QGCentre(QG, QG`X, QG`H);
    Z := ModuleInit (Matrix( [ ElementToSequence(Zbasis[i]) : i in [1..#Zbasis] ] ) );
    F := ModuleConductor(QG, QG`MaxOrd, AssOrd);
    CentralF := ModuleIntersection(QG, F, Z);
    Ideals, InfinitePlaces := ComputeGi(QG, CentralF, QG`OC);
    Ideals;

    nb := Enumeration(QG, Ideals, wnb, Ath);
    TestBasisOverMaxOrd(QG, AssOrd, Ath, [nb]);
    nb*Ath`inv;


    return [* L, theta, Ath, AssOrd, nb, h, QG *];
end intrinsic;

intrinsic AssociatedOrder(f :: RngUPolElt) -> Rec
{}
    local L, G, Aut, h, OL, theta, Ath, QG, AssOrd;

    L := SplittingField(f);
    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);
    return AssOrd;
end intrinsic;


/********************************  S4_nb ***************************************
 *
 * f is a polynomial of degree 4 and group S4. 
 * This function uses the recursiv method described in Section 7 of the paper with H.Johnston.
 *
 *******************************************************************************************/

intrinsic S4_nb(f :: RngUPolElt) -> List
{}
    local L, G, Aut, h, OL, theta, Ath, QG, AssOrd, ZG, rho, M, C, rep, W, MaxOrd, MaxAth,
          wnb, Zbasis, Z, F, CentralF, Ideals, InfinitePlaces, nb;

    L := SplittingField(f);
    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);
    ZG := ModuleInit( IdentityMatrix(Rationals(), #G) );  /* the integral group ring */

    ModuleCompare(QG, ZG, AssOrd);
    IsOrder(QG, AssOrd);

    rho := RegularRep(QG);
    M := ZGModuleInit(Ath`hnf, rho);

    C := LocallyFreeClassgroup(G, true);
    print "cl(ZG) = ", C`cl;  /* We never use this information */
    QG := C`QG;
    rep, W := S4_rep(G);
    AddAttribute(AlgGrp, "W");
    AddAttribute(AlgGrp, "Phi");
    AddAttribute(AlgGrp, "PhiInv");
    QG`W := W;
    QG`Phi := ComputePhiMatrix(QG, rep, W);
    QG`PhiInv := (QG`Phi)^-1;

    MaxOrd := MaximalOrder(QG, AssOrd);
    while not IsVeryNiceOrder(QG, MaxOrd) do
        print "Search a very nice maximal order";
        MaxOrd := MaximalOrder(QG, AssOrd);
    end while;
    MaxAth := ZGModuleInit( ModuleProd(QG, MaxOrd, Ath)`hnf, rho );
    Ath := ZGModuleInit( Ath`hnf, rho );
    isfree2, w := IsLocallyFree(QG, AssOrd, Ath, 2);
    isfree3, w := IsLocallyFree(QG, AssOrd, Ath, 3);

    if isfree2 and isfree3 then
        print "OL is locally free over the associated order";
        b, w := VeryNiceBasis(QG, MaxAth);
        TestBasisOverMaxOrd(QG, MaxOrd, MaxAth, w)  ;
        wnb := w[1];

        Zbasis := QGCentre(QG, QG`X, QG`H);
        Z := ModuleInit (Matrix( [ ElementToSequence(Zbasis[i]) : i in [1..#Zbasis] ] ) );
        F := ModuleConductor(QG, MaxOrd, AssOrd);
        CentralF := ModuleIntersection(QG, F, Z);
        Ideals, InfinitePlaces := ComputeGi(QG, CentralF, QG`OC);
        Ideals;
    
        lambda := QG ! ElementToSequence(wnb);
        Alpha := WedIso(lambda);
        ZBasisAth := [ WedIso(QG!y) : y in RowSequence(Ath`hnf) ];
        ZBasisAlpha := ComputeZBasisAlpha(QG, Alpha);
        S, H, U := S4ComputeS(QG, Alpha, ZBasisAth);
        QMatAlg := MatrixAlgebra(Rationals(), 24);
        Sinv := (QMatAlg!H)^-1;
        Sinv;
        assert TestS(ZBasisAlpha, ZBasisAth, H, U);
        EnumData := PreComputeS4EnumData(QG, Ideals);
        NBElt := Matrix(Rationals(), [[]]);

        bool, nb := RecursiveS4Enumeration(EnumData, Sinv, NBElt, 1);

        /* Reconstruction of the normal basis element */ U := [* *];
        k := 1;
        for i:=1 to 5 do
            R := CoefficientRing(QG`W[i]);
            n := Degree(QG`W[i]);
            v := [R!(nb[1,j]) : j in [k..k+n^2-1]];
            Append(~U, QG`W[i] ! v);
            k := k+n^2;
        end for;
        Beta := [* U[i] * Alpha[i] : i in [1..5] *];
        lambda := WedIsoInv(QG, Beta);
    
        return [* L, theta, Ath, AssOrd, ElementToSequence(lambda), h, QG *];
    else
        print "OL is not locally free over the associated order";
        return [* *];
    end if;
end intrinsic;


/********************************  Abelian_nb ***************************************
 *
 * Finds a generator over AssOrd for L := NumberField(f). It is assumed that L/Q is abelian.
 *
 *******************************************************************************************/

intrinsic Abelian_nb(f :: RngUPolElt) -> List
{}
    local L, G, Aut, h, OL, theta, Ath, QG, AssOrd, ZG, rho, M, C, rep, W, MaxOrd, MaxAth,
          wnb, Zbasis, Z, F, CentralF, Ideals, InfinitePlaces, nb;

    L := SplittingField(f);
    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);
    ZG := ModuleInit( IdentityMatrix(Rationals(), #G) );  /* the integral group ring */

    ModuleCompare(QG, ZG, AssOrd);
    IsOrder(QG, AssOrd);

    rho := RegularRep(QG);
    M := ZGModuleInit(Ath`hnf, rho);

    C := LocallyFreeClassgroup(G, true);
    print "cl(ZG) = ", C`cl;  /* We never use this information */

    QG := C`QG;
    rep, W := Abelian_rep(G);
    AddAttribute(AlgGrp, "W");
    AddAttribute(AlgGrp, "Phi");
    AddAttribute(AlgGrp, "PhiInv");
    QG`W := W;
    QG`Phi := ComputePhiMatrix(QG, rep, W);
    QG`PhiInv := (QG`Phi)^-1;

    MaxOrd := MaximalOrder(QG, AssOrd);
    while not IsVeryNiceOrder(QG, MaxOrd) do
        print "Search a very nice maximal order";
        MaxOrd := MaximalOrder(QG, AssOrd);
    end while;
    MaxAth := ZGModuleInit( ModuleProd(QG, MaxOrd, Ath)`hnf, rho );
    Ath := ZGModuleInit( Ath`hnf, rho );
    b, w := VeryNiceBasis(QG, MaxAth);
    TestBasisOverMaxOrd(QG, MaxOrd, MaxAth, w)  ;
    wnb := w[1];

    Zbasis := QGCentre(QG, QG`X, QG`H);
    Z := ModuleInit (Matrix( [ ElementToSequence(Zbasis[i]) : i in [1..#Zbasis] ] ) );
    F := ModuleConductor(QG, QG`MaxOrd, AssOrd);
    CentralF := ModuleIntersection(QG, F, Z);
    Ideals, InfinitePlaces := ComputeGi(QG, CentralF, QG`OC);
    Ideals;

    nb := Enumeration(QG, Ideals, wnb, Ath);
    TestBasisOverMaxOrd(QG, AssOrd, Ath, [nb]);
    nb*Ath`inv;


    return [* L, theta, Ath, AssOrd, nb, h, QG *];
end intrinsic;




/********************************  IsVeryNiceOrder ***************************************
 *
 * Checks whether the Z-order A is contained in \bigoplus M_{n_i}(OF_i).
 * We assume that QG has QG`W etc. initialized.
 *
 *******************************************************************************************/

intrinsic IsVeryNiceOrder(QG :: AlgGrp, A :: Rec) -> BoolElt
{}
    local k, i, j,Ms, BasisA, b, lambda, M;

    BasisA := [QG ! b : b in RowSequence(A`hnf)];
    Ms := [WedIso(lambda) : lambda in BasisA ];
    for k:=1 to #Ms do
        for M in Ms[k] do
            for i:=1 to Nrows(M) do
	        for j:=1 to Ncols(M) do
                    if M[i,j] ne 0 then             /* necessary because of a bug in MAGMA */
	                if not IsIntegral(M[i, j]) then
		            return false;
		        end if;
                    end if;
	        end for;
	    end for;
        end for;
    end for;
    return true;
end intrinsic;


/********************************  QIdempotents ***************************************
 *
 * Computes the irreducible idempotents of QG.
 * We assume that QG has QG`W etc. initialized.
 *
 *******************************************************************************************/

intrinsic QIdempotents(QG :: AlgGrp) -> SeqEnum
{}
    local J, idem, i;

    idem := [];
    for i:= 1 to #QG`W do
        J := [* ScalarMatrix(M, 0) : M in QG`W *];
        J[i] := ScalarMatrix(QG`W[i], 1);
	Append(~idem, WedIsoInv(QG, J));
    end for;
    return idem;
end intrinsic;

/******************************** S4ComputeS ***************************************
 *
 * This function computes the matrix H (notation as in Section 7). Alpha is the generator
 * over the maximal order as an element of the Wedderburn decomposition, hence it is a list
 * of matrices over the rationals.  ZBasisAth is the Z-basis of Ath, again as elements in
 * the Wedderburn decomposition.
 *
 *******************************************************************************************/

intrinsic S4ComputeS(QG :: AlgGrp, Alpha :: List, ZBasisAth :: SeqEnum) -> AlgMatElt, AlgMatElt,AlgMatElt
{}
    local S, row, c, comp, n, B, Q, k, l, b, v, j, H, U;

    S := [];
    for row := 1 to #Group(QG) do
        c := [];
	for comp:=1 to #QG`X do
        // for comp:=#QG`X to 1 by -1  do
	    n := Ncols(Alpha[comp]);
	    B := ZBasisAth[row][comp];
	    Q := []; b := [];
	    for k:=1 to n do
	        for l:=1 to n do
		    Append(~b, Rationals()! B[k,l]);
	            v := [Rationals() ! 0 : i in [1..n^2]];
		    for j:=1 to n do
		        v[(k-1)*n+j] := Alpha[comp][j,l];
		    end for;
		    Append(~Q, v);
		end for;
	    end for;
	    Q := Matrix(Q);
	    b := Matrix(n^2,1,b);
	    c := c cat ElementToSequence(Q^-1 * b);
	end for;
	Append(~S, c);
    end for;
    S := MakeMatrixIntegral(  Matrix(S)  );
    H, U := HermiteForm(S);
    return S, H, U;
end intrinsic;

	    
/******************************** S4ComputeS ***************************************
 *
 * Just used for testing.
 *
 *******************************************************************************************/

intrinsic TestS(ZBasisAlpha :: SeqEnum, ZBasisAth :: SeqEnum, H :: AlgMatElt, U :: AlgMatElt) -> BoolElt
{}
    local n, r, row, e, f, col;
    
    n := #ZBasisAlpha;
    r := #ZBasisAlpha[1];
    for row := 1 to n do
        e := [* ZeroMatrix(Rationals(), Nrows(ZBasisAlpha[1][comp]), Nrows(ZBasisAlpha[1][comp])) : comp in [1..r] *];
	f := [* ZeroMatrix(Rationals(), Nrows(ZBasisAlpha[1][comp]), Nrows(ZBasisAlpha[1][comp])) : comp in [1..r] *];

        for col := 1 to n do
	    for comp := 1 to r do
	       e[comp] := e[comp] + H[row][col]*ZBasisAlpha[col][comp]; 
	       f[comp] := f[comp] + U[row][col]*ZBasisAth[col][comp];
	    end for;
	end for;
        for comp := 1 to r do
            if e[comp] ne f[comp] then
                return false;
            end if;
        end for; 
    end for;
    return true;
end intrinsic;


intrinsic ComputeZBasisAlpha(QG :: AlgGrp, Alpha :: List)-> SeqEnum
{}
    local ZBasisAlpha, comp, n, i, j, e, k;

    ZBasisAlpha := [ ];
    for comp := 1 to #QG`X do
    // for comp := #QG`X to 1 by -1 do
        n := Nrows(Alpha[comp]);
        for i := 1 to n do
            for j:= 1 to n do
                e := [* ZeroMatrix(Rationals(), Nrows(Alpha[i]), Nrows(Alpha[i])) : i in [1..#QG`X] *];
                for k:=1 to n do
                    e[comp][i,k] := Alpha[comp][j, k];
                end for;
                Append(~ZBasisAlpha, e);
            end for;
        end for;
    end for;
    return ZBasisAlpha;
end intrinsic;

/******************************** PreComputeS4EnumData ***************************************
 *
 * This function computes the sets U_i of Corollary 2.4 and some additional technical data.
 * This data is used in each step of the recursive enumeration.
 *
 *******************************************************************************************/

intrinsic PreComputeS4EnumData(QG :: AlgGrp, Ideals :: List) -> List
{}
    local n, m, a, b,R, M, X, S, E, Q, L, W, rho, phi, F, BasisF, psi, K, U, pi,
          U1, U2, u1, u2, v, lambda, beta, T;

    EnumData := [* *];
    for i:=1 to #Ideals do
        n := Degree(QG`W[i]);
	QMatAlg := MatrixAlgebra(Rationals(), n);
        OF := MaximalOrder( CoefficientRing(QG`W[i]) );
        R, f := quo< OF | Ideals[i] >;
        B := BasisMatrix(Ideals[i]);
        T := CartesianProduct( [ [1..Abs(B[i,i])] : i in [1..Nrows(B)] ] );
        EltsOfR := [ f( OF ! [x : x in t] ) : t in T ];
        M := MatrixRing(R, n);
        X := CartesianPower(EltsOfR, n);
        X := [ Matrix(n, 1, [y: y in x]) : x in X ];
        print "i = ", i, "    #X = ", #X;
        S := Sym(#X);
        E := [QMatAlg ! e : e in GlGenerators(QG`W[i])];
        Q := [ S ! [Index(X, M!g*x) : x in X] : g in E ];
        L := sub< S | Q >;
        W := WordGroup(L);
        rho :=InverseWordMap(L);
        W := Codomain(rho);
        phi := hom< W->L | Q >;
        F := FreeGroup(#Generators(W));
        BasisF := [F.i : i in [1..#Generators(F)]];
        psi := hom< W->F | BasisF >;
        V := [E[1]^0];
        Append(~EnumData, [* L, psi, rho, E *] );
        print "#L = ", #L;
    end for;
    return EnumData;
end intrinsic;


/******************************** RecursiveS4Enumeration ***************************************
 *
 * This function does the recursive enumeration described in Section 7.
 * EnumData is precomputed and contains data needed to compute the elemnts of the sets U_i.
 * T is the matrix H^-1 (notation of Section 7).
 * NBElt is the normal basis element constructed so far in the recursion.
 * depth is the depth of the recursion.
 *
 *******************************************************************************************/

intrinsic RecursiveS4Enumeration(EnumData :: List, T :: Mtrx, NBElt :: Mtrx, depth :: RngIntElt) -> BoolElt, Mtrx
{}
    local counter, L, psi, rho, E, p, rows, P, y, v, w, i, U, nb;

    L := EnumData[depth][1];
    psi := EnumData[depth][2];
    rho := EnumData[depth][3];
    E := EnumData[depth][4];

    print "depth = ", depth, "    NBElt = ", NBElt;
    if depth eq 1 then 
        p := 1; 
    else 
       p := &+[Nrows(EnumData[i][4][1])^2 : i in [1..depth-1]] + 1;
    end if;
    rows := Nrows( E[1]  )^2;
    P := Submatrix(T, 1, p, p+rows-1, rows);
    // print "P = ", P;

    counter := 0;
    for y in L do
        counter := counter + 1;
        // if counter mod 10000 eq 0 then printf "[%o, %o]  ", depth, counter; end if;
        if counter mod 1000000 eq 0 then print "depth = ", depth, " counter = ", counter; end if;
        v := ElementToSequence(psi(rho(y)));
        if #v gt 0 then
            w := [];
            for i in v do
                if i ge 0 then
                    Append(~w, E[i]);
                else
                    Append(~w, E[-i]^-1);
                end if;
            end for;
            U := &*w;
	else
	    U :=  E[1]^0;
        end if;
	nb := HorizontalJoin(NBElt, Matrix(Rationals(), [ ElementToSequence(U) ]));

        if IsIntegralMatrix(nb*P) then
        // if IsIntegral(nb, P) then
	    if depth lt 5 then
                printf "depth = %o, counter = %o\n", depth, counter;
	        NBEltFound, nb := RecursiveS4Enumeration(EnumData, T, nb, depth+1);
		if NBEltFound eq true then
		    return true, nb;
		end if;
            else
	       return true, nb;
	    end if;
	end if;
    end for;
    return false, NBElt;
end intrinsic;

intrinsic IsIntegral(v :: Mtrx, Q :: Mtrx) -> BoolElt
{}
    local a, col, i;

    for col := 1 to Ncols(Q)  do
        a := &+[v[1,i]*Q[i, col] : i in [1..Nrows(Q)] ];
        if not IsIntegral(a) then
	    return false;
	end if;
    end for;
    return true;
end intrinsic;
    

intrinsic IsIntegralMatrix(Q :: Mtrx) -> BoolElt
{}
    local i, j;

    for i:=1 to Nrows(Q) do
        for j:= 1 to Ncols(Q) do
	    if not IsIntegral(Q[i,j]) then
	        return false;
	    end if;
	end for;
    end for;
    return true;
end intrinsic;


intrinsic S4RandomUnits(QG :: AlgGrp, AssOrd :: Rec, EnumData :: List) -> List
{}
    local counter, L, psi, rho, E, p, q, U, Y, Units, y,v, T, lambda;

    L := [ EnumData[i][1] :i in [1..5]];
    psi := [ EnumData[i][2] :i in [1..5]];
    rho := [ EnumData[i][3] :i in [1..5]];
    E := [* EnumData[i][4] :i in [1..5]*];

    D,p,q := DirectProduct(L);


    Units := [];
    for counter := 1 to 10^6 do
        if counter mod 1000 eq 0 then print "  ", counter; end if;
        U := [* *];
        Y := [* *];
        for comp:=1 to 5 do
            y := Random(L[comp]);
            Append(~Y, y);
            v := ElementToSequence(psi[comp](rho[comp](y)));
            if #v gt 0 then
                w := [];
                for i in v do
                    if i ge 0 then
                        Append(~w, E[comp][i]);
                    else
                        Append(~w, E[comp][-i]^-1);
                    end if;
                end for;
                T := &*w;
                Append(~U, T);
            else
                Append(~U, E[comp][1]);
            end if;
        end for;
        lambda := WedIsoInv(QG, U);
        if IsModuleElt(QG, lambda, AssOrd) then
            print "Unit found !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
            Append(~Units,  &*[p[i](Y[i]) : i in [1..5]]);
        end if;
    end for;

    return Units;
end intrinsic;