declare verbose INB, 3;
SetVerbose("epsConj", 1);

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(1) : 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;

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

    G := Domain(h);
    fac := [y[1] : y in Factorization(p*OL)];
    P := fac[1];
    
    Fq, mFq := ResidueClassField(P);
    for b in Fq do
        D := Matrix([ ElementToSequence( b^(p^i) ) : i in [0..Degree(Fq)-1] ]);
	if Determinant(D) ne 0 then
	    alpha := (mFq^-1)(b);
	    break;
	end if;
    end for;
    m := [OL!0 : y in fac];
    m[1] := alpha;
    theta := ChineseRemainderTheorem(m, fac);
    return theta;
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;
    print "Compute the matrix A";
    A := Matrix( [ElementToSequence( OL ! h(g)(theta)  ) : g in G] );
    b := ZeroMatrix(Rationals(), 1, n);
    Atheta := [];
    print "Solve systems of linear equations";
    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;
    print "Init Atheta, then return";
    Atheta := ModuleInit(Matrix(Atheta));
    return Atheta;
end intrinsic;


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

    G := Domain(h);
    n := #G;
    S := MatrixAlgebra(Rationals(), n) ! BasisMatrix(I);
    A := Matrix( [ElementToSequence( OL ! h(g)(theta)  ) : g in G] );
    A := A * S^-1;
    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;


intrinsic ComputeAtheta(OL :: RngOrd, I :: RngOrdFracIdl, h :: Map, delta :: FldNumElt) -> Rec
{}
    local n, A, b, i, x, Atheta;;

    G := Domain(h);
    n := #G;
    S := MatrixAlgebra(Rationals(), n) ! BasisMatrix(I);

    D := Denominator(delta);
    A := 1/D * Matrix( [ElementToSequence( OL ! (D * h(g)(delta))  ) : g in G] );
    A := A * S^-1;
    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;

/******************************** Q8FindOFBasis ***************************************
 *
 * 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(Q), if ind = 1,...,4, 
 * QG`W[ind] = [ H_8, ...] if ind = 5.
 * V is a free OF-module.
 *
 * Computes an OF-basis of V, if ind = 1,..,4, or a Delta-basis of V, if ind = 5 and where Delta
 * is the maximal order in the quaternion algebra H_8.
 *
 * Attention: One should always use the more general intrinsic QuaternionFindOFBasis
 *            The name is misleading and should be changed!!!!!!!!!!!!!!!!!
 *
 *******************************************************************************************/

intrinsic Q8FindOFBasis(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;

    
    if ind ne 5 then
        F := CoefficientRing(QG`W[ind]);
    else
        F := CoefficientRing(QG`W[ind, 1]);
    end if;

    n := Integers() ! (Nrows(V`hnf) / Degree(F));
    Fbasis := Basis(F);
    Vbasis := [ Vector(Rationals(), c) : c in RowSequence(V`hnf) ];
    if ind ne 5 then
        BasisOverF := [ Vbasis[1] ];

        J := [* Zero(QG`W[i]) : i in [1..4] *];
        Append(~J, Zero(QG`W[5,1]) );
        FbasisInQG := [];
        for b in Fbasis do
                J[ind] := ScalarMatrix(QG`W[ind], b);
                lambda := Q8WedIsoInv(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(QG`W[i]) : i in [1..4] *];
        Append(~J, Zero(QG`W[5,1]));
        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 :=  Q8WedIsoInv(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 := Q8WedIsoInv(QG, J);
            Append(~BasisOverOF, GroupRingAction(V, 1/d * lambda, gen));
        end for;
        return BasisOverOF;
    else
       D := QG`W[5,1];
       Delta:=MaximalOrder(D);

       J := [* Zero(QG`W[i]) : i in [1..4] *];
       Append(~J, Zero(QG`W[5,1]) );
       DbasisInQG := [];
       for b in Basis(D) do
                J[ind] :=  b;
                lambda := Q8WedIsoInv(QG, J);
                Append(~DbasisInQG, lambda);
       end for;
       V_Qbasis := [ GroupRingAction(V, lambda, Vbasis[1]) : lambda in DbasisInQG ];
       A := Matrix(V_Qbasis);
       b := Solution(A, Vbasis);
       lambdas := [ D ! ElementToSequence(t) : t in b ];
       I := lideal<Delta | lambdas>;
       UnitIdeal := lideal<Delta | [One(Delta)]>;
       iota, gen := LeftIsomorphism(UnitIdeal, I);
       J[ind] :=  gen;
       lambda := Q8WedIsoInv(QG, J);
       gen := GroupRingAction(V, lambda, Vbasis[1]);
       return [ gen ];
    end if;

end intrinsic;


/******************************** QuaternionFindOFBasis ***************************************
 *
 * 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(Q), if it is a component
 * of Schurindex 1, otherwisw QG`W[ind] = [quaternion algebra D, some map, order Delta in D].
 * We do not allow matrix rings over D !!!!!!!!!!!!!!!!!!!!!!!!
 *
 * V is either a free OF-module or a free Delta-module. We do not check this freeness assumption, when
 * we compute OF-basis, but implicitly assume that OF has class number 1.
 * This should be changed.
 *
 * Computes an OF- or Delta-basis of V.
 *
 * Attention: We implicitly assume that OF has class number 1.
 *            The name is misleading and should be changed!!!!!!!!!!!!!!!!!
 *
 *******************************************************************************************/

intrinsic QuaternionFindOFBasis(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;


    if Type(QG`W[ind]) eq AlgMat then
        F := CoefficientRing(QG`W[ind]);
    else
        F := CoefficientRing(QG`W[ind, 1]);
    end if;

    n := Integers() ! (Nrows(V`hnf) / Degree(F));
    Fbasis := Basis(F);
    Vbasis := [ Vector(Rationals(), c) : c in RowSequence(V`hnf) ];
    if Type(QG`W[ind]) eq AlgMat then
        BasisOverF := [ Vbasis[1] ];

	J := [* *];
	for i:=1 to #QG`W do
	    if Type(QG`W[i]) eq AlgMat then
	        Append(~J, Zero(QG`W[i]) );
	    else
	        Append(~J, Zero(QG`W[i, 1]) );
	    end if;
	end for;
        FbasisInQG := [];
        for b in Fbasis do
                J[ind] := ScalarMatrix(QG`W[ind], b);
                lambda := QuaternionWedIsoInv(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 := [* *];
        for i:=1 to #QG`W do
            if Type(QG`W[i]) eq AlgMat then
                Append(~J, Zero(QG`W[i]) );
            else
                Append(~J, Zero(QG`W[i, 1]) );
            end if;
        end for;
        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 :=  QuaternionWedIsoInv(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 := QuaternionWedIsoInv(QG, J);
            Append(~BasisOverOF, GroupRingAction(V, 1/d * lambda, gen));
        end for;
        return BasisOverOF;
    else
       D := QG`W[ind,1];
       Delta:= QG`W[ind, 3];

       J := [* *];
       for i:=1 to #QG`W do
            if Type(QG`W[i]) eq AlgMat then
                Append(~J, Zero(QG`W[i]) );
            else
                Append(~J, Zero(QG`W[i, 1]) );
            end if;
       end for;
 
       Dbasis := [b*w : b in Basis(F), w in Basis(D)];
       DbasisInQG := [];
       for b in Dbasis do
                J[ind] :=  b;
                lambda := QuaternionWedIsoInv(QG, J);
                Append(~DbasisInQG, lambda);
       end for;
       V_Qbasis := [ GroupRingAction(V, lambda, Vbasis[1]) : lambda in DbasisInQG ];
       A := Matrix(V_Qbasis);
       b := Solution(A, Vbasis);
       lambdas := [ &+[ t[i]*Dbasis[i] : i in [1..#Dbasis] ]    : t in b ];
       I := lideal<Delta | lambdas>;
       UnitIdeal := lideal<Delta | [One(Delta)]>;
       assert IsLeftIsomorphic(UnitIdeal, I);
       if Degree(CoefficientRing(D)) eq 1 then     /* due to stupid MAGMA routines */
           _, gen := IsPrincipal(I);
       else
           _, gen := IsLeftIsomorphic(UnitIdeal, I);
       end if;
       J[ind] :=  gen;
       lambda := QuaternionWedIsoInv(QG, J);
       gen := GroupRingAction(V, lambda, Vbasis[1]);
       return [ gen ];
    end if;

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;


/******************************** Q8VeryNiceMoritaBasis ***************************************
 *
 * 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(Q), if ind = 1,..,4 and
 * QG`W[ind] = [H_8, some map, maximal order in H_8]
 * P is a M_n(OF)-module.
 *
 * Computes an M_n(OF)-basis of V, if one exists and ind = 1,..,4. Otherwise computes a Delta basis,
 * where Delta denotes the maximal order in H_8 containing 1,i,j,k.
 *
 * Attention: One should always use QuaternionVeryNiceMoritaBasis !!!!!!!!!!!!!!!!!!!!!!!!!!!
 *
 *******************************************************************************************/

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

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


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

    J := [* Zero(QG`W[i]) : i in [1..4] *];
    Append(~J, Zero( QG`W[5, 1] ));
    if ind ne 5 then
        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] := Q8WedIsoInv(QG, J);
            end for;
        end for;
    else
        E := ZeroMatrix(QG, 1, 1);
	J[ind] := One(QG`W[5,1]);
	E[1, 1] := Q8WedIsoInv(QG, J);
    end if;

    BasisP := [Vector(Rationals(),c) : c in RowSequence(P`hnf)];
    gens := [ GroupRingAction(P, E[1,1], b) : b in BasisP ];
    V := ZGSubModule(P, gens);
    OFBasisV := Q8FindOFBasis(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;



/******************************** QuaternionVeryNiceMoritaBasis ***************************************
 *
 * 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(Q), if it is a component
 * with Schurindex 1, otherwise QG`W[ind] = [quaternion field D, some map, maximal order Delta in D]
 * We do not allow matrix rings over quaternion fields.
 * P is a M_n(OF)-module.
 *
 * Computes an M_n(OF)-basis of V, if one exists,  respectively a Delta basis, if one exists.
 *
 * Attention: Implicitly we probably assume that that OF has trivial class number.
 *
 *******************************************************************************************/

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

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


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

    J := [* *];
    for i:=1 to #QG`W do
        if Type(QG`W[i]) eq AlgMat then
	    Append(~J, Zero(QG`W[i]) );
	else 
	    Append(~J, Zero(QG`W[i, 1]));
	end if;
    end for;
    if Type(QG`W[ind]) eq AlgMat then
        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] := QuaternionWedIsoInv(QG, J);
            end for;
        end for;
    else
        E := ZeroMatrix(QG, 1, 1);
        J[ind] := One(QG`W[ind,1]);
        E[1, 1] := QuaternionWedIsoInv(QG, J);
    end if;

    BasisP := [Vector(Rationals(),c) : c in RowSequence(P`hnf)];
    gens := [ GroupRingAction(P, E[1,1], b) : b in BasisP ];
    V := ZGSubModule(P, gens);
    OFBasisV := QuaternionFindOFBasis(QG, ind, V);
    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
        /* 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;

/******************************** Q8VeryNiceBasis ***************************************
 *
 * 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 = Z + Z + Z + Z + Delta, where Delta is the unique
 * maximal order in H_8 containing 1,i,j,k.
 *
 * Computes MaxOrd-basis of M, if one exists.
 * 
 * Attention:One should always us the more general routine QuaternionVeryNiceBasis
 *
 *******************************************************************************************/

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

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

    E := Q8Idempotents(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 5 do
        // print i, "-te Komponente";
        /* Berechne die Basis in der i-ten Komponente */
        isfree, v := Q8VeryNiceMoritaBasis(QG, i, C[i]);
        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;


/******************************** QuaternionVeryNiceBasis ***************************************
 *
 * 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, where the components are
 * either of the form M_{n_i}(OF_i) or Delta_i, where Delta_i is a maximal order in a totally
 * positive quaternion field. 
 *
 * 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 QuaternionVeryNiceBasis(QG :: AlgGrp, M :: Rec) -> BoolElt, List
{}

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

    E := QuaternionIdempotents(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 := QuaternionVeryNiceMoritaBasis(QG, i, C[i]);
        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. Works with any order, not just with maximal orders.
 *
 *******************************************************************************************/

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;

intrinsic pCubedRep(QG :: AlgGrp) -> List, SeqEnum
{}
    local Gl, G, chi, iota, p, G1, G2, A1, A2, a, b, c, Ma, Mb, Mc, Kchi, zeta, e, n, T, Omega, X, H, iso;

    G := Group(QG);
    p := Factorization(#G)[1,1];

    if Exponent(G) eq p^2 then
        A1<a,b> := Group< a,b | a^(p^2), b^p, b^(-1)*a*b = a^(1+p)>;
        G1, f := PermutationGroup(A1);
        a := f(a);
        b := f(b);

        e := Exponent(G1); n := #G1;
        T := CharacterTable(G1);
        Omega := CompOmega(G1);      /* Omega = Galois group of Q(zeta_e) / Q */
        X := OmegaRepresentatives(T, Omega);
        H := ComputeKi(G1, 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]];
        A := [* *];  B := [* *];
        for i:=1 to #chi do
            if Degree(chi[i]) eq 1 then
                Append( ~A, Gl[i] ! [iota[i](chi[i](a))] );
                Append( ~B, Gl[i] ! [iota[i](chi[i](b))] );
            end if;
           if Degree(chi[i]) eq p then
                Kchi := H[i][1];
                zeta := Kchi.1;
                Mb := DiagonalMatrix(Kchi, [zeta^(-j) : j in [0..p-1]]);
                Ma := ZeroMatrix(Kchi, p,p);
                for j:=2 to p do Ma[j,j-1] := Kchi!1; end for;
                Ma[1,p] := zeta;
                Append( ~A, Gl[i] ! Ma );
                Append( ~B, Gl[i] ! Mb );
            end if;
        end for;

        _, iso := IsIsomorphic(G, G1);
        return [* iso * hom<G1 -> Gl[i] | [a->A[i], b->B[i]]> : i in [1..#X] *], W;
     end if;
    if Exponent(G) eq p then
        A2<a,b,c> := Group< a,b,c | a^p, b^p, c^p, a*b = b*a*c, a*c = c*a, b*c = c*b>;
        G2, f := PermutationGroup(A2);
        a := f(a);
        b := f(b);
        c := f(c);

        e := Exponent(G2); n := #G2;
        T := CharacterTable(G2);
        Omega := CompOmega(G2);      /* Omega = Galois group of Q(zeta_e) / Q */
        X := OmegaRepresentatives(T, Omega);
        H := ComputeKi(G2, 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]];
        A := [* *];  B := [* *]; C:=[* *];
        for i:=1 to #chi do
            if Degree(chi[i]) eq 1 then
                Append( ~A, Gl[i] ! [iota[i](chi[i](a))] );
                Append( ~B, Gl[i] ! [iota[i](chi[i](b))] );
                Append( ~C, Gl[i] ! [iota[i](chi[i](c))] );
            end if;
           if Degree(chi[i]) eq p then
                Kchi := H[i][1];
                zeta := Kchi.1;
                Mb := DiagonalMatrix(Kchi, [zeta^(-j) : j in [0..p-1]]);
                Mc := DiagonalMatrix(Kchi, [zeta : j in [0..p-1]]);
                Ma := ZeroMatrix(Kchi, p,p);
                for j:=2 to p do Ma[j,j-1] := Kchi!1; end for;
                Ma[1,p] := Kchi!1;
                Append( ~A, Gl[i] ! Ma );
                Append( ~B, Gl[i] ! Mb );
                Append( ~C, Gl[i] ! Mc );
            end if;
        end for;
       _, iso := IsIsomorphic(G, G2);
        return [* iso * hom<G2 -> Gl[i] | [a->A[i], b->B[i], c->C[i]]> : i in [1..#X] *], W;
     end if;
end intrinsic;

intrinsic G63_rep(QG :: AlgGrp) -> List, SeqEnum
{}
    local Gl, G, chi, iota, p, G1, G2, A1, A2, a, b, c, Ma, Mb, Mc, Kchi, zeta, e, n, T, Omega, X, H, iso;

    G := Group(QG);

    G1 := [S : S in SmallGroups(63) | not IsAbelian(S) and Exponent(S) eq 21][1];
    a := G1.1; b := G1.2; c := G1.3;
    e := Exponent(G1); n := #G1;
    T := CharacterTable(G1);
    Omega := CompOmega(G1);      /* Omega = Galois group of Q(zeta_e) / Q */
    X := OmegaRepresentatives(T, Omega);
    H := ComputeKi(G1, 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]];
    A := [* *];  B := [* *]; C:=[* *];

    U := [A`subgroup : A in AbelianSubgroups(G1) | A`order eq 21][1];    /* U = <bc> */
    TU := CharacterTable(U);
    IndChars := [Induction(chi,G1) : chi in TU];

    for i:=1 to #chi do
        if Degree(chi[i]) eq 1 then
            Append( ~A, Gl[i] ! [iota[i](chi[i](a))] );
            Append( ~B, Gl[i] ! [iota[i](chi[i](b))] );
            Append( ~C, Gl[i] ! [iota[i](chi[i](c))] );
        end if;
        if Degree(chi[i]) eq 3 then
	    Kchi := H[i, 1];
	    d := Integers() ! Degree(chi[i]);
	    i0 := Index(IndChars, chi[i]);
	    rho := MonomialRepresentation(TU[i0], U, G1);
	    found := false;
	    while not found do
                vprintf INB, 1: "Search an integral monomial represenentation";
	        small_rho := WriteRepresentationOver(rho, Kchi);
		found := IsIntegralRepresentation(small_rho);
	    end while;
	    Ma := small_rho(a); Mb := small_rho(b); Mc := small_rho(c);
            Append( ~A, Gl[i]!Ma );
            Append( ~B, Gl[i]!Mb );
            Append( ~C, Gl[i]!Mc );
        end if;
    end for;
    _, iso := IsIsomorphic(G, G1);
    return [* iso * hom<G1 -> Gl[i] | [a->A[i], b->B[i], c->C[i]]> : i in [1..#X] *], W;
end intrinsic;

intrinsic IsIntegralRepresentation(rho :: Map) -> BoolElt
{}
    local i;

    for g in Generators(Domain(rho)) do
        if not &and[IsIntegral(z) : z in Eltseq(rho(g))] then
	    return false;
	end if;
    end for;
    return true;
end intrinsic;

intrinsic MonomialRepresentation(psi :: AlgChtrElt, U :: GrpPC, G1 :: GrpPC) -> Map
{}
    local Q, mQ, Qelts, R, d, Kpsi, Gl, W, Gens, images, g, sigma, D, i, j, rho;

    Q, mQ := quo< G1 | U>;
    Qelts := [q : q in Q];
    R := [(mQ^-1)(q) : q in Qelts];
    d := #R;
    Kpsi := CoefficientField(psi);
    Gl :=  GL(d,  Kpsi);
    W := MatrixRing(Kpsi, d);

    Gens := [g : g in Generators(G1)];
    images := [];
    for i:=1 to #Gens do
        g := Gens[i];
        sigma := Sym(d) ! [Index(Qelts, mQ( g*R[i] )) : i in [1..#Qelts]];
	D := ZeroMatrix(Kpsi, d, d);
	for j:=1 to d do
            D[j^sigma, j] := psi( R[j^sigma]^(-1)*g*R[j] );
	end for;
	Append(~images, D);
    end for;
    rho := hom<G1 -> Gl | [Gens[i]->images[i] : i in [1..#Gens]]>;
    return rho;
end intrinsic;



intrinsic GenDn_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);

    Us := NormalSubgroups(G : IsAbelian:=true, OrderEqual:=Integers()!(#G/2));
    for v in Us do
        V := v`subgroup;
        GenDn, tau, jot := GenDihedralGroup(V);
        b, iso := IsIsomorphic(G, GenDn); 
        if b then
            U := sub< GenDn | [jot(g) : g in Generators(V)]>;
            break; 
        end if;
    end for;
    e := Exponent(GenDn); n := #GenDn;
    T := CharacterTable(GenDn);
    TU := CharacterTable(U);
    Omega := CompOmega(GenDn);      /* Omega = Galois group of Q(zeta_e) / Q */
    X := OmegaRepresentatives(T, Omega);
    H := ComputeKi(GenDn, 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]];
    GenDnChars := [X[i][1] : i in [1..#X]];
    iota := [H[i][3] : i in [1..#X]];

    sigmas := [g : g in Generators(U) ];
    Sigmas := [ [* *] : s in sigmas ];  Tau := [* *];
    for i:=1 to #GenDnChars do
        if Degree(GenDnChars[i]) eq 1 then
            for j:=1 to #sigmas do
                Append( ~Sigmas[j], Gl[i] ! [iota[i](GenDnChars[i](sigmas[j]))] );
            end for;
            Append( ~Tau, Gl[i] ! [iota[i](GenDnChars[i](tau))] );
        end if;
        if Degree(GenDnChars[i]) eq 2 then
            print "i = ", i;
	    psi := GenDnChars[i];
	    chi := [eta : eta in TU | InnerProduct( eta, Restriction(psi, U) ) eq 1][1];
	    print "Order(chi) = ", Order(chi);
	    f := CyclotomicPolynomial(Order(chi));
	    Qpsi := CharacterField(psi);
	    Qchi := CharacterField(chi);
	    g := [g[1] : g in Factorization( PolynomialRing(Qpsi) ! f) | Degree(g[1]) eq 2] [1];
	    F<z> := NumberField(g);
	    _, emb := IsSubfield(Qchi, F);
            for j:=1 to #sigmas do
                print "j = ", j;
	        v1 := ElementToSequence(emb(chi(sigmas[j])) * 1); 
		v2 := ElementToSequence(emb(chi(sigmas[j])) * z);
                M := Matrix([ [v1[1], v2[1]],
                              [v1[2], v2[2]]
                            ]);
                Append( ~Sigmas[j], Gl[i] ! M );
            end for;
	    z1 := [z1[1] : z1 in Roots(PolynomialRing(F) ! g) | z1[1] ne z] [1];
	    v1 := ElementToSequence(F!1);
	    v2 := ElementToSequence(z1);
            N := Matrix([ [v1[1], v2[1]],
	                  [v1[2], v2[2]]
                       ]);
            Append( ~Tau, Gl[i] ! N );
        end if;
    end for;

    b, iso := IsIsomorphic(G, GenDn);
    return [* iso * hom<GenDn -> Gl[i] | 
                    [sigmas[j]->Sigmas[j,i] : j in [1..#sigmas]] cat [tau->Tau[i]]> : 
		                i in [1..#X] *], W;
end intrinsic;



/******************************** Q8_rep ***************************************
 *
 * Computes the representations of G, which is assumed to be isomorphic to Q8, the
 * quaternion  group of order 8.
 *
 *******************************************************************************************/

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

    G := Group(QG);

    Q8 := QuaternionGroup(2);
    e := Exponent(Q8); n := #Q8;
    T := CharacterTable(Q8);
    Omega := CompOmega(Q8);      /* Omega = Galois group of Q(zeta_e) / Q */
    X := OmegaRepresentatives(T, Omega);
    H := ComputeKi(Q8, 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] *];
    D<x,y,z> := QuaternionAlgebra<Rationals() | -1, -1>;
    Z := QuaternionOrder(Basis(D));
    UnitsZ, h := UnitGroup(Z);
    ux := [u : u in UnitsZ | h(u) eq x] [1];
    uy := [u : u in UnitsZ | h(u) eq y] [1];
    W[5] := [* D, h *];
    Gl[5] := UnitsZ;
    chi := [X[i][1] : i in [1..#X]];
    iota := [H[i][3] : i in [1..#X]];
    sigma := Q8.1;
    tau := Q8.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
            Append( ~Sigma, ux );
            Append( ~Tau, uy );
        end if;
    end for;

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

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

    G := Group(QG);

    C2 := CyclicGroup(2);
    Q8 := QuaternionGroup(2);
    C2Q8 := DirectProduct(C2, Q8);
    e := Exponent(Q8); n := #Q8;
    T := CharacterTable(C2Q8);
    Omega := CompOmega(C2Q8);      /* Omega = Galois group of Q(zeta_e) / Q */
    X := OmegaRepresentatives(T, Omega);
    H := ComputeKi(C2Q8, 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] *];
    D<x,y,z> := QuaternionAlgebra<Rationals() | -1, -1>;
    Z := QuaternionOrder(Basis(D));
    UnitsZ, h := UnitGroup(Z);
    minus1 := [u : u in UnitsZ | h(u) eq -1] [1];
    ux := [u : u in UnitsZ | h(u) eq x] [1];
    uy := [u : u in UnitsZ | h(u) eq y] [1];
    W[9] := [* D, h *];	W[10] := [* D, h *];
    Gl[9] := UnitsZ;    Gl[10] := UnitsZ;
    chi := [X[i][1] : i in [1..#X]];
    iota := [H[i][3] : i in [1..#X]];
    eta := C2Q8.1;
    sigma := C2Q8.2;
    tau := C2Q8.3;
    Eta := [* *]; Sigma := [* *];  Tau := [* *];
    for i:=1 to #chi do
        if Degree(chi[i]) eq 1 then
            Append( ~Eta, Gl[i] ! [iota[i](chi[i](eta))] );
            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
	    if chi[i](eta) eq 2 then   /* eta is mapped to 1 in D, iff chi(eta) = 2 */
                Append( ~Eta, Id(UnitsZ)  );
	    else                       /* eta is mapped to -1 in D, iff chi(eta) = -2 */
	        Append( ~Eta, minus1  );
	    end if;
            Append( ~Sigma, ux );
            Append( ~Tau, uy );
        end if;
    end for;

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



/******************************** Quaternion_rep ***************************************
 *
 * Computes the representations of G, which is assumed to be isomorphic to a generalized
 * quaternion group. This intrinsic uses the description of the Wederburn decomposition in
 * Curtis/Reiner, Part I. 
 *
 *******************************************************************************************/

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

    G := Group(QG);
    m := Integers() ! (#G / 4);

    Q4m := QuaternionGroup(m);

    if IsEven(m) then
        r := Valuation(m, 2);
	m0 := Integers() ! ( m / 2^r );
	ds := [ 2^(r+1)*d : d in Divisors(m0) ];
    else
        ds := [2*d : d in Divisors(m) | d gt 1 ];
    end if;

    
    e := Exponent(Q4m); n := #Q4m;
    T := CharacterTable(Q4m);
    Omega := CompOmega(Q4m);      /* Omega = Galois group of Q(zeta_e) / Q */
    X := OmegaRepresentatives(T, Omega);
    H := ComputeKi(Q4m, 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 := Q4m.1;
    tau := Q4m.2;
    Sigma := [* *];  Tau := [* *];
    Gl := [* *];
    W := [* *];
    for i:=1 to #chi do
        if SchurIndex(chi[i]) eq 1 then   /* abelian or dihedral */
	    Append(~Gl, GL( Integers()!Degree(X[i][1]), H[i][1] ) );
	    Append(~W, MatrixRing( H[i][1], Integers()!Degree(X[i][1]) ) );
            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 /* This is a 2-dim. dihedral component */
                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 if;
       if SchurIndex(chi[i]) eq 2 then    /* a quaternionic component */
           Kchi := H[i,1];
           KchiX<q> := PolynomialRing(Kchi);
           d := Selectd(ds,Kchi);
           g := [g[1] : g in Factorization(KchiX ! CyclotomicPolynomial(d)) | Degree(g[1]) eq 2][1];
           disc := Discriminant(g);
           /* Construct a quaternion order which contains a d-th root of 1 */
           if Degree(Kchi) eq 1 then
	       D<x,y,z> := QuaternionAlgebra<Rationals() | Rationals() ! disc, -1>;
	       Delta := MaximalOrder(D);
               zeta := (-Rationals()!Coefficient(g,1) + x) / 2;
	   else
               D<x,y,z> := QuaternionAlgebra<Kchi | disc, -1>;
               zeta := (-Coefficient(g,1) + x) / 2;
               Delta := MaximalOrder(Order(zeta));
	   end if;
	   U,h := NormOneGroup(Delta:ModScalars:=false);
	   Append(~Gl, U );
	   Append(~W, [* D, h *] );
	   V := [u :  u in U | Order(u) eq d];
	   ux := [v : v in V | y*h(v) eq Conjugate(h(v))*y][1];
	   uy := [u :  u in U | h(u) eq y][1];
	   Append( ~Sigma, ux );
	   Append( ~Tau, uy );
       end if;
    end for;

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



intrinsic Selectd(ds :: SeqEnum, Kchi) -> RngIntElt
{}
    local d, C, Cplus;


    for d in ds do
        if EulerPhi(d) eq Degree(Kchi)*2 then
	    if Degree(Kchi) eq 1 then return d; end if;
	    C<z> := CyclotomicField(d);
            Cplus := sub<C | z+z^-1>;
	    if IsIsomorphic(Kchi, Cplus)  then return d; end if;
	end if;
    end for;
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;


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

    Z := [];

    for g in Group(QG) do
        c := [];
        for i := 1 to 4 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;
        d := W[5, 2]( rep[5](g) );
        c cat:= ElementToSequence(d);
        Append(~Z, c);
        // Z cat:= c;
    end for;
    // Phi := Transpose( Matrix(Z) );
    Phi := Matrix(Z) ;

    return Phi;
end intrinsic;

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

    Z := [];

    for g in Group(QG) do
        c := [];
        for i := 1 to #rep do
           if Type(W[i]) eq AlgMat then		/* Schurindex = 1 */
               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;
	   else
	       d := ElementToSequence( W[i, 2]( rep[i](g) ) );
	       d := [ElementToSequence(x) : x in d];
	       for v in d do
	           c cat:= ElementToSequence(v);
	       end for;
	   end if;
        end for;
        Append(~Z, c);
    end for;
    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;

intrinsic QuaternionWedIso(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
        if Type(QG`W[i]) eq AlgMat then     /* Schurindex equal 1 */
            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);
        else 				   /* Case of a pos. def. quaternion algebra */	
	    K := CoefficientRing(QG`W[i, 1]);
	    length := Degree(K);
	    d := [];
	    for i:=1 to 4 do
	        Append(~d, K ! [c[first+t] : t in [0..length-1]]);
		first := first + length;
	    end for;
	    D := QG`W[i,1] ! d;
	    Append(~T, D);
        end if;
    end for;
    return T;
end intrinsic;

	

intrinsic Q8WedIso(lambda :: AlgGrpElt) -> List
{}
    local QG, c, d, 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 := [* *];
    for i:=1 to 4 do
        D := ZeroMatrix(Rationals(), 1, 1);
        D[1,1] := c[i];
        Append(~T, D);
    end for;
    d := [c[i] : i in [5..8]];
    Append(~T, QG`W[5,1] ! d);
    
    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;


intrinsic Q8WedIsoInv(QG :: AlgGrp, T :: List) -> AlgGrpElt
{}
    local c, i, k, l, lambda;

    c := [];
    for i:=1 to 4 do
       c cat:= ElementToSequence(T[i][1, 1]);
    end for;
    c cat:= ElementToSequence(T[5]);
    lambda := QG ! ElementToSequence( Vector(c)*QG`PhiInv );
    return lambda;
end intrinsic;


intrinsic QuaternionWedIsoInv(QG :: AlgGrp, T :: List) -> AlgGrpElt
{}
    local c, v, i, k, l, lambda;

    c := [];
    for i:=1 to #QG`W do
       if Type(QG`W[i]) eq AlgMat then 
           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;
       else
           v := ElementToSequence(T[i]);
           for k:=1 to 4 do
	       c cat:= ElementToSequence(v[k]);
	   end for;
       end if;
    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, BoolElt
{}
    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, found;

    U := [* *];
    for i:=1 to #Ideals do
        n := Degree(QG`W[i]);
	// OF := MaximalOrder( CoefficientRing(QG`W[i]) );
        OF := Order( Ideals[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]);
        print "Compute Q";
        Q := [ S ! [Index(X, M!g*x) : x in X] : g in E ];
        print "Q computed";
        L := sub< S | Q >;
        print "Compute WordGroup W";
        W := WordGroup(L);
        print "W computed";
        rho :=InverseWordMap(L);
        print "rho computed";
        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, true;
        end if;
        counter := counter+1;
    end for;
    return beta, false;
end intrinsic;



/********************************  QuaternionEnumeration ***************************************
 *
 * QG is a group ring with QG`W etc.initialized.
 * We assume that G is a generalized quaternion group. It is only reasonably tested for Q_8 and Q_12.
 *
 * 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 QuaternionEnumeration(QG :: AlgGrp, Ideals :: List, wnb :: ModTupFldElt, Ath :: Rec) -> List, BoolElt
{}
    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
        if Type(QG`W[i]) eq AlgMat then
            V := CompEnumerationSet(QG`W[i], Ideals[i]);
            Append(~U, V);
        else
	    V := CompEnumerationSet(QG`W[i], Ideals[i]);
            Append(~U, V);
        end if;
    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 := QuaternionWedIsoInv(QG, [* U[i][y[i]] : i in [1..#U] *]);
        if counter mod 10000 eq 0 then printf "%o  ", counter; end if;
        beta := GroupRingAction(Ath, lambda, wnb);
        if IsZGModuleElt(ElementToSequence(beta), Ath) then
             print "HURRA HURRA***************************************************";
             return beta, true;
        end if;
        counter := counter+1;
    end for;
    return beta, false;
end intrinsic;


/**************************** CompEnumerationSet *******************************************
 *
 *  Computes the enumeration set if the Wedderburn component is of the form M_n(OF).
 *
*/

intrinsic CompEnumerationSet(W :: AlgMat, I :: RngOrdIdl) -> SeqEnum
{}
    local n, OF, R, f, B, T, EltsOfR, M, X, S, E, Q, L, WG, phi, F, 
          BasisF, psi, V, counter, v, y, i, x, w;

    n := Degree(W);
    OF := MaximalOrder( CoefficientRing(W) );
    R, f := quo< OF | I >;
    B := BasisMatrix(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 ];
    S := Sym(#X);
    E := GlGenerators(W);
    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;
           Append(~V, T);
       end if;
       counter := counter+1;
    end for;
    return V;
end intrinsic;


/**************************** CompEnumerationSet *******************************************
 *
 *  Computes the enumeration set if the Wedderburn component is a positively definite quaternion algebra.
 *  This is only reasonably tested for Q_8 and Q_12. For all other generalized quaternion groups the
 *  coefficient fields are extensions of Q. This is implemented, but because of lack of examples not tested. 
 *  It runs for the only Q_16 extension we were able to construct. This extension is not tame, so we
 *  cannot predict whether OLis free over AssOrd or not. Need several tamely ramified Q_16 examples.
 *
*/

intrinsic CompEnumerationSet(W :: List, I :: RngOrdIdl) -> SeqEnum
{}
    local n, OF, R, f, B, T, EltsOfR, M, X, S, E, Q, L, WG, phi, F,
          U_Delta, h1, h2, U_OF, BasisF, psi, V, counter, v, y, i, x, w, alpha;

    D := W[1];
    Delta := W[3];
    U_Delta, h := UnitGroup(Delta);
    // E := [Delta ! h(u) : u in U_Delta]; /* elements of Delta are wrt. a pseudo basis */
    E := [D ! h(u) : u in U_Delta]; /* elements of Delta are wrt. a pseudo basis */

    OF := Order( I );
    U_OF, h := UnitGroup(OF);

    if Ngens(U_OF) eq 1 then         /* only finitely many units */
        return E;
    else                             /* This should technically work, but is not tested !!!!!! */
        n := 4;
        R, f := quo< OF | I >;
        B := BasisMatrix(I);
        T := CartesianProduct( [ [1..Abs(B[i,i])] : i in [1..Nrows(B)] ] );
        EltsOfR := [ f( OF ! [x : x in t] ) : t in T ];
        Q := [];
	BasisDelta, alpha := NormalizeBasis(Delta);           /* Work with an OF-basis, not just a pseudo basis */
        for a in BasisDelta do
            for b in BasisDelta do
	        v := ElementToSequence( Delta ! (a*b) );
                Q cat:= [OF ! (v[i] / alpha[i]) : i in [1..#v]];
            end for;
        end for;
        Q := [f(x) : x in Q];
        DeltaBar := Algebra<R, 4 | Q>;

        X := CartesianPower(EltsOfR, n);
        X := [ DeltaBar ! [y: y in x] : x in X ];
        S := Sym(#X);
        E cat:= [Delta ! [h(u), 0, 0, 0] : u in Generators(U_OF)]; /* We assume that the first basis element is 1 */
        
        /* The next step has to be improved because this takes much to long. For example, for 
           a Q_20 extension one has #E = 62, #X = 160000.
        */ 
        Q := [ S ! [Index(X, ( DeltaBar ! [f(c) : c in NormalizedCoeff(OF, Delta!g, alpha)] ) * 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;
               Append(~V, T);
           end if;
           counter := counter+1;
        end for;
        return V;
    end if;
end intrinsic;


/**************************************************************************************
 *
 *  Delta is given by a OF-pseudo basis. We always assume that OF is a principal ideal domain.
 *  We return a OF-basis BasisDelta and the list of generators of the ideals in the pseudo basis.
*/

intrinsic NormalizeBasis(Delta :: AlgAssVOrd) -> SeqEnum, SeqEnum
{}
    local BasisDelta, alpha, s, gen;

    BasisDelta := [];
    alpha := [];
    for s in PseudoBasis(Delta) do
        assert IsPrincipal(s[1]);
        _, gen := IsPrincipal(s[1]);
        Append(~BasisDelta, gen * s[2]);
	Append(~alpha, gen);
    end for;
    return BasisDelta, alpha;
end intrinsic;

/**************************************************************************************
 *
 *  If g is an element in Delta, then it is given with respect to the pseudo basis. This intrinsic
 *  returns the coefficients of g with respect to the OF-basis previously computed with
 *  NormalizeBasis.
*/

intrinsic NormalizedCoeff(OF :: RngOrd, g :: AlgAssVOrdElt, alpha :: SeqEnum) -> SeqEnum
{}
    local v, i;

    v := ElementToSequence(g);
    return [OF ! (v[i] / alpha[i]) : i in [1..#v] ];

end intrinsic;


/************************************************************************************************
 *
 * No longer need. Use QuaternionEnumeration instead.
 *
*/
intrinsic Q8Enumeration(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 4 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;

    D := QG`W[5,1];
    Delta := MaximalOrder(D);
    U_Delta, h := UnitGroup(Delta);
    V := [D ! h(u) : u in U_Delta];
    Append(~U, V);


    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 := Q8WedIsoInv(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 true, beta;
        end if;
        counter := counter+1;
    end for;
    return false, 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,BoolElt
{}
    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, true;
		    end if;
		    counter := counter+1;
	        end for;
	    end for;
	end if;
    end for;
    return beta, false;
end intrinsic;


intrinsic INB(K :: FldNum : WeakINB:=false) -> List, BoolElt
{}
    local implemented, found, X, A4, S4, C2Q8;

    /* 
    G, Aut, h := AutomorphismGroup(K);
    h := map<Domain(h)->Codomain(h) | g:->h(g^-1)>;
    */

    A := ArtinRepresentations(K);
    L :=K`artinrepdata`K;
    G := K`artinrepdata`G;
    act := K`artinrepdata`act;
    h := map<G->Parent(act(G.1)) | g:->act(g^-1)>;
  
    implemented := false;
    found := false;

    A4 := AlternatingGroup(4);
    S4 := SymmetricGroup(4);
    C2Q8 := DirectProduct( CyclicGroup(2), QuaternionGroup(2) );

    if IsAbelian(G) then
        X, found := Abelian_nb(L, G, h : WeakINB);
	implemented := true;
    elif IsDihedral(G) then
        X, found := Dn_nb(L, G, h : WeakINB);
	implemented := true;
    elif IsGenDihedral(G) then
        X, found := GenDn_nb(L, G, h : WeakINB);
	implemented := true;
    elif IsQuaternion(G) then
        X, found := Quaternion_nb(L, G, h : WeakINB);
        implemented := true;
    elif IsIsomorphic(G, A4) then
        X, found := A4_nb(L, G, h : WeakINB);
        implemented := true;
    elif IsIsomorphic(G, S4) then
        X, found := S4_nb(L, G, h : WeakINB);
        implemented := true;
    elif IsIsomorphic(G, C2Q8 : WeakINB) then
        X, found := C2Q8_nb(L, G, h);
        implemented := true;
    end if;

    if not implemented then
        print "Algorithm only implemented for abelian, dihedrald or quaternion groups, A_4 and S_4. Sorry!";
        return [**], found;
    else
        return X, found;
    end if;
end intrinsic;
 
      
intrinsic IsQuaternion(G :: Grp) -> BoolElt
{}
    local D;

    if not IsDivisibleBy(#G, 4) then
        return false;
    end if;

    D := QuaternionGroup( Integers() ! (#G/4) );
    return IsIsomorphic(G, D);
end intrinsic;




/********************************  A4_nb ***************************************
 *
 * Finds a generator over AssOrd for an A4 extension.
 *
 *******************************************************************************************/

intrinsic A4_nb(L :: FldNum, G :: GrpPerm, h :: Map : WeakINB:=false) -> List, BoolElt
{}
    local OL, theta, Ath, QG, AssOrd, ZG, rho, M, C, rep, W, MaxOrd, MaxAth,
          wnb, Zbasis, Z, F, CentralF, Ideals, InfinitePlaces, nb, found;

    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];
    if WeakINB then
        return [* L, theta, Ath, AssOrd, wnb, h, QG *], true;
    end if;

    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);

    nb, found := A4_Enumeration(QG, Ideals, wnb, Ath);
    // nb, found := 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 *], found;
end intrinsic;


intrinsic pCubed_nb(L :: FldNum, G :: GrpPerm, AL :: RngOrdFracIdl, h :: Map : WeakINB:=false) -> List, BoolElt
{}
    local OL, theta, Ath, QG, AssOrd, ZG, rho, M, C, rep, W, MaxOrd, MaxAth,
          wnb, Zbasis, Z, F, CentralF, Ideals, InfinitePlaces, nb, found;

    OL := MaximalOrder(L);

    theta := NormalBasisElement(OL, h);  
    Ath := ComputeAtheta(OL, AL, 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 := pCubedRep(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 := VeryNiceMaximalOrder(QG);
    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];
    if WeakINB then
        return [* L, theta, Ath, AssOrd, wnb, h, QG *], true;
    end if;


    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);

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

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


intrinsic G63_nb(L :: FldNum, G :: GrpPerm, h :: Map : WeakINB:=false) -> List, BoolElt
{}
    local OL, theta, Ath, QG, AssOrd, ZG, rho, M, C, rep, W, MaxOrd, MaxAth,
          wnb, Zbasis, Z, F, CentralF, Ideals, InfinitePlaces, nb, found, p;

    OL := MaximalOrder(L);

    p := 7;
    P := Factorisation(p*OL)[1,1];
    // AL := p^-1*P; /* AL is the square root of the inverse different. */
    // AL := P^(-13);

    DL := Different(OL);
    AL := Sqrt(DL)^-1;
    // AL := p^2*AL;

    vprint INB, 1:  "Compute a normal basis element";
    theta := NormalBasisElement(OL, h, 5); 
    // theta := NormalBasisElement(OL, h);
    print "Compute Atheta";
    Ath := ComputeAtheta(OL, AL, 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, rho);

    // C := LocallyFreeClassgroup(G, true);
    // print "cl(ZG) = ", C`cl;  /* We never use this information */
    // QG := C`QG;
    print "Compute representations";
    rep, W := G63_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 := VeryNiceMaximalOrder(QG);
    vprint INB, 1: "Compute MaxAt";
    MaxAth := ZGModuleInit( ModuleProd(QG, MaxOrd, Ath), rho );
    Ath := ZGModuleInit( Ath, rho );
    vprint INB, 1:  "Compute a basis over the maximal order";
    b, w := VeryNiceBasis(QG, MaxAth);
    // TestBasisOverMaxOrd(QG, MaxOrd, MaxAth, w)  ;
    wnb := w[1];
    if WeakINB then
        return [* L, theta, Ath, AssOrd, wnb, h, QG, MaxOrd, MaxAth *], true;
    end if;


    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);

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

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




intrinsic VeryNiceMaximalOrder(QG :: AlgGrp) -> Rec
{}
    local MaxOrdBasis, ind, F, OF, OFBasis, J, i, j, b, lambda, M, MaxOrd;
    
    MaxOrdBasis := [];
    for ind:=1 to #QG`W do
        F := CoefficientRing(QG`W[ind]);
        OF := MaximalOrder( F );
        OFbasis := [F!b : b in Basis(OF)];
        J := [* Zero(w) : w in QG`W *];
        for i:=1 to Degree(QG`W[ind]) do
            for j:=1 to Degree(QG`W[ind]) do
                for b in OFbasis do
                    J[ind][i,j] := b;
                    lambda := WedIsoInv(QG, J);
                    if not &and[IsIntegral(z) : z in Coefficients(63*lambda)] then 
                        // print "ind, i, j = ", ind, i, j;
                    end if;
                    Append(~MaxOrdBasis, lambda);
                    J[ind][i,j] := Zero(F);
                end for;
            end for;
        end for;
    end for;
    M :=[Coefficients(lambda) : lambda in MaxOrdBasis];
    MaxOrd := ModuleInit(Matrix(M));
    return MaxOrd;
end intrinsic;





/********************************  Dn_nb ***************************************
 *
 * Finds a generator over AssOrd for a dihedral extension.
 *
 *******************************************************************************************/

intrinsic Dn_nb(L :: FldNum, G :: GrpPerm, h :: Map : WeakINB:=false) -> List, BoolElt
{}
    local OL, theta, Ath, QG, AssOrd, ZG, rho, M, C, rep, W, MaxOrd, MaxAth,
          wnb, Zbasis, Z, F, CentralF, Ideals, InfinitePlaces, nb, found;

    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];
    if WeakINB then
        return [* L, theta, Ath, AssOrd, wnb, h, QG *], true;
    end if;


    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);

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


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

intrinsic GenDn_nb(L :: FldNum, G :: GrpPerm, h :: Map : WeakINB:=false) -> List, BoolElt
{}
    local OL, theta, Ath, QG, AssOrd, ZG, rho, M, C, rep, W, MaxOrd, MaxAth,
          wnb, Zbasis, Z, F, CentralF, Ideals, InfinitePlaces, nb, found;

    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 := GenDn_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];
    if WeakINB then
        return [* L, theta, Ath, AssOrd, wnb, h, QG *], true;
    end if;

    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);

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


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





/********************************  Quaternion_nb ***************************************
 *
 * Finds a generator over AssOrd for a quaternion extension.
 *
 *******************************************************************************************/

intrinsic Quaternion_nb(L :: FldNum, G :: GrpPerm, h :: Map : WeakINB:=false) -> List, BoolElt
{}
    local found, OL, theta, Ath, QG, AssOrd, ZG, rho, M, C, rep, W, MaxOrd, MaxAth,
          wnb, Zbasis, Z, F, CentralF, Ideals, InfinitePlaces, nb;

    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 */
    rho := RegularRep(QG);
    P := ZGModuleInit(Ath`hnf, rho);
    

    if IsTamelyRamified(L) then
        isfree := Zero(C`cl) eq ClassGroupLog(C, P);
        print "---------------------------------------->>>>>>>>  OL is free? ", isfree;
    end if;
    
    rep, W := Quaternion_rep(QG);
    QG`W := W;
    QG`Phi := QuaternionComputePhiMatrix(QG, rep, W);
    QG`PhiInv := (QG`Phi)^-1;

    MaxOrd := MaximalOrder(QG, AssOrd);
    while not QuaternionIsVeryNiceOrder(QG, MaxOrd) do
        print "Search a very nice maximal order";
        MaxOrd := MaximalOrder(QG, AssOrd);
    end while;
    for i:=1 to #QG`W do
        if Type(QG`W[i]) ne AlgMat then
	    if Degree( CoefficientRing(QG`W[i,1]) ) eq 1 then
	        Delta := MaximalOrder(QG`W[i,1]);
	    else
	        Delta:=Order([QuaternionWedIso( QG!v )[i] :v in RowSequence(MaxOrd`hnf)]);
	        assert IsMaximal(Delta);
	    end if;
	    Append(~QG`W[i],Delta);
	end if;
    end for;
   
    MaxAth := ZGModuleInit( ModuleProd(QG, MaxOrd, Ath)`hnf, rho );
    Ath := ZGModuleInit( Ath`hnf, rho );
    b, w := QuaternionVeryNiceBasis(QG, MaxAth);

    TestBasisOverMaxOrd(QG, MaxOrd, MaxAth, w)  ;
    wnb := w[1];
    if WeakINB then
        return [* L, theta, Ath, AssOrd, wnb, h, QG *], true;
    end if;


    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);

    nb, found := QuaternionEnumeration(QG, Ideals, wnb, Ath);
    TestBasisOverMaxOrd(QG, AssOrd, Ath, [nb]);
    nb*Ath`inv;
    if IsTamelyRamified(L) then assert found eq isfree; end if;


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

/********************************  C2Q8_nb ***************************************
 *
 * Finds a generator over AssOrd for a C2 x Q8 - extension.
 *
 *******************************************************************************************/

intrinsic C2Q8_nb(L :: FldNum, G :: GrpPerm, h :: Map : WeakINB:=false) -> List, BoolElt
{}
    local found, OL, theta, Ath, QG, AssOrd, ZG, rho, M, C, rep, W, MaxOrd, MaxAth,
          wnb, Zbasis, Z, F, CentralF, Ideals, InfinitePlaces, nb;

    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 */
    rho := RegularRep(QG);
    P := ZGModuleInit(Ath`hnf, rho);


    if IsTamelyRamified(L) then
        isfree := Zero(C`cl) eq ClassGroupLog(C, P);
        print "---------------------------------------->>>>>>>>  OL is free? ", isfree;
    end if;

    rep, W := C2Q8_rep(QG);
    QG`W := W;
    QG`Phi := QuaternionComputePhiMatrix(QG, rep, W);
    QG`PhiInv := (QG`Phi)^-1;

    MaxOrd := MaximalOrder(QG, AssOrd);
    while not QuaternionIsVeryNiceOrder(QG, MaxOrd) do
        print "Search a very nice maximal order";
        MaxOrd := MaximalOrder(QG, AssOrd);
    end while;
    for i:=1 to #QG`W do
        if Type(QG`W[i]) ne AlgMat then
            if Degree( CoefficientRing(QG`W[i,1]) ) eq 1 then
                Delta := MaximalOrder(QG`W[i,1]);
            else
                Delta:=Order([QuaternionWedIso( QG!v )[i] :v in RowSequence(MaxOrd`hnf)]);
                assert IsMaximal(Delta);
            end if;
            Append(~QG`W[i],Delta);
        end if;
    end for;

    MaxAth := ZGModuleInit( ModuleProd(QG, MaxOrd, Ath)`hnf, rho );
    Ath := ZGModuleInit( Ath`hnf, rho );
    b, w := QuaternionVeryNiceBasis(QG, MaxAth);

    TestBasisOverMaxOrd(QG, MaxOrd, MaxAth, w)  ;
    wnb := w[1];
    if WeakINB then
        return [* L, theta, Ath, AssOrd, wnb, h, QG *], true;
    end if;


    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);

    nb, found := QuaternionEnumeration(QG, Ideals, wnb, Ath);
    TestBasisOverMaxOrd(QG, AssOrd, Ath, [nb]);
    nb*Ath`inv;
    // if IsTamelyRamified(L) then assert found eq isfree; end if;

    return [* L, theta, Ath, AssOrd, nb, h, QG *], found;
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(L :: FldNum, G :: GrpPerm, h :: Map : WeakINB:=false) -> List, BoolElt
{}
    local OL, theta, Ath, QG, AssOrd, ZG, rho, M, C, rep, W, MaxOrd, MaxAth,
          wnb, Zbasis, Z, F, CentralF, Ideals, InfinitePlaces, nb, found;

    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];
        if WeakINB then
            return [* L, theta, Ath, AssOrd, wnb, h, QG *], true;
        end if;


        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);
    
        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);
        if bool eq false then
            return [**], false;
        end if;

        /* 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 *], bool;
    else
        print "OL is not locally free over the associated order";
        return [* *], false;
    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(L :: FldNum, G :: GrpPerm, h :: Map : WeakINB:=false) -> List, BoolElt
{}
    local OL, theta, Ath, QG, AssOrd, ZG, rho, M, C, rep, W, MaxOrd, MaxAth,
          wnb, Zbasis, Z, F, CentralF, Ideals, InfinitePlaces, nb, found;

    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];
    if WeakINB then
        return [* L, theta, Ath, AssOrd, wnb, h, QG *], true;
    end if;


    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);

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


    return [* L, theta, Ath, AssOrd, nb, h, QG *], found;
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;


/********************************  QuaternionIsVeryNiceOrder ***************************************
 *
 * We assume that G is a generalized quaternion group.
 * Checks whether the components of the the Z-order A which are of Schurindex 1 are 
 * contained in M_{n_i}(OF_i).
 * No check for the quaternionic components is performed. (We don't need this because in quaternionic
 * components we construct the maximal order Delta from MaxOrd.
 *
 * We assume that QG has QG`W etc. initialized.
 *
 *******************************************************************************************/

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

    BasisA := [QG ! b : b in RowSequence(A`hnf)];
    Ms := [QuaternionWedIso(lambda) : lambda in BasisA ];
    for k:=1 to #Ms do
       for M in Ms[k] do
           if Type(M) eq AlgMatElt then
                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 if;
       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;


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

    idem := [];
    W := QG`W;
    for i:= 1 to 5 do
        J := [* ScalarMatrix(W[i], 0) : i in [1..4] *];
        Append(~J, Zero(W[5,1]));
        if i ne 5 then
             J[i] := ScalarMatrix(QG`W[i], 1);
        else
             J[i] := One(W[5,1]);
        end if;
        Append(~idem, Q8WedIsoInv(QG, J));
    end for;
    return idem;
end intrinsic;


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

    idem := [];
    W := QG`W;
    for i:= 1 to #W do
        J := [* *];
	for j:=1 to #W do
	    if Type(W[j]) eq AlgMat then  
	       if j eq i then 
	           Append(~J, ScalarMatrix(QG`W[j], 1));
	       else
	           Append(~J, ScalarMatrix(QG`W[j], 0));
               end if;
	    else
	       if j eq i then
	           Append(~J, One(W[j,1]));
	       else
	           Append(~J, Zero(W[j,1]));
	       end if;
	    end if;
	end for;
        Append(~idem, QuaternionWedIsoInv(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;

