/*********************  Changes since 25.9.2017 **************************

26.9.2017 In Quotient replaced Extension by LiftCharacter since Extension does no longer exist.
          Fixed a bug: in order to identify the inflation chi of a character of the quotient only the fixed
          representative of the orbit was considered. Now we identify the orbit X[i] such that chi is in X[i].

	  In addition, we then set psi := X[i][2] and compute iota: Kchi -> Kchi such that iota(psi) = chi. 
          Using iota in the following guarantees compatabilty of the representations of DT for G and its quiotient Q = G/U.

27.9.2017 Inserted a second form of the intrinsic Approximation with a different list of parameters; needed for
          computations of local Artin symbols.

28.9.2017 Replace print by vprintf - this has just begun and needs to be completed.

**************************** end of Changes **********************************/



// verbose comments
declare verbose RelKTh, 5;
SetVerbose("RelKTh", 0);



declare attributes AlgGrp : MaxOrd, X, Mods, Reps, H, Z, OC, B, Psi, Brauer,  W, Phi, PhiInv;



/* The Rec ZModule is the basic structure for modules over the integers.
   We view Z-modules as element in Q^n, so each element can be represented
   by a row vector of length n with coefficients in Q. The matrix hnf is
   a matrix with rational coefficients, each row representing a module element.
   We always assume hnf to be in row Hermite Normal Form. denom is the denominator of
   hnf, and inv the inverse of hnf.
*/

ZModule := recformat<hnf : Mtrx,
                     denom : RngIntElt,
                     inv : Mtrx
                    >;

/* The Rec ZGModule is the basic structure for modules over ZG.
   We view Z-modules as subsets of  Q^n = Q v_1 \oplus ... \oplus Q v_n =: V,
   so each element can be represented
   by a row vector of length n with coefficients in Q. The matrix hnf is
   a matrix with rational coefficients, each row representing a module element.
   We always assume hnf to be in row Hermite Normal Form. denom is the denominator of
   hnf, and inv the inverse of hnf.
   phi : G -> Gl_n(Q) is a homomorphism which gives the G-module structure. The j-th
   column in phi(g) represents the g(v_j).
   If (a_1, ..., a_n) is a row vector representing an element a in V, then g(a) is
   represented by the row vector (a_1, ..., a_n)*Transpose(phi(g)).
*/
ZGModule := recformat<hnf : Mtrx,
                      denom : RngIntElt,
                      inv : Mtrx,
                      phi : Map
                    >;


/* FG-modules are given by a M over a finite field F and a homomorphism phi:G->Gl_n(F) */

FGModule := recformat<M : ModTupFld,
                     phi : HomGrp
                    >;


/* BrauerStructure consists of 3 components: chi is a character of G, Uphis is a list
of triples [U, phi, c_{U,phi}] which gives the Brauer induction of chi, iota is a hom from
the field of values of chi (as given by the character table) to K_chi.
*/

BrauerStructure := recformat<chi : AlgChtrElt,
                             Uphis : SeqEnum,
                             iota : Map
                            >;


/******************************* LocFreeClassGrp *****************************************
 *
 * This is the basic structure which discribes the locally free class group.
 *
********************************************************************************/

LocFreeClassGrp := recformat<cl : GrpAb,
                             f : Map,
                             D : GrpAb,
                             rcgp : GrpAb,
                             m : List,
                             InfinitePlaces : SeqEnum,
                             Ideals : List,
                             QG : AlgGrp
                            >;


/******************************* DTGrp *****************************************
 *
 * This is the basic structure which discribes the relative torsion DT(Z_p[G]).
 *
********************************************************************************/
DTGrp := recformat<DT : GrpAb,       /* DT(Z_p[G]) as abelian group                         */
                   OCmodG : GrpAb,   /* (OC / \frg_p)^\times as abelian group               */
                   m : List,         /* a list of maps as given by MultiplicativeGroup      */
                   Gp : List,        /* the list of ideals (\frg_{1,p}, ..., \frg_{r,p})    */
                   f : Map,          /* the map from (OC / \frg_p)^\times to DT             */
                   p : RngIntElt,
                   QG : AlgGrp
                  >;


/******************************* K0RelGrp *****************************************
 *
 * This is the basic structure which discribes the relative group K_0(Z_p[G], Q_p).
 * Most record fields are as in DTGrp.
 *
********************************************************************************/

K0RelGrp := recformat<DT : GrpAb,
                      PrimeIdeals : List,   /* a List of SeqEnum of prime ideals above p */
                      PrimeElts : List,     /* corresponding List of SeqEnum of uniformizing elements */
                      OCmodG : GrpAb,
                      m : List,
                      Gp : List,
                      f : Map,
                      p : RngIntElt,
                      QG : AlgGrp
                    >;


/******************************* K0RelElt *****************************************
 *
 * This is the basic structure which discribes an element of the relative group K_0(Z_p[G], Q_p).
 * A and B are rec ZGModule and Theta is a matrix in Gl_m(QG) which describes an isomorphism
 * A \tensor Q --> B \tensor Q.
 *
********************************************************************************/

K0RelElt := recformat<A : Rec,
                      B : Rec,
                      Theta : Mtrx
                    >;



/* ****************************** RelativeGroup ******************************
 *
 * RelativeGroup(QG::AlgGrp, p :: RngIntElt) -> Rec: Computes the group K_0(Z_p[G], Q_p)
 * and enough additional data to solve the discrete logarithm problem.
 *
 * It is assumed that QG has already assigned all the values computed for example in
 * LocallyFreeClassgroup. This version of RelativeGroup is for example used in ClassGroupLog.
 *
  ************************************************************************* */


intrinsic RelativeGroup(QG :: AlgGrp, p :: RngIntElt) -> Rec
{}

    local e, n,  E, Z, ZG, CentralF, F, Ideals, InfinitePlaces;

    e := Exponent(Group(QG)); n := #Group(QG);
    E := QIdempotents(QG, QG`X, QG`H);     /* Q-rational idempotents */

    Z := ModuleInit (Matrix( [ ElementToSequence(QG`OC[i]) : i in [1..#QG`OC] ] ) );

    ZG := ModuleInit( IdentityMatrix(Rationals(), n) );  /* the integral group ring */

    CentralF := CentralConductor(QG);
    F := ModuleProd(QG, CentralF, ZG);

    // print "Compute the ideals g_i";
    Ideals, InfinitePlaces := ComputeGi(QG, CentralF, QG`OC);

    // print "Compute ResidueClassGroup";
    // OCmodG, m := ResidueClassGroup(QG, Ideals);

    if not assigned(QG`Brauer) then
        print "Compute Brauer induction";
        QG`Brauer := InitBrauer(QG`H);
    end if;

    // print "Compute the relative group";
    K :=  K0Rel(QG, ZG, F, Ideals, p);
    K`QG := QG;

   return K;
end intrinsic;

intrinsic RelativeGroup(QG :: AlgGrp, p :: RngIntElt, neglect :: SetEnum) -> Rec
{}

    local e, n,  E, Z, ZG, CentralF, F, Ideals, InfinitePlaces;

    e := Exponent(Group(QG)); n := #Group(QG);
    E := QIdempotents(QG, QG`X, QG`H);     /* Q-rational idempotents */

    Z := ModuleInit (Matrix( [ ElementToSequence(QG`OC[i]) : i in [1..#QG`OC] ] ) );

    ZG := ModuleInit( IdentityMatrix(Rationals(), n) );  /* the integral group ring */

    CentralF := CentralConductor(QG);
    F := ModuleProd(QG, CentralF, ZG);

    // print "Compute the ideals g_i";
    Ideals, InfinitePlaces := ComputeGi(QG, CentralF, QG`OC);

    for i in neglect do
        Ideals[i] := 1 * QG`H[i][4];
    end for;
    

    // print "Compute ResidueClassGroup";
    // OCmodG, m := ResidueClassGroup(QG, Ideals);

    if not assigned(QG`Brauer) then
        // print "Compute Brauer induction";
        QG`Brauer := InitBrauer(QG`H);
    end if;

    // print "Compute the relative group";
    K :=  K0Rel(QG, ZG, F, Ideals, p);
    K`QG := QG;

   return K;
end intrinsic;


/* ****************************** RelativeGroup ******************************
 *
 * RelativeGroup(G::Grp, p :: RngIntElt) -> Rec: Computes the group K_0(Z_p[G], Q_p)
 * and enough additional data to solve the discrete logarithm problem.
 *
  ************************************************************************* */

intrinsic RelativeGroup(G :: Grp, p :: RngIntElt) -> Rec
{}

    local e, n, QG, T, Omega, E, Zbasis, Z, ZG, CentralF, F, Ideals,
          InfinitePlaces; 

    e := Exponent(G); n := #G;
    QG := GroupAlgebra(Rationals(), G); 

    vprintf RelKTh, 1:  "Compute Character table\n";
    T := CharacterTable(G);

    vprintf RelKTh, 1:  "Compute Omega orbits\n";
    Omega := CompOmega(G);      /* Omega = Galois group of Q(zeta_e) / Q */
    QG`X := OmegaRepresentatives(T, Omega);

    vprintf RelKTh, 1:  "Compute character fields\n";
    QG`H := ComputeKi(G, QG`X);

    QG`B := ComputeB(QG, QG`H);

    vprintf RelKTh, 1:  "Compute data for Psii\n";
    QG`Psi := ComputePsi(QG);

    vprintf RelKTh, 1:  "Compute Q-rational idempotents\n";
    E := QIdempotents(QG, QG`X, QG`H);     /* Q-rational idempotents */
    
    vprintf RelKTh, 1:  "Compute the maximal order of the centre of QG\n";
    Zbasis := QGCentre(QG, QG`X, QG`H);  /* NEVER CHANGE THIS LIST: it is used in the computation of PsiInv */
    QG`OC := Zbasis;
    Z := ModuleInit (Matrix( [ ElementToSequence(Zbasis[i]) : i in [1..#Zbasis] ] ) );

    ZG := ModuleInit( IdentityMatrix(Rationals(), n) );  /* the integral group ring */

    CentralF := CentralConductor(QG);
    F := ModuleProd(QG, CentralF, ZG);

    // print "Compute the ideals g_i";
    Ideals, InfinitePlaces := ComputeGi(QG, CentralF, Zbasis);

    vprintf RelKTh, 1:  "Compute Brauer induction\n";
    QG`Brauer := InitBrauer(QG`H);

    // print "Compute ResidueClassGroup";
    // OCmodG, m := ResidueClassGroup(QG, Ideals);

    vprintf RelKTh, 1:  "Compute the relative group\n";
    K :=  K0Rel(QG, ZG, F, Ideals, p);
    K`QG := QG;

   return K;
end intrinsic;

			
/* ****************************** LocallyFreeClassgroup ******************************
 * LocallyFreeClassgroup(G::Grp, WithMaxOrd::BoolElt) -> GrpAb, GrpAb: Computes the locally free classgroup of
 * the order ZG. As well as the kernel group D(ZG).
 * If WithMaxOrd eq true then we compute a maximal order and compute F = conductor of ZG in the maximal order.
 * Otherwise we use F = g * ZG, where g denotes the central conductor. g is computed using
 * [Curtis/Reiner I, (27.13)].
 *
  ************************************************************************* */


intrinsic LocallyFreeClassgroup(G :: Grp, WithMaxOrd :: BoolElt) -> Rec
{}
    local e, n, QG, T, Omega, X, H, E, Zbasis, Z, ZG, MaxOrd, F, CentralF, Ideals,
          InfinitePlaces, rcgp, cl, f, UseK1, AmodF, exact, ps, qs, K1AmodF, C;

    e := Exponent(G); n := #G;
    QG := GroupAlgebra(Rationals(), G); 

    vprintf RelKTh, 1: "Compute Character table\n";
    T := CharacterTable(G);

    vprintf RelKTh, 1: "Compute Omega orbits\n";
    Omega := CompOmega(G);      /* Omega = Galois group of Q(zeta_e) / Q */
    X := OmegaRepresentatives(T, Omega);
    QG`X := X;

    vprintf RelKTh, 1: "Compute character fields\n";
    H := ComputeKi(G, X);
    QG`H := H;

    QG`B := ComputeB(QG, H);

    vprintf RelKTh, 1: "Compute data for Psii\n";
    QG`Psi := ComputePsi(QG);

    vprintf RelKTh, 1: "Compute Q-rational idempotents\n";
    E := QIdempotents(QG, X, H);     /* Q-rational idempotents */
    
    vprintf RelKTh, 1: "Compute the maximal order of the centre of QG\n";
    Zbasis := QGCentre(QG, X, H);  /* NEVER CHANGE THIS LIST: it is used in the computation of PsiInv */
    QG`OC := Zbasis;
    Z := ModuleInit (Matrix( [ ElementToSequence(Zbasis[i]) : i in [1..#Zbasis] ] ) );

    ZG := ModuleInit( IdentityMatrix(Rationals(), n) );  /* the integral group ring */

    if WithMaxOrd then
        vprintf RelKTh, 1: "Compute a maximal order in QG containing ZG\n";
        MaxOrd := MaximalOrder(QG, ZG);
        QG`MaxOrd := MaxOrd;

        vprintf RelKTh, 1: "Compute the conductor F and the central conductor g\n";
        F := ModuleConductor(QG, MaxOrd, ZG);                /* the conductor */       
        CentralF := ModuleIntersection(QG, F, Z);            /* the central conductor */  
    else
        CentralF := CentralConductor(QG);
	F := ModuleProd(QG, CentralF, ZG);
    end if;

    vprintf RelKTh, 1: "Compute the ideals g_i\n";
    Ideals, InfinitePlaces := ComputeGi(QG, CentralF, Zbasis);

    vprintf RelKTh, 1: "Compute the ray class groups\n";
    IgmodPg, K, m := RayClassGroup(QG, Ideals, InfinitePlaces);

    vprintf RelKTh, 1: "Compute K_1(ZG / F)\n";
    UseK1 := true;
    K1AmodF, exact, ps, qs := Oumf(QG, ZG, F, Ideals, UseK1);


    /* Ps, Qs will not be needed. In earlier versions it was used for testing. */
    Ps := [ModuleSum(QG, F, ModuleProd(QG, p, ZG)) : p in ps];
    Qs := [ModuleSum(QG, F, ModuleProd(QG, q, ZG)) : q in qs];

    vprintf RelKTh, 1: "Compute the relation matrix \n";
    clA, D, f := RelationMatrix(QG, IgmodPg, K, m, K1AmodF);

/*
    QG`Brauer := InitBrauer(QG`H);
    OCmodG, m := ResidueClassGroup(QG, Ideals);
    K0Rels := [ K0Rel(QG, ZG, F, Ideals, p) : p in PrimeDivisors(#G) ];
*/
    C := rec< LocFreeClassGrp | >;
    C`cl := clA;
    C`f := f;
    C`D := D;
    C`rcgp := IgmodPg;
    C`m := m;
    C`InfinitePlaces := InfinitePlaces;
    C`Ideals := Ideals;
    C`QG := QG;

    return C;

end intrinsic;


/* ****************************** LocallyFreeClassgroup ******************************
 *
 * LocallyFreeClassgroup(QG::AlgGrp, WithMaxOrd::BoolElt) -> GrpAb, GrpAb: 
 * Computes the locally free classgroup of
 * the order ZG. As well as the kernel group D(ZG).
 * If WithMaxOrd eq true then we compute a maximal order and compute F = conductor of ZG in the maximal order.
 * Otherwise we use F = g * ZG, where g denotes the central conductor. g is computed using
 * [Curtis/Reiner I, (27.13)].
 *
 * It is assumed that QG has already assigned all the values computed in LocallyFreeClassGroup(G, WithMaxOrd)
 * or in RelativeGroup(G, p).
 *
  ************************************************************************* */

intrinsic LocallyFreeClassgroup(QG :: AlgGrp, WithMaxOrd :: BoolElt) -> Rec
{}
    local n, H, E, Zbasis, Z, ZG, F, CentralF, Ideals,
          InfinitePlaces, rcgp, cl, f, UseK1, AmodF, exact, Ps, Qs, K1AmodF, C;

    n := #Group(QG);
    // E := QIdempotents(QG, X, H);     /* Q-rational idempotents */

    Z := ModuleInit (Matrix( [ ElementToSequence(QG`OC[i]) : i in [1..#QG`OC] ] ) );

    ZG := ModuleInit( IdentityMatrix(Rationals(), n) );  /* the integral group ring */

    if WithMaxOrd then
        print "Compute the conductor F and the central conductor g";
        F := ModuleConductor(QG, QG`MaxOrd, ZG);                /* the conductor */
        CentralF := ModuleIntersection(QG, F, Z);            /* the central conductor */
    else
        CentralF := CentralConductor(QG);
        F := ModuleProd(QG, CentralF, ZG);
    end if;

    print "Compute the ideals g_i";
    Ideals, InfinitePlaces := ComputeGi(QG, CentralF, QG`OC);

    print "Compute the ray class groups";
    IgmodPg, K, m := RayClassGroup(QG, Ideals, InfinitePlaces);

    print "Compute K_1(ZG / F)";
    UseK1 := true;
    K1AmodF, exact, ps, qs := Oumf(QG, ZG, F, Ideals, UseK1);


    print "Compute the relation matrix ";
    clA, D, f := RelationMatrix(QG, IgmodPg, K, m, K1AmodF);

    C := rec< LocFreeClassGrp | >;
    C`cl := clA;
    C`f := f;
    C`D := D;
    C`rcgp := IgmodPg;
    C`m := m;
    C`InfinitePlaces := InfinitePlaces;
    C`Ideals := Ideals;
    C`QG := QG;

    return C;

end intrinsic;


/* ****************************** CompOmega ******************************
 * CompOmega(G::Grp) -> SeqEnum: Just computes (Z / e)^* as a list of integers, where
 * e denotes the exponent of G.
 *
  ************************************************************************* */

intrinsic CompOmega(G::Grp) -> SeqEnum
{ Computes the Galois group of Q(zeta_e)/Q as a list of integers, where
    e is the exponent of G}

    local e, Omega, Q, h;

    e := Exponent(G);
    Q := quo<Integers() | e>;
    Q,h := MultiplicativeGroup(Q);
    Omega := [Integers()!h(g) : g in Q];
    return Omega;
end intrinsic;


/* ****************************** OmegaRepresentatives ******************************
 *
 * OmegaRepresentatives(T::SeqEnum, Omega::SeqEnum) -> SeqEnum:
 * Computes a list X of Omega orbits of absolute irreducible characters.
 * Each X[i] is a list of characters. X[i][1] should be considered as our choice of
 * fixed representative of the orbit X[i].
 *
************************************************************************************* */

intrinsic OmegaRepresentatives(T::SeqEnum, Omega::SeqEnum) -> SeqEnum
{Computes the Omega orbits of the absolute irreducible characters}

    local X, flag, i, chi, j, k, psi;
    
    X := [];
    flag := [false : i in T];
    i := 1;
    while i le #T do
	chi := T[i];
	Xi := [];
	for j in Omega do
	    psi := GaloisConjugate(chi, j);
	    for k in [i..#T] do
		if not flag[k] and psi eq T[k] then
		   flag[k] := true;
		   Xi := Append(Xi, T[k]);
		end if;
	    end for;
	end for;
	while i le #T and flag[i] do
	    i := i+1;
	end while;
	X := Append(X, Xi);
    end while;
    return X;
end intrinsic;


intrinsic OmegaOrbitsOfReps(IrrMod::SeqEnum, Omega::SeqEnum) -> List, List, SeqEnum
{Computes the Omega orbits of the absolute irreducible characters}

    local X, flag, i, chi, j, k, psi;

    X := [ ];
    Reps := [* *];
    Mods := [* *];
    flag := [false : i in IrrMod];
    i := 1;
    IrrRep := [ Representation(V) : V in IrrMod ];
    while i le #IrrRep do
        chi := Character(IrrRep[i]);
        Xi := [ ];
        Ri := [* *];
        Vi := [* *];
        for j in Omega do
            psi := GaloisConjugate(chi, j);
            for k in [i..#IrrRep] do
                if not flag[k] and psi eq Character( IrrRep[k] ) then
                   flag[k] := true;
                   Append(~Ri, IrrRep[k]);
                   Append(~Vi, IrrMod[k]);
                   Append(~Xi, Character(IrrRep[k]) );
                end if;
            end for;
        end for;
        while i le #IrrRep and flag[i] do
            i := i+1;
        end while;
        Append(~X, Xi);
        Append(~Reps, Ri);
        Append(~Mods, Vi);
    end while;
    return Mods, Reps, X;
end intrinsic;




/* ****************************** ComputeKi ******************************
 *
 * ComputeKi(G::Grp, X::SeqEnum) -> List: Computes the character fields together with
 * additional data. For each Omega orbit X[i] we compute a List H[i] with
 *    H[i][1]    Q(chi), the field generated by the values of chi = X[i][1]
 *    H[i][2]    the character chi
 *    H[i][3]    iota : C -> K(chi), where C denotes the cyclotomic field used by Magma
 *               for the values of chi.
 *    H[i][4]    the maximal order of Q(chi)
 *
**************************************************************************** */


intrinsic ComputeKi(G::Grp, X::SeqEnum) -> List
{Computes the character fields. For each Q-rational character outputs
 [Kchi, chi, iota, maximal order], where iota: CoefficientField(chi) -> Kchi
}

    local H, i, chi, Qz, Qchi, i1;

    H := [**];

    for i in [1..#X] do
        chi :=  X[i,1];
        Qz := CoefficientField(chi);
        // Qz := SimpleExtension( CoefficientField(chi) );
        Qchi, i1  := sub< Qz | [Qz!chi[j] : j in [1..#chi]]>;
        H := Append(H, [* Qchi, chi, (i1^-1), MaximalOrder(Qchi) *]);
     end for;
     return H;
end intrinsic;


intrinsic ComputeCharacterFields(QG::AlgGrp) -> List
{Computes the character fields. For each Q-rational character outputs
 [Kchi, chi, iota, maximal order], where iota: CoefficientField(chi) -> Kchi
}

    local H, i, chi, Qz, Qchi, i1;

    H := [**];

    for i in [1..#QG`X] do
        chi :=  QG`X[i,1];
        Qz := CoefficientField(chi);
        Qchi, i1  := sub< Qz | [Qz!chi[j] : j in [1..#chi]]>;
        H := Append(H, [* Qchi, chi, (i1^-1), MaximalOrder(Qchi) *]);
     end for;
     return H;
end intrinsic;



/* ****************************** QIdempotents ******************************
 * 
 * QIdempotents(QG::AlgGrp, X::SeqEnum, H::List) -> SeqEnum: Computes the
 * Q-rational idempotents as a SeqEnum of AlgGrpElt's.
 *
****************************************************************************** */

intrinsic QIdempotents(QG::AlgGrp, X::SeqEnum, H::List) -> SeqEnum 
{Computes the Q-rational idempotents}

    local G, K, i, chi, psi, Kchi, i1, i2, v, e, E;

    E := [];
    K := CoefficientRing(QG);
    G := Group(QG);
    for i in [1..#X] do
        chi := X[i,1];
        psi := &+X[i];
        Kchi := H[i,1];
        iota := H[i,3];
        v := [K ! iota(psi(g^-1)) : g in G];
        v := [Degree(chi) / #G * a : a in v];
        e := QG ! v;
        Append(~E, e);
    end for;
    return E;
end intrinsic;

/*
intrinsic QIdempotents(QG::AlgGrp) -> SeqEnum
{Computes the Q-rational idempotents}

    local G, K, i, chi, psi, Kchi, i1, i2, v, e, E;

    E := [];
    K := CoefficientRing(QG);
    G := Group(QG);
    for i in [1..#QG`X] do
        chi := QG`X[i,1];
        psi := &+QG`X[i];
        Kchi := QG`H[i,1];
        iota := QG`H[i,3];
        v := [K ! iota(psi(g^-1)) : g in G];
        v := [Degree(chi) / #G * a : a in v];
        e := QG ! v;
        Append(~E, e);
    end for;
    return E;
end intrinsic;
*/

/* ****************************** ComputeB ******************************
 *
 * ComputeB(QG :: AlgGrp, H :: List) -> List: Computes data used by the
 * faster version of Psii.
 *
****************************************************************************** */


intrinsic ComputeB(QG :: AlgGrp, H :: List) -> List
{}
    local Bs, i, chi, Kchi, iota, B, g;

    Bs := [* *];
    for i:=1 to #H do
        chi := H[i,2];
        Kchi := H[i,1];
        iota := H[i,3];

        B :=[];
        for g in Group(QG) do
            Append(~B, ElementToSequence( iota(chi(g)) ));
        end for;
        Append(~Bs, Matrix(B));
    end for;
    return Bs;
end intrinsic;
        

/* ****************************** Psii ******************************
*
* Psii(x::FldNumElt, i::RngIntElt, QG::AlgGrp) -> AlgGrpElt: a much faster version
* Psii which assumes that QG`Psi and QG`B have already been computed. QG`Psi is computed with the
* function ComputePsi and QG`B with ComputeB.
*
************************************************************************ */

intrinsic Psii(x::FldNumElt, i::RngIntElt, QG::AlgGrp) -> AlgGrpElt
{Computes the image of x \in Kchi_i under the map Kchi_i --> Z(QG)}

    local chi, Kchi, K, iota, G, a, B, g, y, v, h, sum, j, t, psi;


    a := ElementToSequence(x);
    y:=Solution(QG`B[i], Vector(a));

    lambda := &+[y[j] * QG`Psi[i,j] : j in [1..Ncols(QG`Psi)]];
    return lambda;
end intrinsic;

/* ****************************** ComputePsi ******************************
*
* ComputePsi(QG::AlgGrp) -> Mtrx: Pre-computes data for Psii. For big examples
* this takes a lot of time.
*
************************************************************************ */

intrinsic ComputePsi(QG::AlgGrp) -> Mtrx
{}

    local chi, Kchi, K, iota, G, a, B, g, y, v, h, sum, j, t, psi;


    G := Group(QG);
    K := CoefficientRing(QG);

    Gelts := [g : g in G];
    
    M := [];
    for i:=1 to #QG`X do
        // print "In Zeile ", i, " von ", #QG`X;
        chi := QG`H[i,2];
        Kchi := QG`H[i,1];
        iota := QG`H[i,3];
        Zeile := [];
        for j:=1 to #Gelts do
	    v := [];
            for h in G do
                t := &+[psi(h^-1) * psi(Gelts[j]) : psi in QG`X[i]];
		Append(~v, K ! Degree(chi) / #G * iota(t));
            end for;
	    lambda := QG ! v;
	    Append(~Zeile, lambda);
	end for;
	Append(~M, Zeile);
    end for;

    return Matrix(M);
end intrinsic;


/* ****************************** Psi ************************************************
 *
 * Psi(x::FldNumElt, X::SeqEnum, H::List, QG::AlgGrp) -> AlgGrpElt: Computes the image
 * of x \in K(chi_1) \oplus ... \oplus K(chi_r) in Z(QG) under the map oplus_{i=1}^r Psii.
 *
************************************************************************************ */

intrinsic Psi(x::List, X::SeqEnum, H::List, QG::AlgGrp) -> AlgGrpElt
{ Computes the image of x \in Kchi_1 \opluss Kchi_r ---> Z(QG) }

    local i, lambda;
    
    lambda := QG ! 0;
    for i in [1..#x] do
	lambda := lambda + Psii(x[i], i, QG);
    end for;
    return lambda;
end intrinsic;


/* ****************************** QGCentre ************************************************
 *
 * QGCentre(QG::AlgGrp, X::SeqEnum, H::List) -> SeqEnum: Computes a SeqEnum of elements
 * of QG which constitutes a Z-basis of the maximal order in Z(QG).
 *
******************************************************************************************* */

intrinsic QGCentre(QG::AlgGrp, X::SeqEnum, H::List) -> SeqEnum
{ Computes a K-basis for the centre Z = Z(QG) of QG }
    local C, i, basis, b;

    C := [];
    for i in [1..#H] do
	basis := Basis(H[i,4]);
	for b in basis do
	    Append(~C, Psii(H[i,1] ! b, i, QG));
	end for;
    end for;
    return C;
end intrinsic;

intrinsic QGCentre(QG::AlgGrp) -> SeqEnum
{ Computes a K-basis for the centre Z = Z(QG) of QG }
    local C, i, basis, b;

    C := [];
    for i in [1..#QG`H] do
        basis := Basis(QG`H[i,4]);
        for b in basis do
            Append(~C, Psii(QG`H[i,1] ! b, i, QG));
        end for;
    end for;
    return C;
end intrinsic;



/* ****************************** CentralConductor ************************************************
 *
 * CentralConductor(QG::AlgGrp) -> Rec: Computes the central conductor
 * of ZG using [Curtis/Reiner I, (27.13)].
 *
******************************************************************************************* */

intrinsic CentralConductor(QG::AlgGrp) -> Rec
{ Computes a Z-basis for the central conductor of ZG }
    local C, i, basis, b, CentralF;

    H := QG`H;
    n := #Group(QG);
    C := [];
    for i in [1..#H] do
        basis := Basis(n / Degree(H[i,2]) * Different(H[i,4])^-1);
        for b in basis do
            Append(~C, Psii(H[i,1] ! b, i, QG));
        end for;
    end for;
    CentralF := ModuleInit (Matrix( [ ElementToSequence(C[i]) : i in [1..#C] ] ) );
    return CentralF;
end intrinsic;




/* ****************************** IstZentral ************************************************
 *
 * IstZentral(lambda::AlgGrpElt) -> BoolElt: Checks whether lambda is central
 *
******************************************************************************************* */
	
intrinsic IstZentral(lambda::AlgGrpElt) -> BoolElt
{ Returns true, if lambda is central }

    local g;

    for g in Group( Parent(lambda) ) do
	if g*lambda ne lambda*g then
	   return false;
	end if;
    end for;
    return true;
end intrinsic;


/* ****************************** CentralEltToBasis ******************************************
 * 
 * CentralEltToBasis(lambda::AlgGrpElt, Z::SeqEnum) -> SeqEnum: Given lambda \in Z(QG) this
 * computes the coefficients with respect to the Q-basis Z of Z(QG).

********************************************************************************************* */

intrinsic CentralEltToBasis(lambda::AlgGrpElt, Z::SeqEnum) -> SeqEnum
{ Returns the representation of the central element lambda with respect to
    the K-basis Z of Z(QG) }

    local a, B, z, y;

    assert IstZentral(lambda);
    a := Vector( ElementToSequence(lambda) );
    B := [];
    for z in Z do
	Append(~B, ElementToSequence(z));
    end for;
    B := Matrix(B);
    y := Solution(B, a);
    return y;
end intrinsic;


/* ****************************** CentralEltToBasis ******************************************
 *
 * PsiInv(lambda::AlgGrpElt, H::List, Z::SeqEnum) -> List: Computes the inverse of Psi.
 * It is important that Z is the Z-basis computed by QGCentre.
 *
********************************************************************************************* */

intrinsic PsiInv(lambda::AlgGrpElt, H::List, Z::SeqEnum) -> List
{ Computes the Inverse of Psi : Kchi_1 \oplus Kchi_r --> Z(QG) }

     local result, y, k, i, basis, alpha, b;

     assert IstZentral(lambda);
     result := [* *];
     y := CentralEltToBasis(lambda, Z);
     k := 1;
     for i in [1..#H] do
	 basis := Basis(H[i,4]);
	 alpha := 0;
	 for b in basis do
	     alpha := alpha + y[k]*( H[i, 1] ! b);
	     k := k+1;
	 end for;
	 Append(~result, alpha);
     end for;
     return result;
end intrinsic;

/* ****************************** DihedralS ******************************************
 *
 * DihedralS(G::Grp, H::List) -> AlgMatElt: We consider the Wedderburn decomposition as
 * described in CR: 
 *        Psi : QG --> oplus_{i=1}^r explicitly given matrix rings. 
 * Let w_1, ..., w_n be a "canonically" ordered Q-basis of this sum of matrix rings.
 * Then this function computes the matrix S such that
 *
 *    Psi( (g_1, ..., g_n)~ ) = S * (w_1, ..., w_n)~
 *
****************************************************************************************/


intrinsic DihedralS(G::Grp, H::List) -> AlgMatElt
{Computes the matrix representation of the Wedderburn decomposition }

    local v, sigma, tau, r, first_row, row, col, S, 
	  i, j, kchi, chi, iota, M, f, g, T, z;

    v := [g : g in Generators(G)];
    if Order(v[1]) eq #G/2 then
        sigma := v[1]; tau := v[2];
    else
        sigma := v[2]; tau := v[1];
    end if;

    r := #H;
    first_row := 1;
    S := ZeroMatrix(Rationals(), #G, #G);
    for i in [1..r] do
	Kchi := H[i, 1];
        chi := H[i, 2];
	iota := H[i, 3];

        /* Tsigma and Ttau define the Q-rational representation corresponding to chi */
        if Degree(chi) eq 1 then
	    Tsigma := Matrix(1,1, [ Rationals() ! chi(sigma) ]);
	    Ttau := Matrix(1,1, [ Rationals() ! chi(tau) ]);
	else
	    Tsigma := Matrix(2,2, [0,1,-1, iota(chi(sigma))] );
	    Ttau := Matrix(2,2, [1, 0, iota(chi(sigma)), -1]);
	end if;

        /* f is the Q-rational representation corresponding to chi */
        M := MatrixGroup<Integers() ! Degree(chi), Kchi | [Tsigma, Ttau]>;
	f := hom<G -> M | sigma -> Tsigma, tau -> Ttau>;

        col := 1;
	for g in G do
	    row := first_row;
	    T := ElementToSequence( f(g) ) ;
	    for z in T do
		v := ElementToSequence(z);
		for j in [1..#v] do
		    S[row, col] := Rationals() ! v[j];
		    row := row+1;
		end for;
	    end for;
	    col := col+1;
	end for;

        first_row := Integers() ! (first_row + Degree(chi)^2 * Degree(Kchi));
    end for;
    return S;
end intrinsic;

	   

/* ****************************** DihedralMaximalOrder ******************************************
 *
 * DihedralMaximalOrder(G::Grp, H::List) -> AlgMatElt: This computes the maximal order
 * naturally given by the results in CR: just replace the matrix rings over the fields by 
 * matrix rings over the corresponding maximal orders. 
 *
 * We represent the Z-basis of this maximal order in terms of the basis w_1, ..., w_n already 
 * used in DihedralS. Multiplying this matrix by S^-1 means applying the inverse of the "big"
 * Wedderburn isomorphism.
 *
********************************************************************************************** */

intrinsic DihedralMaximalOrder(G::Grp, H::List) -> AlgMatElt
{}

    local S, first_row, W, col, row, i, j, k, Kchi, chi, iota, d;


    S := DihedralS(G, H)^-1;

    first_row := 1;
    W := Matrix(#G, #G, [Rationals() ! 0 : i in [1..#G^2]]);
    col := 1;
    for i in [1..#H] do
        Kchi := H[i, 1];
        chi := H[i, 2];
        iota := H[i, 3];
        d := Integers() ! Degree(chi)^2;

        OKchi := MaximalOrder(Kchi);
        basis := Basis(OKchi);

        for j in [1..d] do
            row := first_row;
            for z in basis do
                v := ElementToSequence(z);
                for k in [1..#v] do
                    W[row, col] := Rationals() ! v[k];
                    row := row+1;
                end for;
                col := col+1;
                row := first_row;
            end for;
            first_row := Integers() ! first_row + Degree(Kchi);
        end for;
    end for;

    return Transpose( S*W );  /* Need a transpose because we always consider the rows as group
				 ring elements */
end intrinsic;


/* ****************************** AbelianMaximalOrder ******************************************
 *
 * AbelianMaximalOrder(QG :: AlgGrp) -> Rec: This computes the maximal order in QG when 
 * G is abelian.
 *
********************************************************************************************** */


intrinsic AbelianMaximalOrder(QG :: AlgGrp) -> Rec
{}

    local W, i, basis, z, lambda;


    W := [];
    for i in [1..#QG`H] do
        basis := Basis(MaximalOrder(QG`H[i,1]));

        for z in basis do
            lambda := Psii(QG`H[i, 1] ! z, i, QG);
            Append(~W, ElementToSequence(lambda));
        end for;
    end for;

    return ModuleInit( Matrix(W) );
end intrinsic;


/* ****************************** ComputeGi ******************************************
 *
 * ComputeGi(QG :: AlgGrp, CentralF :: Rec, Zbasis :: SeqEnum) -> List, List:
 * Computes  the List of central conductors [\frg_i : i = 1,...,r] and a list
 * [ InfinitePlaces_i : i = 1, ...., r ]. InfinitePlaces_i is the full list of all
 * embeddings, if the according component is quaternion. Otherwise it is an empty list.
 *
 * CentralF is the conductor \frF as a ZModule, Zbasis must be computed by QGCentre
 * in order to ensure that PsiInv works corectly.
 *
***************************************************************************************/

intrinsic ComputeGi(QG :: AlgGrp, CentralF :: Rec, Zbasis :: SeqEnum) -> List, List
{}

    local I, v, a, alpha, i, ArchPlaces, H, X, x, places, psi, CR;

    H := QG`H;
    X := QG`X;
    I := [* 0*H[i, 4] : i in [1..#H] *];
    for v in RowSequence(CentralF`hnf) do
        a := QG ! v;
        alpha := PsiInv(a, H, Zbasis);
        for i in [1..#H] do
           I[i] := I[i] + alpha[i]*H[i,4];
        end for;
    end for;
    ArchPlaces := [];
    for x in X do
        places := [];
	psi := AdamsOperation(x[1], 2);
	CR := Parent(psi);
	if InnerProduct(Id(CR), psi) eq -1 then
	    places := [i : i in [1..#x]];
	end if;
	Append(~ArchPlaces, places);
    end for;

    return I, ArchPlaces;
end intrinsic;

/* ****************************** ReducedNorm ******************************************
 *
 * ReducedNorm(QG :: AlgGrp, alpha :: AlgGrpElt) -> List: Computes a list of
 * ideals [* c_i : i = 1,...,r *], where c_i is the principal ideal generated
 * by nr(e_{chi_i}*alpha) in the i-th Wedderburn component of Z(QG).
 *
***************************************************************************************/

intrinsic ReducedNorm(QG :: AlgGrp, alpha :: AlgGrpElt) -> List
{}

    local i, redNorms, Kchi, iota, OKchi, chi, e, M, g, a, B, C, id;

    RedNorms := [* *];
    for i:=1 to #QG`H do
	Kchi := QG`H[i, 1];
	OKchi := QG`H[i, 4];
	chi := QG`H[i, 2];
	iota := QG`H[i, 3];
	KchiG := GroupAlgebra(Kchi, Group(QG));

	e := KchiG ! [Degree(chi) / #Group(QG) * iota(chi(g^-1)) : g in Group(QG)];
	M := Matrix( [ ElementToSequence(e*b) : b in Basis(KchiG) ] );
	M := EliminateZeroRows( EchelonForm(M) );
        

	lambda := [ KchiG ! a : a in RowSequence(M) ];
	a := KchiG ! [Kchi!c : c in Coefficients(alpha)];
	B := Matrix( [ ElementToSequence(a*b) : b in lambda  ]  );
	C := Solution(M, B);
	Id := Root( Determinant(C)*QG`H[i,4], Integers() ! Degree(chi) );
	Append(~RedNorms, Id);
    end for;

    return RedNorms;

end intrinsic;


/* ****************************** RelationMatrix ******************************************
 *
 * RelationMatrix(QG :: AlgGrp, AmodF :: SeqEnum, rcgp :: List) -> Mtrx: Given the ray
 * class groups in rcgp and a list of generators of (ZG / F)^times we compute the
 * relations given by these generators.
 *
***************************************************************************************/

intrinsic RelationMatrix(QG :: AlgGrp, AmodF :: SeqEnum, rcgp :: List) -> Mtrx
{}

    local i;

    R := [];
    orders := [];
    for i:=1 to #rcgp do
        // orders cat:= [Order(g) : g in Generators(rcgp[i,1])];
        orders cat:= [ Order(rcgp[i,1].j) : j in [1..NumberOfGenerators(rcgp[i,1])] ];
    end for;
    
    /* The trivial relations given by the orders of the generators of rcgp */
    if #orders ne 0 then
        D := DiagonalMatrix(orders);
    else
        D := ZeroMatrix(Integers(), 0, 0);
    end if;
	
    count := 1;
    vprintf RelKTh, 1: "We will compute %o relations\n", #AmodF;
    for alpha in AmodF do
        if count mod 10 eq 0 then
            vprintf RelKTh, 1: "%o relations computed\n", count;
        end if;
	row := [];
	Id := ReducedNorm(QG, alpha);
	for i:=1 to #Id do
	    f := rcgp[i, 2]^-1;
	    row cat:= ElementToSequence( f(Id[i]) );
	end for;
	Append(~R, row);
        count +:= 1;
    end for;

    return VerticalJoin(D, Matrix(R));
end intrinsic;


/* ****************************** RelationMatrix ******************************************
 *
 * RelationMatrix(QG :: AlgGrp, IgmodPg :: GrpAb, K :: GrpAb, m :: List, K1AmodF :: SeqEnum) -> GrpAb, GrpAb: 
 * Given the ray class groups in IgmodPg, K = ker(IgmodPg -> I1modP1) and a list of generators of K1AmodF 
 * we compute the relations given by these generators. As quotients we obtain cl(A) and D(A).
 *
***************************************************************************************/

intrinsic RelationMatrix(QG :: AlgGrp, IgmodPg :: GrpAb, K :: GrpAb, m :: List, K1AmodF :: SeqEnum) -> GrpAb, GrpAb, Map
{}

    local R, count, alpha, row, Id, i, clA, D, v;

    R := [];
    count := 1;
    vprintf RelKTh, 1: "We will compute %o relations\n", #K1AmodF;
    for alpha in K1AmodF do
        if not IsInvertible(alpha) then
            print "alpha not invertible !!!!!!!!!!!!!!!!!!!!!!!!  count = ", count;
        end if;
        if count mod 10 eq 0 then
            vprintf RelKTh, 1: "%o relations computed\n", count;
        end if;
        row := [];
        Id := ReducedNorm(QG, alpha);
        for i:=1 to #Id do
	    v := ElementToSequence( (m[i]) (Id[i]) );
	    if #v eq 0 then
	        v := [0];
	    end if;
            row cat:= v;
        end for;
        Append(~R, row);
        count +:= 1;
    end for;

    clA, f := quo< IgmodPg | [IgmodPg ! a : a in R] >;
    D := sub< clA |  [ f( IgmodPg ! gen) : gen in Generators(K) ] >;

    return clA, D, f;
end intrinsic;


/* ****************************** TestMetacyclic ******************************************
 *
 * TestMetacyclic(p :: RngIntElt, q :: RngIntElt, r :: RngIntElt) -> BoolElt
 * This function tests whether our computations agree with the result of
 * [CR, Theorem 50.25] for metacyclic groups. The function only compares the orders
 * of the groups which occur in [CR, Theorem 50.25].
 *
***************************************************************************************/

intrinsic TestMetacyclic(p :: RngIntElt, q :: RngIntElt, r :: RngIntElt) -> BoolElt
{}

    local G, clZG, C, z, beta, f, S, clS, d0, H, clZH;

    G := MetacyclicGroup(p, q, r);
    clZG := LocallyFreeClassgroup(G);
    C<z> := CyclotomicField(p);
    beta := &+[z^(r^i) : i in [1..q]];
    f := MinimalPolynomial(beta);
    S := MaximalOrder(f);
    clS := ClassGroup(S);
    if IsEven(q) then
        d0 := q div 2;
    else
        d0 := q;
    end if;
    H := CyclicGroup(q);
    clZH := LocallyFreeClassgroup(H);
    return  Order(clZG) eq d0*Order(clS)*Order(clZH);
end intrinsic;


/* ****************************** RayClassGroup ******************************************
 *
 * RayClassGroup(QG :: AlgGrp, Ideals :: List, InfinitePlaces :: SeqEnum) -> GrpAb, GrpAb, List
 * Given a list of ideals (\frakg_1, ..., \frakg_r) in Ideals and a list of infinite places
 * (\infty_1, ..., \infty_r) we compute the ray classgroup I_\frakg / P_\frakg^+ and the
 * kernel of the map I_\frakg / P_\frakg^+ --> I_1 / P_1^+.  
 * m is a list of maps which is used to compute the discrete logarithm in P_\frakg^+ --> I_1 / P_1^+.
 * In each component it is the inverse of the map returned by the MAGMA function RayClassGroup
 *
***************************************************************************************/

intrinsic RayClassGroup(QG :: AlgGrp, Ideals :: List, InfinitePlaces :: SeqEnum) -> GrpAb, GrpAb, List
{}
    local i, ordcl, ordcl1, m, images, cl, f, cl1, f1, Bilder, IgmodPg, I1modP1, s, v, pi, K;

    ordcl := [];
    ordcl1 := [];
    m := [* *];
    images := [];

    for i:=1 to #Ideals do
       if #InfinitePlaces[i] gt 0 then
            cl, f := RayClassGroup(Ideals[i], InfinitePlaces[i]);
            cl1, f1 := RayClassGroup(1*QG`H[i, 4], InfinitePlaces[i]);
        else
            cl, f := RayClassGroup(Ideals[i]);
            cl1, f1 := RayClassGroup(1*QG`H[i, 4]);
        end if;
	Append(~m, f^-1);
	if NumberOfGenerators(cl) eq 0 then
	    Append(~images, [ [0] ]);
	    Append(~ordcl, [ 1 ]);
	else
            if NumberOfGenerators(cl1) eq 0 then
                Append(~images, [[0] : j in [1..NumberOfGenerators(cl)] ] );
	    else
	        Append(~images, [ ElementToSequence( (f*(f1^-1))(cl.j) ) : j in [1..NumberOfGenerators(cl)] ]);
	    end if;
	    Append(~ordcl, [Order(cl.j) : j in [1..NumberOfGenerators(cl)]]); 
	end if;
        if NumberOfGenerators(cl1) eq 0 then
            Append(~ordcl1, [ 1 ]);
        else
            Append(~ordcl1, [Order(cl1.j) : j in [1..NumberOfGenerators(cl1)]]);
        end if;
    end for;
    IgmodPg := AbelianGroup( &cat ordcl );
    I1modP1 := AbelianGroup( &cat ordcl1 );
    s := 0;
    Bilder := [];
    for i:=1 to #Ideals do
        v := [0 : i in [1..NumberOfGenerators(I1modP1)]];
	for b in images[i] do
	    for j := 1 to #b do
	       v[s+j] := b[j];
	    end for;
	    Append(~Bilder, I1modP1 ! v);
	end for;
	s := s + #images[i,1];
    end for;
    pi := hom< IgmodPg -> I1modP1 | Bilder >;
    K := Kernel(pi);

    return IgmodPg, K, m;
end intrinsic;


MAX_ENUMERATION := 10000;
TRIALS := 5;

/*
*       NOTE: The notation in this file is bad. ZG may be any order in QG, not necessarily the
*       integral group ring.
*/



/* ****************************Oumf*************************************
*
*       Oumf(QG::AlgGrp, ZG::Rec, F::Rec, Ideals::List, K1Version::BoolElt) -> SeqEnum, BoolElt: 
*       This function should usually be used with K1Version := true. K1Version := false is only
*       kept for testing.
*       If K1Version := false, this function
*       tries to compute a system of generators for (ZG / F)^\times, where F is
*       a full two-sided ideal in ZG,in our application always the conductor.
*       The computation follows the approach described in the fax to Robert.
*       If G is not abelian and the pieces ZG / (\frp ZG + F) are two big, i.e. > MAX_ENUMERATION, then
*       we only randomly choose TRIALS many elements. In this case the functions returns
*       false as its second return value. If G is abelian, we know that ZG / (\frp ZG + F) is a field,
*       so we can find the generator with a probabilistic algorithm.
*
*       If K1Version := true, then we compute generators for K_1(ZG / F) as described
*       in the paper. This always returns an exact answer.
*
*       Note: the notation here is bad. ZG may be any order in QG, not necessarily the
*       integral group ring.
*
********************************************************************** */

intrinsic Oumf(QG::AlgGrp, ZG::Rec, F::Rec, Ideals::List, K1Version :: BoolElt) -> SeqEnum, BoolElt
{}
    local  i, j, k, m, n, r, primefact, ei, P, P_i,Ps,erz,  ps, qs, pj, prodaqj, 
           Amod, p_pot, neu, AmodF, bj, zz, Fcon,
           H, X, OC, exact, ex;

    H := QG`H; X := QG`X; OC := QG`OC;  // Faulheit, wieder aendern
        
    r := #Ideals; n := #Group(QG);
    s := 0;
    Ps := [];    eis := [];

    /* Compute a list of matrices P_i, one for each Wedderburn component, such
       that the rows of P_i represent the maximal order in the centre of QG.
    */
    erz := [**];  ei := [* H[i,1] ! 0 : i in [1..r] *];
    for j in [1..r] do
	basis := Basis(H[j,4]);
	P_i := ZeroMatrix(Rationals(), #basis, #Group(QG));
	row := 1;
	for k in [1..#basis] do
	    s := s+1;
	    ei[j] := H[j,1] ! basis[k];
            v := ElementToSequence( Psi(ei, X, H, QG) );
	    for col in [1..n] do
		P_i[row, col] := v[col];
	    end for;
	    row := row + 1;
	end for;
	ei[j] := H[j,1] ! 0;
	Append(~erz, P_i);
    end for;


    /* Compute the list of ideals \frp \sseq Z(QG) together with exponents e(p) such that
       CentralF = \prod \frp^e(p).
    */

    for i in [1..r] do
	primefact := Factorization(Ideals[i]);
	P := ZeroMatrix(Rationals(), 0, n);
	P_i := ZeroMatrix(Rationals(), 0, n);
	ei := [* H[j,1] ! 0 : j in [1..r] *];
	for j in [1..r] do
            if j ne i then
	        P_i := VerticalJoin(P_i, erz[j]);	
	    end if;
	end for;

	for pe in primefact do
	    P := P_i;
	    basis := [ H[i,1] ! v : v in Basis(pe[1]) ];
	    for alpha in basis do
		ei[i] := alpha;
		v := Matrix(Rationals(), 1, n, ElementToSequence( Psi(ei, X, H, QG) ));
		P := VerticalJoin(P, v);
	    end for;
	    P := ModuleInit(P);
	    Append(~Ps, P);
	    Append(~eis, pe[2]);
	end for;
    end for;

    /* Compute the intersections with ZG, find those which coincide after intersection, and compute
       the corresponding primary ideals.
    */
    ps := []; qs := [];
    pj := ModuleIntersection(QG, Ps[1], ZG);
    Append(~ps, pj);
    Append(~qs, ModuleIntersection(QG, ZG, ModulePow(QG, Ps[1], eis[1])));

    for i in [2..#eis] do
	// print "i = ", i;
	pj := ModuleIntersection(QG, Ps[i], ZG);
	j := 1; k := #ps;
	while j le k do
	    // print "j = ", j;
	    if ModuleCompare(QG,ps[j],pj) then
		// print "equality!!!!!!!!!!!!!!";
		qs[j] := ModuleIntersection(QG, qs[j], ModulePow(QG, Ps[i], eis[i]), OC); 
	        j := k+2;
	    end if;
	    j := j+1;
	end while;
	if j eq k+1 then
	    Append(~ps, pj);
	    Append(~qs, ModuleIntersection(QG, ZG, ModulePow(QG, Ps[i], eis[i])));
	end if;
    end for;

    /* Qs is a list of two sided ideals of ZG such that

                  1) Each two of the list are coprime, i.e.their sum equals ZG
		  2) Their intersection = product = F
    */
    Qs := [ModuleSum(QG, F, ModuleProd(QG, q, ZG)) : q in qs];
    m := #ps;
    AmodF := [];
    
    /* Try to compute a system of generators for (ZG / F)^times.
    */
    exact := true;
    for j in [1..m] do
	k, p_pot := K_and_p_pots(QG, ps[j], qs[j], ZG);
        prodaqj, ex :=  GroupExt(QG, ZG, F, k, p_pot, K1Version);
        exact := exact and ex;
        bj := ZG;
        for i:=1 to m do
            if i ne j then
                bj := ModuleIntersection(QG,bj, Qs[i]);
            end if;
        end for;
        AmodF cat:= ChinRev(QG, ZG, Qs[j], bj, prodaqj);
    end for;

    AmodF := MakeEltsQGInvertible(QG, AmodF, F); /* It may happen that elements alpha in AmodF
                                                    (which by construction are invertible in ZG / F)
                                                    are not invertible in QG.
						    Since F is full, one always finds gamma in F
						    such that alpha + gamma is invertible in QG */
    return AmodF, exact, ps, qs, Qs;
end intrinsic;




/* ****************************Oumf*************************************
*
*       Oumf(QG::AlgGrp, ZG::Rec, F::Rec, Ideals::List,  p :: RngIntElt) -> SeqEnum, BoolElt: 
*
*       We compute generators for K_1(ZpG / F_p) as described
*       in the paper. For this, only the p-part of the ideals in Ideals is considered.
*
*       Note: the notation here is bad. ZG may be any order in QG, not necessarily the
*       integral group ring.
*
********************************************************************** */

intrinsic Oumf(QG::AlgGrp, ZG::Rec, F::Rec, Ideals::List, p :: RngIntElt) -> SeqEnum
{}
    local  i, j, k, m, n, r, primefact, ei, P, P_i,Ps,erz,  ps, qs, pj, prodaqj, 
           Amod, p_pot, neu, AmodF, bj, zz, Fcon,
           H, X, OC, exact, ex;

    H := QG`H; X := QG`X; OC := QG`OC;  // Faulheit, wieder aendern
        
    r := #Ideals; n := #Group(QG);
    s := 0;
    Ps := [];    eis := [];

    /* Compute a list of matrices P_i, one for each Wedderburn component, such
       that the rows of P_i represent the maximal order in the centre of QG.
    */
    erz := [**];  ei := [* H[i,1] ! 0 : i in [1..r] *];
    for j in [1..r] do
	basis := Basis(H[j,4]);
	P_i := ZeroMatrix(Rationals(), #basis, #Group(QG));
	row := 1;
	for k in [1..#basis] do
	    s := s+1;
	    ei[j] := H[j,1] ! basis[k];
            v := ElementToSequence( Psi(ei, X, H, QG) );
	    for col in [1..n] do
		P_i[row, col] := v[col];
	    end for;
	    row := row + 1;
	end for;
	ei[j] := H[j,1] ! 0;
	Append(~erz, P_i);
    end for;


    /* Compute the list of ideals \frp \sseq Z(QG) together with exponents e(\frp) such that
       CentralF_p = \prod \frp^e(\frp).
    */

    for i in [1..r] do
	primefact := Factorization(Ideals[i]);
	P := ZeroMatrix(Rationals(), 0, n);
	P_i := ZeroMatrix(Rationals(), 0, n);
	ei := [* H[j,1] ! 0 : j in [1..r] *];
	for j in [1..r] do
            if j ne i then
	        P_i := VerticalJoin(P_i, erz[j]);	
	    end if;
	end for;

	for pe in primefact do
	    if IsPowerOf( #quo< H[i][4] | pe[1] >, p) then
	        P := P_i;
	        basis := [ H[i,1] ! v : v in Basis(pe[1]) ];
	        for alpha in basis do
		    ei[i] := alpha;
		    v := Matrix(Rationals(), 1, n, ElementToSequence( Psi(ei, X, H, QG) ));
		    P := VerticalJoin(P, v);
	        end for;
	        P := ModuleInit(P);
	        Append(~Ps, P);
	        Append(~eis, pe[2]);
	    end if;
	end for;
    end for;

    /* Compute the intersections with ZG, find those which coincide after intersection, and compute
       the corresponding primary ideals.
    */
    ps := []; qs := [];
/*    pj := ModuleIntersection(QG, Ps[1], ZG);
    Append(~ps, pj);
    Append(~qs, ModuleIntersection(QG, ZG, ModulePow(QG, Ps[1], eis[1])));
*/

    for i in [1..#eis] do
	// print "i = ", i;
	pj := ModuleIntersection(QG, Ps[i], ZG);
	j := 1; k := #ps;
	while j le k do
	    // print "j = ", j;
	    if ModuleCompare(QG,ps[j],pj) then
		// print "equality!!!!!!!!!!!!!!";
		qs[j] := ModuleIntersection(QG, qs[j], ModulePow(QG, Ps[i], eis[i]), OC); 
	        j := k+2;
	    end if;
	    j := j+1;
	end while;
	if j eq k+1 then
	    Append(~ps, pj);
	    Append(~qs, ModuleIntersection(QG, ZG, ModulePow(QG, Ps[i], eis[i])));
	end if;
    end for;

    /* Qs is a list of two sided ideals of ZG such that

                  1) Each two of the list are coprime, i.e.their sum equals ZG
		  2) Their intersection = product = p-part of F
    */
    Qs := [ModuleSum(QG, F, ModuleProd(QG, q, ZG)) : q in qs];
    m := #ps;
    AmodF := [];
    
    /* Try to compute a system of generators for (ZG / F)^times.
    */
    exact := true;
    for j in [1..m] do
	k, p_pot := K_and_p_pots(QG, ps[j], qs[j], ZG);
        prodaqj, ex :=  GroupExt(QG, ZG, F, k, p_pot, true);
        exact := exact and ex;
        bj := ZG;
        for i:=1 to m do
            if i ne j then
                bj := ModuleIntersection(QG,bj, Qs[i]);
            end if;
        end for;
        AmodF cat:= ChinRev(QG, ZG, Qs[j], bj, prodaqj);
    end for;

    AmodF := MakeEltsQGInvertible(QG, AmodF, F); /* It may happen that elements alpha in AmodF
                                                    (which by construction are invertible in ZG / F_p)
                                                    are not invertible in QG.
                                                    Since F is full, one always finds gamma in F
                                                    such that alpha + gamma is invertible in QG. Since
                                                    F \sseq F_p, we are allowed to change alpha in
                                                    this way.
                                                 */
    return AmodF, exact, ps, qs, Qs;
end intrinsic;

intrinsic MakeEltsQGInvertible(QG :: AlgGrp, AmodF:: SeqEnum, F :: Rec) -> SeqEnum
{}
    
    local i, beta;

    for i:=1 to #AmodF do
        if not IsInvertible(AmodF[i]) then
	    print "AmodF[", i, "] not invertible!!!!!!!!!!!!!!!!!!";
	    beta := AmodF[i] + Random(QG, F);
	    while not IsInvertible(beta) do
	        beta :=  AmodF[i] + Random(QG, F);
	    end while;
	    AmodF[i] := beta;
	end if;
    end for;

    return AmodF;
end intrinsic;

	    

/* Returns true if n is a power of the rational prime p. */

intrinsic IsPowerOf(n :: RngIntElt, p :: RngIntElt) -> BoolElt
{}

    f := Factorization(n);
    return #f eq 1 and f[1][1] eq p;

end intrinsic;


/* This function computes data which will be used in GroupExt. q is a p-primary ideal in 
   the centre Z(ZG).
*/
intrinsic K_and_p_pots(QG :: AlgGrp, p :: Rec, q :: Rec, ZG :: Rec) -> RngIntElt, SeqEnum
{}
        local k, p_pot;

        k := 0;
        p_pot := [p];      // p^1 = p;

        while not ModuleIncl(QG, p_pot[k+1], q) do
                Append(~p_pot, ModulePow(QG, p_pot[k+1], 2, ZG));
                k := k+1;
        end while;

        Append(~p_pot,q);      // for computational reasons
        return k, p_pot;
end intrinsic;


/* **************************** GroupExt *************************************
*
*       GroupExt(QG :: AlgGrp, ZG :: Rec, F :: Rec, k :: RngIntElt, p_pot :: SeqEnum, K1Version::BoolElt) -> SeqEnum, BoolElt
*       If K1Version := false this function
*       computes (ZG / qZG+F)^times, where q = p_pot[#p_pot]. Here q is p-primary with p = p_pot[1].
*       If #(ZG / pZG+F) > MAX_ENUMERATE and G is not abelian, then we only compute TRIALS many units of
*       (ZG / pZG+F)^times and return false as second return value.
*
*       If K1Version := true this function computes K_1(ZG / qZG+F).
*       We use the exact sequence
*                (1+pZG+F)/(1+qZG+F) --> K_1(ZG/qZG+F) --> K_1(ZG/pZG+F) -- 0.
*
********************************************************************** */

intrinsic GroupExt(QG :: AlgGrp, ZG :: Rec, F :: Rec, k :: RngIntElt, p_pot :: SeqEnum, K1Version :: BoolElt) -> SeqEnum, BoolElt
{}
       local p, q, a, b, G, cs, gammas, U_as, ns, h, Amod, 
             g0, H, v, i, j, elt, pos, eins, exact;

       p := p_pot[1];
       q := p_pot[#p_pot];

       v := #Group(QG);

       /* Try to compute a system of generators for the right hand side of the ses 
       */

       if not K1Version then
            Amod, exact := AmodPx(QG, ZG, ModuleSum(QG, F, ModuleProd(QG, p, ZG)));
       else
           Amod := ComputeK1(QG, ZG, ModuleSum(QG, F, ModuleProd(QG, p, ZG)));
           exact := true;
       end if;
     
       if k eq 0 then 
	   return Amod, exact; 
	end if;

       a := 1; k +:= 1;

       /* Successively compute generators for the left hand side of the ses
       */
       while a lt k do
                b := Min(a+1, k);
                M1 := ModuleSum(QG, F, ModuleProd(QG, ZG, ModuleSum(QG, q, p_pot[a])));
                M2 := ModuleSum(QG, F, ModuleProd(QG, ZG, ModuleSum(QG, q, p_pot[b])));
                G, orders := MmodN(QG, M1, M2);
                Amod cat:= [ Id(QG) + G[i] : i in [1..#orders] | orders[i] gt 1 ];

                a +:= 1;

       end while;
       return Amod, exact;
end intrinsic;


/* **************************** ComputeK1 *************************************
*
*       ComputeK1(QG :: AlgGrp, ZG :: Rec, p :: Rec ) -> SeqEnum:
*       Computes K_1(ZG/p), where p is assumed to be full two-sided ideal in ZG
*       containing the rational prime c.
*
********************************************************************** */

intrinsic ComputeK1(QG :: AlgGrp, ZG :: Rec, p :: Rec ) -> SeqEnum
{}
    local cn, c, n, FG, I, ZGmodp, f, K1, a, gen;

    /* Compute R := ZG/p as an associative algebra */
    cn := Factorization( Integers() ! Determinant(p`hnf))[1];
    c := cn[1]; n := cn[2];
    FG := GroupAlgebra(GF(c), Group(QG));
    I := [ FG ! a : a in RowSequence(p`hnf) ];
    ZGmodp, f := quo< FG | I >;
    f := f^-1;

    assert #ZGmodp eq (Integers() ! Determinant(p`hnf));

    K1 := K1OfFiniteAlgAss(ZGmodp);

    return [ Lift(QG, f(gen)) : gen in K1 ];

end intrinsic;


/* **************************** K1OfFiniteAlgAss *************************************
*
*       K1OfFiniteAlgAss(R :: AlgAss) -> SeqEnum: Computes generators of K_1 of an associative
*       algebra defined over a finite field GF(c).
*       We use the exact sequence
*                   1+J --> K_1(R) --> K_1(A/J) --> 0.
*
********************************************************************** */

intrinsic K1OfFiniteAlgAss(R :: AlgAss) -> SeqEnum
{}
    local J, RmodJ, f, g, K1OfRmodJ, K1, OneplusJ, Jsquared, Q;

    J := JacobsonRadical(R);
    RmodJ, g := quo< R | J >;

    K1OfRmodJ := K1OfFiniteSemisimpleAlgAss(RmodJ);
    K1 := [ (g^-1)(gen) : gen in K1OfRmodJ ];

    // print "Dimension of J = ", Dimension(J);

    /* Compute 1 + J */ 
    OneplusJ := [];
    f := map< J -> R | x :-> R!x >;
    while Dimension(J) gt 0 do
        Jsquared := J*J;
	Q, g := quo< J | Jsquared >;
	OneplusJ cat:= [ One(R) + ((g^-1)*f)(b) : b in Basis(Q) ]; 
        /* OneplusJ cat:= [ One(R) + (R!b) : b in Basis(Q) ]; */

	J := Jsquared;
	f := map< J -> R | x :-> R!x >;
    end while;

    return K1 cat OneplusJ;

end intrinsic;


/* **************************** K1OfFiniteSemisimpleAlgAss *************************************
*
*       K1OfFiniteSemisimpleAlgAss(Rss :: AlgAss) -> SeqEnum: Computes generators of K_1 of an 
*       associative semisimple algebra defined over a finite field GF(c).
*       We decompose Rss into simple components M_{n_i}(F_i) and use the isomorphism
*                F_i^\times \simeq K_1(F_i) \simeq K_1(M_{n_i}(F_i)
*
********************************************************************** */

intrinsic K1OfFiniteSemisimpleAlgAss(Rss :: AlgAss) -> SeqEnum
{}
    local Comp, C, Z, a, aR, K1;

    assert IsSemisimple(Rss);

    K1 := [];
    ZRss := Centre(Rss);
    CompZRss := PrimitiveIdempotents(ZRss);
    Comp := [ sub< Rss | [ b*(Rss!comp[2]) : b in Basis(Rss)]> : comp in CompZRss ];    
    idem := [ Rss ! comp[2] : comp in CompZRss ];
    
    for i:=1 to #Comp do
        Z := Centre(Comp[i]);
	a := Random(Z);
	while not IsUnit(a) do
	    a := Random(Z);
	end while;
	while Order(a) ne #Z-1 do
            a := Random(Z);
            while not IsUnit(a) do
                a := Random(Z);
            end while;
        end while;

        /* Let R := Comp[i]. Then R is isomorphic to a matrix ring M_n(Z). We have to
           lift the primitive element a in Z along
               R^* ->> K_1(R) -> K_1(M_n(Z)) -> K_1(Z) -> Z^*
           This means we must construct phi:M_n(Z) -> R. Then a lift of a \in Z^* is given, for example,
	   by phi(M), where M ist the identity matrix, the entry (1,1) replaced by a.
	   This will happen in the function K1Lift.
	*/
        
	aR := K1Lift(Comp[i], a);

        b := idem;
        b[i] := Rss!aR;
	Append(~K1, &+b);
    end for;

    return K1;
end intrinsic;


/* **************************** K1Lift *************************************
*
*       K1Lift(R :: AlgAss, a :: AlgAssElt) -> AlgAssElt: R is assumed to be a finite
*       simple associative algebra, and a is an element in the centre Z of R. Then
*       R is a matrix ring over Z. We construct an isomorphism phi:M_n(Z) -> R and return
*       the element phi(M(a)), where M(a) is the identity matrix with the entry at (1,1)
*       replaced by a.  
*
********************************************************************** */

intrinsic K1Lift(R :: AlgAss, a :: AlgAssElt) -> AlgAssElt
{}

    local i, Z, Fp, L, LFbasis, Rbasis, m,n, S, k, row, v, j, b, x;

    Z := Parent(a);
    Fp := BaseRing(Z);

    L := MinimalLeftIdeals(R : Limit:=1)[1];
    LFbasis := FindFBasis(L,Z);
    Rbasis := Basis(R);
    m := Dimension(R);
    n := #LFbasis;
    S := ZeroMatrix(Fp, m*n, m);
    for k:=1 to m do
        row := 1;
	for i:=1 to n do
	    v := ElementToSequence( Rbasis[k] * (R!LFbasis[i]) );
	    for j:=1 to m do
	        S[row, k] := v[j];
		row := row+1;
	    end for;
	end for;
    end for;

    b := ElementToSequence(R ! (a*LFbasis[1]));
    for i:=2 to n do
        b cat:= ElementToSequence(R ! LFbasis[i]);
    end for;
    b := Vector(b);
    x := Solution(Transpose(S), b);

    return R ! x;

end intrinsic;


/* **************************** FindFBasis *************************************
*
*       FindFBasis(L :: AlgAss, Z :: AlgAss) -> SeqEnum: L is assumed to be an R-module,
*       where R is a finite simple algebra over Fp with centre Z=F. L is given by
*       a Fp-basis. This function computes an F-basis of L. 
*
********************************************************************** */

intrinsic FindFBasis(L :: AlgAss, Z :: AlgAss) -> SeqEnum
{}
   
    local i, n, Zbasis, Lbasis, BasisOverF, M, b, N, MN, d;

    n := Integers() ! (Dimension(L) / Dimension(Z));
    Zbasis := Basis(Z);
    Lbasis := Basis(L);
    BasisOverF := [ Lbasis[1] ];
    M := Matrix( [ ElementToSequence(L!(b*Lbasis[1])) : b in Zbasis ] );
    d := 1;
    i := 2;
    while d lt n do
        N := Matrix( [ ElementToSequence(L!(b*Lbasis[i])) : b in Zbasis ] );
	MN := VerticalJoin(M, N);
	if Dimension(Kernel(MN)) eq 0 then
	    M := MN;
	    Append(~BasisOverF, Lbasis[i]);
	    d +:= 1;
	end if;
	i +:= 1;
    end while;
    return BasisOverF;
end intrinsic;

/* **************************** FindFBasis *************************************
*
*       FindFBasis(V :: Rec, M :: AlgMat, Q :: List, psi :: Map, R :: AlgAss) -> SeqEnum:
*       V is a F_p-module (given as a Rec FGModule (misuse of notation)), M = M_m(F) is a
*       matrix algebra and a factor of R, Q a list as computed by CentralSimpleAlgebraRep, 
*       psi is a map FpG -> R, R is a semisimple algebra over a finite field F_p.
*
*       V is given by a F_p-basis. We compute a $F$-basis of V.
*
********************************************************************** */

intrinsic FindFBasis(V :: Rec, M :: AlgMat, Q :: List, psi :: Map, 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( [ GroupRingAction(V, (psi^-1)(R!beta(M!b, Q)), Vbasis[1]) : b in Fbasis ] );

    d := 1;
    i := 2;
    while d lt n do
        N := Matrix( [ GroupRingAction(V, (psi^-1)(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;

/* **************************** AmodPx *************************************
*
*       AmodPx(QG :: AlgGrp, ZG :: Rec, p :: Rec) -> SeqEnum, BoolElt: Tries to
*       compute (ZG / p)^times, where p is a full twosided ideal in ZG which contains
*       a rational prime c, so that we can consider ZG / p as an Z/cZ-algebra.
*
********************************************************************** */

intrinsic AmodPx(QG :: AlgGrp, ZG :: Rec, p :: Rec) -> SeqEnum, BoolElt
{}
 
    if IsAbelian(Group(QG)) then
        return AbelianAmodPx(QG, ZG, p), true;
    end if;
 
    if Determinant(p`hnf) lt MAX_ENUMERATION then
        return EnumerateAmodPx(QG, ZG, p), true; 
    else
        return RandomAmodPx(QG, ZG, p), false; 
    end if;

end intrinsic;


/* **************************** AbelianAmodPx *************************************
*
*       AbelianAmodPx(QG :: AlgGrp, ZG :: Rec, p :: Rec) -> SeqEnum: In the abelian case
*       we know that ZG/p is a field, hence its multiplicative group is cyclic. We
*       randomly choose elements and compute the order.
*
********************************************************************** */

intrinsic AbelianAmodPx(QG :: AlgGrp, ZG :: Rec, p :: Rec) -> SeqEnum
{}
    local cn, c, n, FG, I, ZGmodp, f, d, a;

    print "In AbelianAmodPx Index = ", Factorization( Integers() ! Determinant(p`hnf));

    cn := Factorization( Integers() ! Determinant(p`hnf))[1];
    c := cn[1]; n := cn[2];
    FG := GroupAlgebra(GF(c), Group(QG));
    I := [ FG ! a : a in RowSequence(p`hnf) ];
    ZGmodp, f := quo< FG | I >;
    f := f^-1;

    d := c^n - 1;
    a := Random(ZGmodp);
    while not IsUnit(a) do a:= Random(ZGmodp); end while;
    while ord(a) ne d do
        a := Random(ZGmodp);
	while not IsUnit(a) do a:= Random(ZGmodp); end while;
    end while;

    return [ Lift( QG, f(a) ) ];
end intrinsic;

/* Compute the order of an invertible algebra element, assuming the order to be finite and small.
   Of course, this could be done faster (in particular when used by AbelianAmodPx, where we know
   an upper bound d for the order.
*/
intrinsic ord(a :: AlgAssElt) -> RngIntElt
{}
     local o, b;

     o := 1;
     b := a;
     while b ne Id(Parent(a)) do
         b := a * b;
	 o := o + 1;
     end while;
     
     return o;
end intrinsic;


/* **************************** RandomAmodPx *************************************
*
*       RandomAmodPx(QG :: AlgGrp, ZG :: Rec, p :: Rec) -> SeqEnum: Randomly choose
*       units in $ZG / p$.
*
********************************************************************** */

intrinsic RandomAmodPx(QG :: AlgGrp, ZG :: Rec, p :: Rec) -> SeqEnum
{}
    local i;

    print "In RandomAmodPx Index = ", Factorization( Integers() ! Determinant(p`hnf));

    cn := Factorization( Integers() ! Determinant(p`hnf))[1];
    c := cn[1]; n := cn[2];
    FG := GroupAlgebra(GF(c), Group(QG));
    I := [ FG ! a : a in RowSequence(p`hnf) ];
    ZGmodp, f := quo< FG | I >;
    f := f^-1;

    units := [];
    for i:= 1 to TRIALS do
        a := Random(ZGmodp);
        while not IsUnit(a) do a:= Random(ZGmodp); end while;
        Append( ~units, Lift(QG, f(a)) );
    end for;
    return units;
end intrinsic;


/* **************************** EnumerateAmodPx *************************************
*
*       EnumerateAmodPx(QG :: AlgGrp, ZG :: Rec, p :: Rec) -> SeqEnum: Computes a system
*       of generators for (ZG / p)^times by simple enumeation. This is very time and 
*       memory consuming. This is very,very ad hoc, there should by much smarter ways
*       to do this enumeration in order to get a smaller number of generators.
*
********************************************************************** */

intrinsic EnumerateAmodPx(QG :: AlgGrp, ZG :: Rec, p :: Rec) -> SeqEnum
{}
        local cn, c, n, FG, I, ZGmodp, f, orders, V, U, flag, i, gen, k, q, lambda, v, b, gens;

        print "In EnumerateAmodPx Index = ", Factorization( Integers() ! Determinant(p`hnf));
        
        cn := Factorization( Integers() ! Determinant(p`hnf))[1];
	c := cn[1]; n := cn[2];
	FG := GroupAlgebra(GF(c), Group(QG));
	I := [ FG ! a : a in RowSequence(p`hnf) ];
	ZGmodp, f := quo< FG | I >;
	f := f^-1;

        /* Each element in ZG/p corresponds to an index i.
	   flag[i] = 0, element has not yet been seen
	             1, this is a generator
		     2, unit, but not a generator
		     -1,not a unit
	*/
        flag := [0 : i in [1..c^n]];
	orders := [c : i in [1..n]];
        V :=  [ Id(ZGmodp) ];  /* a list of elements that can be generated by the generators already found */
        flag[ SequenceToIndex( ElementToSequence(Id(ZGmodp) ), orders ) ] := 2;

        for i:=1 to #flag do
            if flag[i] eq 0 then                /* a new element that has to be tested */
                gen := ZGmodp ! IndexToSequence(i, orders);
                if IsUnit(gen) then             /* the element is a unit */
                   flag[i] := 1;
                   lambda := gen*gen;
                   U := [Id(ZGmodp), gen];     /* compute the cyclic subgroup generated by gen */
                   while lambda ne Id(ZGmodp) do
                       Append( ~U, lambda );
                       flag[SequenceToIndex( ElementToSequence(Id(ZGmodp) ), orders )] := 2;
                       lambda := lambda*gen;
                    end while;
                    h := #V;
                    for j:=1 to Min(h, Round(#flag / #U)) do    /* there may not even be a heuristic
		                                                   argument for Min(h, Round(#flag / #U) */
                        for lambda in U do
			    xi := V[j] * lambda;
                            index := SequenceToIndex( ElementToSequence( xi ), orders );
                            if index eq i then
                                Append(~V, xi);
                            end if;
                            if flag[index] eq 0  then
                                flag[index] := 2;
                                Append(~V, xi);
                            end if;
			    xi := lambda * V[j];
                            index := SequenceToIndex( ElementToSequence( xi ), orders );
                            if flag[index] eq 0 then
                                flag[index] := 2;
                                Append(~V, xi);
                            end if;
                        end for;
                    end for;

                 else
                    flag[i] := -1;  // not invertible
                 end if;
             end if;
         end for;
	 
         gens := [];
         for i in [1..#flag] do
             if flag[i] eq 1 then
                Append(~gens,  Lift( QG, f( ZGmodp ! IndexToSequence(i, orders) ) ));
             end if;
         end for;

         return gens, flag;
end intrinsic;



/* No longer needed */
intrinsic MmodN(QG :: AlgGrp,M :: Rec,N :: Rec) -> SeqEnum, SeqEnum
{}
        local X, G, m, i, orders, gammas;

        X := N`hnf * M`hnf^(-1);

        S,B,C := SmithForm(MakeMatrixIntegral(X, 1));

        X := RowSequence( C^-1 * M`hnf );
        gammas := []; orders := [];
	for i in [1..#X] do
	    // if S[i,i] gt 1 then
	        Append(~gammas, QG ! X[i]);
		Append(~orders, S[i,i]);
	    // end if;
	end for;


        return gammas, orders, M`inv * C;
end intrinsic;



/* this function makes the inverse of the isomorphism phi (page 8 of the paper with M.Endres) explicit.
  phi is given by the Chinese Remainder theorem */

intrinsic ChinRev(QG :: AlgGrp, ZG :: Rec, M1 :: Rec, M2 :: Rec, v :: SeqEnum) -> SeqEnum
{}
        local A, B, H, U, a, b;

        A := MakeMatrixIntegral( M1`hnf * ZG`inv );
        B := MakeMatrixIntegral( M2`hnf * ZG`inv );

        H, U := HermiteForm(VerticalJoin(A,B));
        U := Submatrix(U,1,1,1,Ncols(U));

	a := QG ! ElementToSequence( Submatrix(U,1,1,1,Nrows(A)) * A );
	b := QG ! ElementToSequence( Submatrix(U,1,Nrows(A)+1 ,1,Nrows(B)) * B );

        assert a+b eq Id(QG);

	return [ lambda*b + a : lambda in v ];
end intrinsic;


        	
 
/* ****************************** Kern ******************************
 *
 * Kern(A :: AlgAss, idem :: AlgAssElt, a :: AlgAssElt, FX :: RngUPol) -> RngUPol:
 * Computes the kernel of the ring homomorphism
 *		phi : F[X] -> A, f(X) \mapsto idem*f(a)
 * A is an associative algebra over a finite field F,
 * idem is assumed to be an idempotent of A, a an element of A and FX the polynomial
 * ring over the field F.
 *
 * ********************************************************************** */

intrinsic Kern(A :: AlgAss, idem :: AlgAssElt, a :: AlgAssElt, FX :: RngUPol) -> RngUPolElt
{}

    local phi, M, b, T, S, j;

    phi := hom<FX -> A | t :-> idem*Evaluate(t,a)>;
    
/*
    M := [];
    for b in Basis(A) do
        Append(~M, ElementToSequence(b));
    end for;
*/

    M := [ ElementToSequence(b) : b in Basis(A) ];
    
    M :=Matrix(M);

    T := [ElementToSequence(phi(1))];
    S := Kernel(Matrix(T));
    j := 1;
    while Dimension(S) eq 0 do
        Append(~T, ElementToSequence(phi(FX.1^j)));
	S := Kernel(Matrix(T));
	j := j+1;
    end while;

    return FX ! ElementToSequence(Basis(S)[1]);
end intrinsic;


/********************* pRadical *******************************************
 *
 * pRadical(QG :: AlgGrp, A :: Rec, p :: RngIntElt) -> Rec: Computes the p-radical
 * of the order A using [Fr,Algo.5.1, page 45].
 *
 ***************************************************************************/

intrinsic pRadical(QG :: AlgGrp, A :: Rec, p :: RngIntElt) -> Rec
{}

    local F, Abasis, T, b1, b2, Amodp, J, Jbasis, W, b, w;

    // print "In pRadical with p = ", p;

    F := GF(p);
    Abasis := [ QG ! w : w in RowSequence(A`hnf) ];
/*
    T := [];
    for b1 in Abasis do
	for b2 in Abasis do
	    Append(~T, ElementToSequence(b1*b2));
	end for;
    end for;
*/
 
    T := [ ElementToSequence(b1*b2) : b1 in Abasis, b2 in Abasis ];

    T := Matrix(T) * A`inv;
    // print "    Generate Amodp";
    Amodp := AssociativeAlgebra< F, #Group(QG) | ElementToSequence(T) : Check := false >;
    // print "    Compute the Jacobson radical";
    J := JacobsonRadical(Amodp);
    // print "    J computed";
    Jbasis := [Amodp ! b : b in Basis(J)];

    W := [];
    for b in Jbasis do
        w := &+[ (Integers() ! b[i]) * Abasis[i] : i in [1..#Abasis] ];
	Append(~W, ElementToSequence(w));
    end for;
    if #W ne 0 then
        W := VerticalJoin(p*A`hnf, Matrix(W));
    else 
        W := p*A`hnf;
    end if;
    // print "Return from pRadical";
    return ModuleInit(W);
end intrinsic;


/********************* pHeriditaryOrder *******************************************
 *
 * pHeriditaryOrder(QG :: AlgGrp, A :: Rec, p :: RngIntElt) -> Rec: Computes the 
 * p-heriditary order of A using [Fr, Algo. 4.12, page 42]. 
 *
 ***************************************************************************/

intrinsic pHeriditaryOrder(QG :: AlgGrp, A :: Rec, p :: RngIntElt) -> Rec
{}

    local Rad, Aher, Ol;

    // print "In pHeriditaryOrder with p = ", p;

    Rad := pRadical(QG, A, p);
    Aher := A;
    Ol := ModuleConductor(QG, Rad, Rad);
    
    while not ModuleIncl(QG, Ol, Aher) do
        Aher := Ol;
	Rad := pRadical(QG, Aher, p);
	Ol := ModuleConductor(QG, Rad, Rad);
    end while;

    // print "Return from pHerdidtaryOrder";
    return Aher;
end intrinsic;


/********************* pMaximalOrder *******************************************
 *
 * pMaximalOrder(QG :: AlgGrp, A :: Rec, p :: RngIntElt) -> Rec: Computes the
 * p-maximal order of A using [Fr, Algo. 3.16, page 33].
 *
 ***************************************************************************/

intrinsic pMaximalOrder(QG :: AlgGrp, A :: Rec, p :: RngIntElt) -> Rec
{}

    local p_maximal, Amax, Ps, P, Ol;

    // print "In pMaximalOrder with p = ", p;

    p_maximal := false;
    Amax := A;

    while not p_maximal do
        Ps := PrimesOverp(QG, Amax, p);
	// print "		There are ", #Ps, " primes above ", p;
	p_maximal := true;
	for P in Ps do
	    // print "		Compote Ol";
            Ol := ModuleConductor(QG, P, P);
	    // print "		Ol computed";
	    if not ModuleIncl(QG, Ol, Amax) then
	        Amax := Ol;
		p_maximal := false;
		break;
	    end if;
	end for;
    end while;
    // print "Return from pMaximalOrder";
    return Amax;
end intrinsic;


/********************* MaximalOrder *******************************************
 *
 * MaximalOrder(QG :: AlgGrp, A :: Rec) -> Rec: Computes a
 * maximal order of A using a variant of [Fr, Algo. 4.17, page 43].
 * 
 * One should compare the performance of Friedrich's algorithm and the implemented
 * version!!!!!
 *
 ***************************************************************************/

intrinsic MaximalOrder(QG :: AlgGrp, A :: Rec) -> Rec
{}
    local ps, Amax, p;
     
    // print "In MaximalOrder";

    ps := [a[1] : a in Factorization(#Group(QG))];

    Amax := HeriditaryOrder(QG, A);
    for p in ps do
        Amax := pMaximalOrder(QG, Amax, p);
    end for;

    // print "Return from MaximalOrder";

    return Amax;
end intrinsic;
 
/********************* HeriditaryOrder *******************************************
 *
 * HeriditaryOrder(QG :: AlgGrp, A :: Rec) -> Rec: Computes a
 * heriditary order over A using [Fr, Algo. 4.16, page 43].
 *
 ***************************************************************************/
    
intrinsic HeriditaryOrder(QG :: AlgGrp, A :: Rec) -> Rec
{}

    local ps, Aher, p;
    
    // print "In HeriditaryOrder";

    ps := [a[1] : a in Factorization(#Group(QG))];

    Aher := A;
    for p in ps do
        Aher := pHeriditaryOrder(QG, Aher, p);
    end for;

    // print "Return from HeriditaryOrder";
    return Aher;
end intrinsic;

/********************* PrimesOverp *******************************************
 *
 * PrimesOverp(QG :: AlgGrp, A :: Rec, p :: RngIntElt) -> SeqEnum: Computes
 * the list of two-sided prime ideals of the order A over the rational prime p
 * using [Fr, Algo. 5.23, page 52]. Note that two-sided primes are maximal by
 * [Fr, Th. 2.27] or [Re, Th. 22.3].
 *
 * The main problem ist the computation of the simple components of the centre
 * of Ahut = A / \sqrt(pA). Here we use a probabilistic algorithm do to Eberly,
 * which works well in practice.
 *
 ***************************************************************************/

intrinsic PrimesOverp(QG :: AlgGrp, A :: Rec, p :: RngIntElt) -> SeqEnum
{}
    local F, Abasis, T, w, b1, b2, Amodp, pRad, basis, b, Ahut, f, ZAhut, SimpleCompZAhu,
          SimpleCompAhut, Ws, i, C, Ps;

    // print "In PrimesOverp with p = ", p;
    F := GF(p);
    Abasis := [ QG ! w : w in RowSequence(A`hnf) ];
    T := [];
    for b1 in Abasis do
        for b2 in Abasis do
            Append(~T, ElementToSequence(b1*b2));
        end for;
    end for;

    T := Matrix(T) * A`inv;
    Amodp := AssociativeAlgebra< F, #Group(QG) | ElementToSequence(T) : Check := false >;

    pRad := pRadical(QG, A, p);
    basis := [ QG ! w : w in RowSequence(pRad`hnf) ];
    T := [];
    for b in basis do
        Append(~T, ElementToSequence(b));
    end for;

    T := RowSequence( Matrix(T) * A`inv );
    // print "     Compute Ahut";
    Ahut, f := quo< Amodp | [Amodp ! w : w in T] >;
    // print "     Compute ZAhut";
    ZAhut := Centre(Ahut);
    // print "ZAhut computed";
    SimpleCompZAhut := PrimitiveIdempotents(ZAhut);

    /* The simple components of Ahut correspond to the simple components of ZAhut in
       the obvious way. For each simple component comp of ZAhut, comp[2] is  the primitive
       central idempotent.
    */ 
    SimpleCompAhut := [ sub< Ahut | [ b*(Ahut!comp[2]) : b in Basis(Ahut)]> : comp in SimpleCompZAhut ];

    /* We have to lift the maximal ideals in Ahut along A -> Ahut as described in
       [Fr, Algo. 5.23 (4)]. Ws will be a list of matrices, 
       one matrix W for each component C,
       where the rows of W correspond to the lifts of the basis elements of C.
    */
    Ws := [* *];
    for i:=1 to #SimpleCompAhut do
        C := SimpleCompAhut[i];
        basis := [Amodp ! (f^-1)(b) : b in Basis(C)];
        W := [];
        for b in basis do
            w := &+[ (Integers() ! b[i]) * Abasis[i] : i in [1..#Abasis] ];
            Append(~W, ElementToSequence(w));
        end for;
        Append(~Ws, W);
    end for;
    
    /* After these preparations we lift the maximal ideals of Ahut */
    Ps := [];
    for i:=1 to #Ws do
        W := [];
	for j:=1 to #Ws do
	    if j ne i then
	        W cat:= Ws[j];
	    end if;
	end for;
	W := VerticalJoin(pRad`hnf, Matrix(W));
	Append(~Ps, ModuleInit(W));
    end for;
    // print "Return from PrimesOverp";
    return Ps;
    
end intrinsic;

/********************* PrimitiveIdempotents *******************************************
 *
 * PrimitiveIdempotents(B :: AlgAss) -> List: Here B is a semi-simple algebra
 * over a finite field F. This function computes the simple components of B using
 * a probabilistic algorithm of Wayne Eberly, Computations for Algeras and Group 
 * Representations, Phd Thesis,University of Toronto, 1989. The thesis is available
 * via the author's homepage.
 *
 * We return a list [* C_1, .., C_r *], where C_i = [* simple algebra, idempotent , phi *].
 * phi is the inclusion map from the simple algebra to B, the idempotents are NOT
 * the idempotents of the simple algebra, nor of B. They are only needed for computational
 * reasons. 
 * The idempotents of B are obtained by phi(One(simple component));
 *
 ***************************************************************************/

intrinsic PrimitiveIdempotents(B :: AlgAss) -> List
{}
    local SimpleBs, OtherBs, FX, f, a, s, e, h, i, fine_idem, Comp;

    // print "In PrimitiveIdempotents";

    FX := PolynomialRing(CoefficientField(B));
    SimpleBs := [* *]; OtherBs := [* [* B, One(B), map<B->B | x:->x> *]  *];
    while #OtherBs gt 0 do
        B := OtherBs[1][1];   idem := OtherBs[1][2];  phi := OtherBs[1][3];
	if IsSimple(B) then
	    Append(~SimpleBs, OtherBs[1]);
	    Remove(~OtherBs, 1);
	else
            assert IsSemisimple(B);
            f := [];
	    while #f le 1 do
	        a := Random(B);
	        s := Factorization(Kern(B, idem, a, FX));
	        f := [v[1]^v[2] : v in s];
	    end while;
            h := [FX ! 0 : z in f];
            e := [FX ! 0 : z in f];
            for i:= 1 to #f do
                e[i] := FX!1;
                // h[i] := ChineseRemainder(f, e); 
                                                /* In new Magma versions this has to be replaced
                                                   by ChineseRemainderTheorem(e, f)
                                                */
                h[i] := ChineseRemainderTheorem(e, f); 
                e[i] := FX!0;
            end for;
            fine_idem := [Evaluate(g, a) : g in h];
            Comp := [* *];
	    for e in fine_idem do
	        Bsub, psi := sub<B | [e*x : x in Basis(B)]>;
	        Append(~Comp, [* Bsub, e, psi*phi  *]);
	    end for;
	    OtherBs := Comp cat Remove(OtherBs, 1);
	end if;
     end while;
     // print "Return from PrimitiveIdempotents";
     return SimpleBs;
end intrinsic;


/* ****************************** ModuleInit ******************************
 *
 * ModuleInit(W :: Mtrx) -> Rec: Given a matrix W with rational coefficients
 * this functon computes the corresponding Z-Module.
 *
 * ********************************************************************** */

intrinsic ModuleInit(W :: Mtrx) -> Rec
{}
 
 local M, A, H, i, j;


    // print "In ModuleInit";
    M := rec<ZModule | >; 
    A, M`denom := MakeMatrixIntegral(W);
    A := LLL(A);               /* Added 13.10.2016 */
    H := EliminateZeroRows( HermiteForm(A) );

    M`hnf := ZeroMatrix(Rationals(), Nrows(H), Ncols(H));
    for i in [1..Nrows(H)] do
	for j in [1..Ncols(H)] do
	    M`hnf[i, j] := H[i, j] / M`denom;
	end for;
    end for;

    if Ncols(M`hnf) eq Nrows(M`hnf) then
        M`inv := ( M`hnf )^-1;
    end if;

    return M;
end intrinsic;



/* ****************************** ModuleHNF ******************************
 *
 * ModuleHNF(M :: Rec) -> Mtrx: just returns M`hnf, where M is a ZModule
 *
 * ********************************************************************* */

intrinsic ModuleHNF(M :: Rec) -> Mtrx
{}

    return M`hnf;
end intrinsic;


/* ****************************** ModuleProd ******************************
 *
 * ModuleProd(QG :: AlgGrp, M :: Rec, N :: Rec) -> Rec: Computes the product of 
 * two ZModules M and N, where M and N are assumed to be ZModules in QG.
 *
 * ********************************************************************** */

intrinsic ModuleProd(QG :: AlgGrp, M :: Rec, N :: Rec) -> Rec
{}

    local m, n, MN, v, w, a, b;

    m := Nrows(M`hnf); 
    n := Nrows(N`hnf); 
    MN := [];
    for v in RowSequence(M`hnf) do
	a := QG ! v;
        for w in RowSequence(N`hnf) do
	    b := QG ! w;
            Append(~MN, ElementToSequence(a*b));
	end for;
    end for;
    mat := Matrix(MN);
    MN := ModuleInit(mat);
    return MN;
end intrinsic;

/* ****************************** ModulePow ******************************
 *
 * ModulePow(QG :: AlgGrp, M :: Rec , n :: RngIntElt) -> Rec: Computes the nth power of
 * two ZModules M, where M is assumed to be a ZModule in QG. n must be > 0. 
 * Very naive!!
 *
 * ********************************************************************** */

intrinsic ModulePow(QG :: AlgGrp, M :: Rec , n :: RngIntElt) -> Rec
{}

    local M1, i;

    if n eq 0 then 
	return QG`MaxOrd;
    end if;

    if n eq 1 then 
	return M;
    end if;

    M1 := M;

    for i in [2..n] do
       M := ModuleProd(QG , M , M1);
    end for;

    return M;
end intrinsic;

/* ****************************** ModulePow ******************************
 *
 * ModulePow(QG :: AlgGrp, M :: Rec , n :: RngIntElt) -> Rec: Computes the nth power of
 * two ZModules M, where M is assumed to be a ZModule in QG. n must be >= 0. For n=0 we
 * return the module R, which makes sense if M is an ideal in R.
 * Very naive!!
 *
 * ********************************************************************** */

intrinsic ModulePow(QG :: AlgGrp, M :: Rec , n :: RngIntElt, R :: Rec) -> Rec
{}

    local M1, i;

    if n eq 0 then
        return R;
    end if;

    if n eq 1 then
        return M;
    end if;

    M1 := M;

    for i in [2..n] do
       M := ModuleProd(QG , M , M1);
    end for;

    return M;
end intrinsic;


/* ****************************** ModuleSum ******************************
 *
 * ModuleSum(QG :: AlgGrp , M :: Rec, N :: Rec) -> Rec: Compute the sum of two 
 * ZModules M and N. QG is not needed, but all these functions should have the same shape.
 *
 * ********************************************************************* */

intrinsic ModuleSum(QG :: AlgGrp , M :: Rec, N :: Rec) -> Rec
{}

    local A, S;

    A := VerticalJoin(M`hnf, N`hnf);
    S := ModuleInit(A);
    return S;
end intrinsic;


/* ****************************** ModuleScalarProd ******************************
 *
 * ModuleScalarProd(QG:: AlgGrp , c :: AlgGrpElt , M :: Rec) -> Rec: 
 * Computes the left scalar product c*M, where M is assumed to be a ZModule in QG.
 *
 * **************************************************************************** */

intrinsic ModuleScalarProd(QG:: AlgGrp , c :: AlgGrpElt , M :: Rec) -> Rec
{}

    local m, N, v, a, mat;

    m := Nrows(M`hnf);
    N := [];
    for v in RowSequence(M`hnf) do
        a := QG ! v;
        Append(~N, ElementToSequence(c*a));
    end for;
    mat := Matrix(N);
    N := ModuleInit(mat);
    return N;
end intrinsic;


/* ****************************** ModuleScalarProd ******************************
 *
 * ModuleScalarProd(QG:: AlgGrp , M :: Rec, c :: AlgGrpElt) -> Rec:
 * Computes the right scalar product M*c, where M is assumed to be a ZModule in QG.
 *
 * **************************************************************************** */

intrinsic ModuleScalarProd(QG:: AlgGrp , M :: Rec, c :: AlgGrpElt) -> Rec
{}

    local m, N, v, a, mat;

    m := Nrows(M`hnf);
    N := [];
    for v in RowSequence(M`hnf) do
        a := QG ! v;
        Append(~N, ElementToSequence(a*c));
    end for;
    mat := Matrix(N);
    N := ModuleInit(mat);
    return N;
end intrinsic;




/* ****************************** DualModule ******************************
 *
 * DualModule(QG :: AlgGrp , M :: Rec) -> Rec: Computes the dual module of M.
 *
 * We consider the non degenerate, symmetric bilinear form
 *
 * <?,?> : Q[G] x Q[G] -> Q
 *
 * with
 *
 *          { 1, if st=id
 *  <s,t> = {                   s,t in G
 *          { 0, otherwise
 *
 * ********************************************************************** */

intrinsic DualModule(QG :: AlgGrp , M :: Rec) -> Rec
{}

    local n, A, i, v, a, l, g, D, j, y, Mdual;

    assert Nrows(M`hnf) eq Ncols(M`hnf); 
    n := Nrows(M`hnf);
    A := ZeroMatrix(Rationals(), n, n);

    i := 1;
    for v in RowSequence(M`hnf) do
	a := QG ! v;
        l := 1;
	for g in Group(QG) do
	    A[i, l] := Coefficient(a, g^-1);
	    l := l+1;
	end for;
	i := i + 1;
    end for;

    v := Vector(Rationals(), [0 : i in [1..n]]);
    A := Transpose(A);
    D := ZeroMatrix(Rationals(), n, n);

    for j in [1..n] do
       v[j] := Rationals() ! 1;
       y := Solution(A, v);
       for k in [1..n] do
	   D[j,k] := y[k];
       end for;
       v[j] := Rationals() ! 0;
    end for;

    Mdual := ModuleInit(D);
    return Mdual;

end intrinsic;



/* ****************************** ModuleSection ******************************
 *
 * ModuleSection(QG :: AlgGrp , M :: Rec , N :: Rec ) -> Rec: Computes the intersection
 * of the ZModules M and N. Both M and N are assumed to be full Zmodules in QG. 
 *
 * ************************************************************************* */

intrinsic ModuleSection(QG :: AlgGrp , M :: Rec , N :: Rec ) -> Rec
{}

    local S, Mdual, Ndual;

    Mdual := DualModule(QG, M);
    Ndual := DualModule(QG, N);
    S := ModuleSum(QG , Mdual , Ndual);    
    S := DualModule(QG , S);              

    return(S);
end intrinsic;

/* ****************************** ModuleIntersection ******************************
 *
 * ModuleIntersection(QG :: AlgGrp , M :: Rec , N :: Rec ) -> Rec: Computes the intersection
 * of two ZModules using an adaption of [Cohen II, 1.5.1]. One of the modules is assumed to
 * have full rank. 
 * This replaces ModuleSection, since it is more general and probably faster.
 *
 * ************************************************************************* */

intrinsic ModuleIntersection(QG :: AlgGrp , M :: Rec , N :: Rec ) -> Rec
{}

    local d, A, B, C, H, H4, S, i, j;

    assert Nrows(M`hnf) eq #Group(QG) or Nrows(N`hnf) eq #Group(QG);
    d := LCM(M`denom, N`denom);
    if Nrows(M`hnf) eq #Group(QG) then
        A := MakeMatrixIntegral(M`hnf, d);
        B := MakeMatrixIntegral(N`hnf, d);
    else
	B := MakeMatrixIntegral(M`hnf, d);
	A := MakeMatrixIntegral(N`hnf, d);
    end if;


    C := VerticalJoin( HorizontalJoin(A, A),
		       HorizontalJoin(B, ZeroMatrix(Integers(), Nrows(B), Ncols(A)))
		     );
    H := HermiteForm(C);
    H4 := Submatrix(H, Nrows(A)+1, Ncols(A)+1, Nrows(B), Ncols(B));
    H4 := EliminateZeroRows( HermiteForm(H4) );
    S := ZeroMatrix(Rationals(), Nrows(H4), Ncols(H4));
    for i in [1..Nrows(S)] do
	for j in [1..Ncols(S)] do
	    S[i,j] := H4[i,j] / d;
	end for;
    end for;
    return ModuleInit(S);
end intrinsic;

/* ****************************** ModuleIntersection ******************************
 *
 * ModuleIntersection(QG :: AlgGrp , M :: Rec , N :: Rec, OC :: SeqEnum ) -> Rec: 
 * Computes the intersection
 * of two ZModules using an adaption of [Cohen II, 1.5.1]. 
 * Both module are assumed to be full submoudles in the centre of QG
 *
 * ************************************************************************* */

intrinsic ModuleIntersection(QG :: AlgGrp , M :: Rec , N :: Rec, OC :: SeqEnum ) -> Rec
{}

    local d, A, B, C, H, H4, S, i, j, A1, B1, D;

    assert Nrows(M`hnf) eq #OC and  Nrows(N`hnf) eq #OC;
    d := LCM(M`denom, N`denom);
    A := MakeMatrixIntegral(M`hnf, d);
    B := MakeMatrixIntegral(N`hnf, d);

    A1 := ZeroMatrix(Integers(), 0, #OC);
    for a in RowSequence(A) do
        A1 := VerticalJoin(A1, Matrix(Integers(), 1, #OC, 
                                      ElementToSequence( CentralEltToBasis(QG!a, OC) ) ));
    end for;
    B1 := ZeroMatrix(Integers(), 0, #OC);
    for a in RowSequence(B) do
        B1 := VerticalJoin(B1, Matrix(Integers(), 1, #OC, 
                                      ElementToSequence( CentralEltToBasis(QG!a, OC) ) ));
    end for;



    C := VerticalJoin( HorizontalJoin(A1, A1),
                       HorizontalJoin(B1, ZeroMatrix(Integers(), Nrows(B1), Ncols(A1)))
                     );
    H := HermiteForm(C);
    H4 := Submatrix(H, Nrows(A1)+1, Ncols(A1)+1, Nrows(B1), Ncols(B1));
    H4 := EliminateZeroRows( HermiteForm(H4) );
    S := ZeroMatrix(Rationals(), Nrows(H4), Ncols(H4));
    for i in [1..Nrows(S)] do
        for j in [1..Ncols(S)] do
            S[i,j] := H4[i,j] / d;
        end for;
    end for;
    D := Matrix(Rationals(), #OC, #Group(QG), [ElementToSequence(a) : a in OC]);
    return ModuleInit(S*D);
end intrinsic;



/* ****************************** ModuleCompare ******************************
 *
 * ModuleCompare(QG :: AlgGrp , M :: Rec , N :: Rec) -> BoolElt, BoolElt, BoolElt: 
 * Compares two ZModules M and N.
 * It return M = N, M sseq N, N sseq M.
 *
 * ************************************************************************* */

intrinsic ModuleCompare(QG :: AlgGrp , M :: Rec , N :: Rec) -> BoolElt
{}

    return M`hnf eq N`hnf, ModuleIncl(QG, M, N), ModuleIncl(QG, N, M); 

end intrinsic;



/* ****************************** ModuleIncl ******************************
 *
 * ModuleIncl(QG :: AlgGrp , M :: Rec , N :: Rec) -> BoolElt: 
 * Returns M sseq N.
 *
 * ********************************************************************** */

intrinsic ModuleIncl(QG :: AlgGrp , M :: Rec , N :: Rec) -> Rec
{}

 local S;

 S := ModuleSum(QG , M , N);
 return S`hnf eq N`hnf;

end intrinsic;



/* ****************************** ModuleIndex ******************************
 *
 * ModuleIndex(QG :: AlgGrp , M :: Rec , N :: Rec) -> FldRatElt: 
 * returns the index [M : N], where both M and N are assumed to be full ZModules.
 *
 * *********************************************************************** */

intrinsic ModuleIndex(QG :: AlgGrp , M :: Rec , N :: Rec) -> FldRatElt
{}

    return Determinant( N`hnf * M`hnf^-1 );

end intrinsic;


/* ****************************** ModuleConductor ******************************
 *
 * ModuleConductor(QG :: AlgGrp , M :: Rec , N :: Rec) -> Rec:
 * This computes the conductor
 * F :=  {lambda in QG | lambda(M) subseteq N}.
 *
 * We compute the conductor as follows:
 *   F* = MN* ==> F = (MN*)*
 * where * denotes the dual module.
 *
 * *************************************************************************** */

intrinsic ModuleConductor(QG :: AlgGrp , M :: Rec , N :: Rec) -> Rec
{}

 local F;

 F := ModuleProd(QG , M , DualModule(QG, N));
 F := DualModule(QG , F);
 return F;

end intrinsic;


/* ****************************** IsModuleElt ******************************
 *
 * IsModuleElt(QG :: AlgGrp, lambda :: AlgGrpElt, M :: Rec) -> BoolElt:
 * Returns true if lambda in M, otherwise false.
 *
 * *********************************************************************** */

intrinsic IsModuleElt(QG :: AlgGrp, lambda :: AlgGrpElt, M :: Rec) -> BoolElt
{}

    local v, S;

    v := ElementToSequence(lambda);
    v := Matrix(1, #v, v); 
    S := ModuleInit(VerticalJoin(M`hnf, v));
    return S`hnf eq M`hnf;

end intrinsic;


/* ****************************** IsOrder ******************************
 *
 * IsOrder(QG :: AlgGrp, A :: Rec) -> BoolElt:
 * Returns true if A is an order in QG, false otherwise.
 * It does not check whether A is full !!!!!!!!!!!!!!!!!!!!!!!!!
 *
 * *********************************************************************** */

intrinsic IsOrder(QG :: AlgGrp, A :: Rec) -> BoolElt
{}

    local basis, b1, b2;

    basis := [QG ! b : b in RowSequence(A`hnf)];
    for b1 in basis do
        for b2 in basis do
	    if not IsModuleElt(QG, b1*b2, A) then
	        return false;
	    end if;
	end for;
    end for;
    if not IsModuleElt(QG, Id(QG), A) then
        return false;
    end if;

    return true;
end intrinsic;
	        

/* ****************************** IsUnit ******************************
 *
 * IsUnit(QG :: AlgGrp, ZG :: Rec, F :: Rec, a :: AlgGrpElt) -> BoolElt:
 * Returns true if a is a unit of ZG mod F, otherwise false.
 *
 * *********************************************************************** */

intrinsic IsUnit(QG :: AlgGrp, ZG :: Rec, F :: Rec, a :: AlgGrpElt) -> BoolElt
{}

    return ModuleCompare(QG, ZG, ModuleSum(QG, F, ModuleScalarProd(QG, a, ZG)));

end intrinsic;


/* ****************************** Random ******************************
 *
 * Random(QG :: AlgGrp, M :: Rec) -> AlgGrpElt: Returns a random element in the module M.
 *
 * *********************************************************************** */

intrinsic Random(QG :: AlgGrp, M :: Rec) -> AlgGrpElt
{}

    local i, basis, coeff;

    basis := [ QG ! row : row in RowSequence(M`hnf) ];
    coeff := [ Random(10) : b in basis ];
    return QG ! &+[ coeff[i]*basis[i] : i in [1..#basis] ];
end intrinsic;




intrinsic MatrixDenominator(M :: Mtrx) -> RngIntElt
{}

    local d, i, j;

    d := 1;
    for i in [1..NumberOfRows(M)] do
	for j in [1..NumberOfColumns(M)] do
	    if M[i, j] ne 0 then
		d := LCM(d, Denominator(M[i,j]));
	    end if;
	end for;
    end for;
    return d;
end intrinsic;

intrinsic MakeMatrixIntegral(M :: Mtrx) -> Mtrx
{}

    local d, i,j, S;

    d := MatrixDenominator(M);
    S := ZeroMatrix(Integers(), NumberOfRows(M), NumberOfColumns(M));
    for i in [1..NumberOfRows(M)] do
        for j in [1..NumberOfColumns(M)] do
            S[i,j] := Integers() ! (d * M[i,j]);
	end for;
    end for;
    return S, d;
end intrinsic;

intrinsic MakeMatrixIntegral(M :: Mtrx, d :: RngIntElt) -> Mtrx
{}

    local i,j, S;

    S := ZeroMatrix(Integers(), NumberOfRows(M), NumberOfColumns(M));
    for i in [1..NumberOfRows(M)] do
        for j in [1..NumberOfColumns(M)] do
            S[i,j] := Integers() ! (d * M[i,j]);
        end for;
    end for;
    return S, d;
end intrinsic;



intrinsic RowIsZero(M :: Mtrx, i :: RngIntElt) -> BoolElt
{}

    local j;

    for j in [1..Ncols(M)] do
	if M[i,j] ne 0 then
	    return false;
	end if;
    end for;
    return true;
end intrinsic;


intrinsic EliminateZeroRows(H :: Mtrx) -> Mtrx
{}

    local p, i;

    p := Nrows(H);
    for i in [1..Nrows(H)] do
	if RowIsZero(H, i) then
	    p := i-1;
	    break;
	end if;
    end for;
    return Submatrix(H,1,1,p,Ncols(H));
end intrinsic;


intrinsic TupleToIndex(t :: Tup, orders :: SeqEnum) -> RngIntElt
{}

    local basis, i, index;

    basis := 1; index := 1;
    for i:=#orders to 1 by -1 do
	index +:= t[i]*basis;
	basis *:= orders[i];
    end for;
    return index;

end intrinsic;

intrinsic SequenceToIndex(t :: SeqEnum, orders :: SeqEnum) -> RngIntElt
{}

    local basis, i, index;

    basis := 1; index := 1;
    for i:=#orders to 1 by -1 do
        index +:= (Integers() ! t[i]) * basis;
        basis *:= orders[i];
    end for;
    return index;

end intrinsic;



intrinsic ElementToIndex(t :: AlgGrpElt, T :: Mtrx, orders :: SeqEnum) -> RngIntElt
{}

    local b, v;

    b := Vector(ElementToSequence(t)) * T;
    v := < Integers()!b[k] mod orders[k] : k in [1..#orders] >;
    return TupleToIndex(v, orders);

end intrinsic;


intrinsic IndexToSequence(ind :: RngIntElt, orders :: SeqEnum) -> SeqEnum
{}

    local s, i, t;

    s := [];
    i := ind - 1;
    for t in Reverse(orders) do
	Append(~s, i mod t);
        i := (i - (i mod t)) div t;
    end for;
    return Reverse(s);
end intrinsic;

intrinsic Lift(QG :: AlgGrp, a :: AlgGrpElt) -> AlgGrpElt
{}
    return QG ! [ Integers() ! x : x in ElementToSequence(a)];
    
end intrinsic;



intrinsic QuaternionGroup(m :: RngIntElt) -> GrpPerm
{}

    local G, x, index, i, y;

    G := Sym(4*m);
    x := [];
    for index:=1 to 4*m do
        if index gt 2*m then j := 1; else j:=0; end if;
	i := (index - j*2*m) - 1;
	i := (i+1) mod (2*m); 
	Append(~x, 1+i+2*m*j); 
    end for;
    x := G!x;

    y := [];
    for index:=1 to 4*m do
        if index gt 2*m then j := 1; else j:=0; end if;
	i := (index - j*2*m) - 1;

	k := (j+1) mod 4;
	if k ge 2 then
	    i := (-i+m) mod (2*m);    j := k-2;
	else
	    i := (-i) mod (2*m);      j := k;
	end if;
	Append(~y, 1+i+2*m*j);
    end for;
    y := G!y;

    return sub< G | [x, y] >;
end intrinsic;

intrinsic GenDihedralGroup(H :: GrpPerm) -> GrpPerm, GrpPermElt, HomGrp
{}

    local G, elts, gens, g, U, D, b, j, tau;

    G := Sym(#H);
    elts := [g : g in H];
    gens := [ G ! [Index(elts, g*x) : x in elts] : g in Generators(H) ];
    U := sub< G | gens >;
    tau := G ! [Index(elts, x^(-1)) : x in elts];
    Append(~gens, tau);
    D := sub< G | gens >;

    b, j := IsIsomorphic(H, U);
    return D, tau, j;
end intrinsic;


intrinsic IsGenDihedral(G :: GrpPerm) -> BoolElt
{}
    local Us, v, U, GenDn;

    Us := NormalSubgroups(G : IsAbelian:=true, OrderEqual:=Integers()!(#G/2));
    for v in Us do
        U := v`subgroup;
        GenDn := GenDihedralGroup(U);
        if IsIsomorphic(G, GenDn) then
	    return true;
	end if;
    end for;
    return false;
end intrinsic;
  



intrinsic MetacyclicGroup(p :: RngIntElt, q :: RngIntElt, r :: RngIntElt) -> GrpPerm
{}

    local x, y, G, i, j;

    x := [];
    y := [];
    G := Sym(p*q);
    for j:=0 to q-1 do
        for i:=0 to p-1 do
	    Append(~x, j*p + ((i+1) mod p) + 1);
	    Append(~y, ((j+1) mod q)*p + ((i*r) mod p) + 1);
	end for;
    end for;
    x := G!x;
    y := G!y;

    return sub< G | [x, y] >;
end intrinsic;



intrinsic AdamsOperation(chi :: AlgChtrElt, n :: RngIntElt) -> AlgChtrElt
{}
    return (Parent(chi) ! [chi(c[3]^2) : c in Classes(Group(Parent(chi)))]);

end intrinsic;


intrinsic Order(a :: AlgAssElt) -> RngIntElt
{}

    local Eins, n, b;

    assert IsUnit(a);
    Eins := One(Parent(a));
    n := 1; b := a;
    while not b eq Eins do
        n +:= 1;
	b *:= a;
    end while;

    return n;
end intrinsic;


intrinsic IsDihedral(G :: Grp) -> BoolElt
{}
    local D;

    if not IsEven(#G) then
        return false;
    end if;

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




/* ****************************** ZGModuleInit ******************************
 *
 * ZGModuleInit(W :: Mtrx, phi :: Map) -> Rec: Given a matrix W with rational coefficients
 * this functon computes the corresponding ZG-Module.
 *
 * ********************************************************************** */

intrinsic ZGModuleInit(W :: Mtrx, phi :: Map) -> Rec
{}
 
 local M, A, H, i, j;

    M := rec<ZGModule | >; 
    A, M`denom := MakeMatrixIntegral(W);
    A := LLL(A);               /* Added 13.10.2016 */
    // print "HNF computed";
    H := EliminateZeroRows( HermiteForm(A) );
    M`hnf := ZeroMatrix(Rationals(), Nrows(H), Ncols(H));
    for i in [1..Nrows(H)] do
        for j in [1..Ncols(H)] do
            M`hnf[i, j] := H[i, j] / M`denom;
        end for;
    end for;

    if Ncols(M`hnf) eq Nrows(M`hnf) then
        M`inv := ( M`hnf )^-1;
    end if;

    M`phi := phi;

    return M;
end intrinsic;


intrinsic ZGModuleInit(N :: Rec, phi :: Map) -> Rec
{}

    local M;

    M := rec<ZGModule | >;
    M`hnf := N`hnf;
    M`denom := N`denom;
    M`inv := N`inv;
    M`phi := phi;

    return M;
end intrinsic;




/* ****************************** ZGModuleHNF ******************************
 *
 * ZGModuleHNF(M :: Rec) -> Mtrx: just returns M`hnf, where M is a ZModule
 *
 * ********************************************************************* */

intrinsic ZGModuleHNF(M :: Rec) -> Mtrx
{}

    return M`hnf;
end intrinsic;


/* ****************************** ZGModuleSum ******************************
 *
 * ZGModuleSum(M :: Rec, N :: Rec) -> Rec: Compute the sum of two 
 * ZModules M and N. 
 *
 * ********************************************************************* */

intrinsic ZGModuleSum(M :: Rec, N :: Rec) -> Rec
{}

    local A, S;

    assert M`phi eq N`phi;

    A := VerticalJoin(M`hnf, N`hnf);
    S := ZGModuleInit(A, M`phi);
    return S;
end intrinsic;


/* ****************************** ZGModuleIntersection ******************************
 *
 * ZGModuleIntersection(QG :: AlgGrp , M :: Rec , N :: Rec ) -> Rec: Computes the intersection
 * of two ZModules using an adaption of [Cohen II, 1.5.1]. One of the modules is assumed to
 * have full rank. 
 * This replaces ModuleSection, since it is more general and probably faster.
 *
 * ************************************************************************* */

intrinsic ZGModuleIntersection(M :: Rec , N :: Rec ) -> Rec
{}

    local d, A, B, C, H, H4, S, i, j, n;

    assert M`phi eq N`phi;
    n := Dimension( Codomain(M`phi) );
    assert Nrows(M`hnf) eq n  or Nrows(N`hnf) eq n;
    d := LCM(M`denom, N`denom);
    if Nrows(M`hnf) eq n then
        A := MakeMatrixIntegral(M`hnf, d);
        B := MakeMatrixIntegral(N`hnf, d);
    else
	B := MakeMatrixIntegral(M`hnf, d);
	A := MakeMatrixIntegral(N`hnf, d);
    end if;


    C := VerticalJoin( HorizontalJoin(A, A),
		       HorizontalJoin(B, ZeroMatrix(Integers(), Nrows(B), Ncols(A)))
		     );
    H := HermiteForm(C);
    H4 := Submatrix(H, Nrows(A)+1, Ncols(A)+1, Nrows(B), Ncols(B));
    H4 := EliminateZeroRows( HermiteForm(H4) );
    S := ZeroMatrix(Rationals(), Nrows(H4), Ncols(H4));
    for i in [1..Nrows(S)] do
	for j in [1..Ncols(S)] do
	    S[i,j] := H4[i,j] / d;
	end for;
    end for;
    return ZGModuleInit(S, M`phi);
end intrinsic;



/* ****************************** ZGModuleCompare ******************************
 *
 * ZGModuleCompare(M :: Rec , N :: Rec) -> BoolElt, BoolElt, BoolElt: 
 * Compares two ZGModules M and N.
 * It returns M eq N, M sseq N, N sseq M.
 *
 * ************************************************************************* */

intrinsic ZGModuleCompare(M :: Rec , N :: Rec) -> BoolElt, BoolElt, BoolElt
{}

    if M`phi eq N`phi then
        return M`hnf eq N`hnf, ZGModuleIncl(M, N), ZGModuleIncl(N, M); 
    else
        return false,false, false;
    end if;

end intrinsic;



/* ****************************** ZGModuleIncl ******************************
 *
 * ZGModuleIncl(M :: Rec , N :: Rec) -> BoolElt: 
 * Returns M sseq N.
 *
 * ********************************************************************** */

intrinsic ZGModuleIncl(M :: Rec , N :: Rec) -> Rec
{}

 local S;

 S := ZGModuleSum(M , N);
 return M`phi eq N`phi and S`hnf eq N`hnf;

end intrinsic;



/* ****************************** ZGModuleIndex ******************************
 *
 * ZGModuleIndex(M :: Rec , N :: Rec) -> FldRatElt: 
 * returns the index [M : N], where both M and N are assumed to be full ZGModules.
 *
 * *********************************************************************** */

intrinsic ZGModuleIndex(M :: Rec , N :: Rec) -> FldRatElt
{}

    return Determinant( N`hnf * M`hnf^-1 );

end intrinsic;


/* ****************************** IsZGModuleElt ******************************
 *
 * IsZGModuleElt(lambda :: AlgGrpElt, M :: Rec) -> BoolElt:
 * Returns true if lambda in M, otherwise false.
 *
 * *********************************************************************** */

intrinsic IsZGModuleElt(lambda :: SeqEnum, M :: Rec) -> BoolElt
{}

    local v, S;

    v := ElementToSequence(lambda);
    v := Matrix(1, #v, v); 
    S := ZGModuleInit(VerticalJoin(M`hnf, v), M`phi);
    return S`hnf eq M`hnf;

end intrinsic;


/* Computes a map phi : G -> Gl_n(Q), n = #G, which is the regular representation with
   respect to the canonical basis v_1 = g_1, ..., v_n = g_n, G = \{g_1, ...., g_n \}.
*/

intrinsic RegularRep(QG :: AlgGrp) -> Map
{}
    local G, gens, T, g, t, h, Gl, S;

    G := Group(QG);

    gens := [g : g in Generators(G)];
    if #gens eq 0 then
        gens := [ Id(G) ];
    end if;
    T := [];
    for g in gens do
        t := [];
        for h in Basis(QG) do
            Append(~t, ElementToSequence(g * h));
        end for;
        Append(~T, Matrix(t));
    end for;

    Gl := GL(#G, Rationals());
    phi := hom<G->Gl | [gens[i]->Transpose(T[i]) : i in [1..#gens]]>;

    return phi;
end intrinsic;

/* Computes the Swan module < r, \sum_{g in G} g >_ZG. phi is assumed to be the regular
   representation (e.g. computed by RegularRep).
*/
intrinsic SwanModule(r :: RngIntElt, phi :: Map) -> Rec
{}
    local G, W, i, S, gens, T, g, t, v, Gl;

    n := Dimension(Codomain(phi));
    W := [[Rationals()!1 : i in [1..n]]];
    W := VerticalJoin(Matrix(W), ScalarMatrix(Rationals(), n, r));
    S := ZGModuleInit(W, phi);

    return S;
end intrinsic;


/* Computes the Swan module < r, \sum_{g in G} g >_ZG. phi is assumed to be the regular
   representation (e.g. computed by RegularRep).
*/
intrinsic DeltaGH(QG :: AlgGrp, H :: GrpPerm, l :: RngIntElt, phi :: Map) -> Rec
{}
    local n, W1, W2, W, S;

    n := Dimension(Codomain(phi));
    W1 := ScalarMatrix(Rationals(), n, l);
    W2 := [ElementToSequence( QG!(g*h) - QG!g ) : g in RightTransversal(Group(QG), H), h in H];
    W := VerticalJoin(Matrix(W1), Matrix(W2));
    S := ZGModuleInit(W, phi);

    return S;
end intrinsic;



/* Computes ZG as ZG-module. phi is assumed to be the regular
   representation (e.g. computed by RegularRep).
*/
intrinsic RegularModule(phi :: Map) -> Rec
{}
    local G, gens, T, g, t, h, Gl, S;

    n := Dimension(Codomain(phi));
    S := ZGModuleInit( ScalarMatrix(Rationals(), n, 1), phi);
    return S;
end intrinsic;



/* Computes the ZG-Module ZG*lambda. phi is assumed to be the regular
   representation (e.g. computed by RegularRep).
*/
intrinsic ZGlambda(phi :: Map, lambda :: AlgGrpElt) -> Rec
{}
    local G, W, g, S;

    G := Domain(phi);
    W := Matrix( [ ElementToSequence(g*lambda) : g in G] );
    S := ZGModuleInit( W, phi);
    return S;
end intrinsic;

intrinsic ZGSubModule(M :: Rec, gens :: SeqEnum) -> Rec
{}
    local W, S;

    W := Matrix( gens );
    S := ZGModuleInit( W, M`phi);
    return S;
end intrinsic;



/* Computes the ZG-module generated by the elements in gens.
*/
/*
intrinsic ZGModule(phi :: Map, gens :: SeqEnum) -> Rec
{}
    local G, L, x, g, S, W;

    G := Domain(phi);
    L := [];
    for x in gens do
        for g in G do
	    Append(~L, Matrix([x])*phi(g));
	end for;
    end for;
    W := Matrix(L);
    S := ZGModuleInit(W, phi);
    return S;
end intrinsic;

*/

/* Computes the FG-sub module generated by the elements in gens.
*/
intrinsic FGSubModule(S :: Rec, gens :: SeqEnum) -> Rec
{}
    local G, L, x, g, Q, ZG;

    G := Domain(S`phi);
    ZG := GroupAlgebra(Integers(), G);
    L := [];
    for x in gens do
        for g in G do
            Append(~L, GroupRingAction(S, ZG!g, S`M ! Matrix([x])));
        end for;
    end for;
    Q := rec<FGModule | >;
    Q`M := sub<S`M | L >;
    Q`phi := S`phi;
    return Q;
end intrinsic;


/* Computes the F-submodule of S generated by the elements in gens. Although we return
   a Rec FGModule, this is no FG-module in general.
*/
intrinsic FSubModule(S :: Rec, gens :: SeqEnum) -> Rec
{}
    local G, L, Q;

    G := Domain(S`phi);
    L := [];
    Q := rec<FGModule | >;
    Q`M := sub<S`M | gens >;
    Q`phi := S`phi;
    return Q;
end intrinsic;


/* S is a FG-module and J a left ideal in FG. We compute the FG-submodule J*S.
*/
intrinsic FGSubModule(S :: Rec, J :: AlgGrpSub) -> Rec
{}
    local G, gens, lambda, a,  Q;

    G := Domain(S`phi);

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

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


/* U is a FG-submodule of S. We compute the quotient module S/U.
*/
intrinsic FGQuotientModule(S :: Rec, U :: Rec) -> Rec, Map
{}
    local G, F, FG, Q, Gl, gens, g, T, t;

    G := Domain(S`phi);
    F := CoefficientField(S`M);
    FG := GroupAlgebra(F, G);

    Q := rec<FGModule | >;
    Q`M, tau := quo<S`M | U`M >;
    Gl := GL(Dimension(Q`M), F);

    gens := [g : g in Generators(G)];
    if #gens eq 0 then 
        gens := [ Id(G) ];
    end if;
    T := [];
    for g in gens do
        t := [ tau( GroupRingAction(S, FG!g, (tau^-1)(b)) ) : b in Basis(Q`M)   ];
        Append(~T, Matrix(t));
    end for;
    Q`phi := hom<G->Gl | [gens[i]->Transpose(T[i]) : i in [1..#gens]]>;
    
    return Q, tau;
end intrinsic;


/* Computes lambda(a) for a module element a. a is assumed to be a row vector representing
   an element of S with respect of the underlying vector space. This works for ZG-modules and
   FG-modules.
*/
intrinsic GroupRingAction(S :: Rec, lambda :: AlgGrpElt, a :: ModTupRngElt) -> ModTupRngElt
{}
    local G, g;
    
    G := Group(Parent(lambda));
    return &+[ Coefficient(lambda, g) * a * Transpose(S`phi(g)) : g in G];
end intrinsic;


intrinsic LinExtChar(chi :: AlgChtrElt, lambda :: AlgGrpElt) -> FldCycElt
{}
    local G, g;

    G := Group(Parent(lambda));
    return &+[Coefficient(lambda, g) * chi(g) : g in G];
end intrinsic;


/* Computes the direct sum of the ZGmodules M and N, both given as Rec ZGModule.
*/
intrinsic ZGModuleDirectSum(M :: Rec, N :: Rec) -> Rec
{}
    local G, D, Gl, g, W, phi;

    G := Domain(M`phi);
    
    n := Dimension(Codomain(M`phi)) + Dimension(Codomain(N`phi));
    W := DiagonalJoin(M`hnf, N`hnf);
    Gl := GL(n, Rationals());
    phi := hom<G->Gl | [ g -> DiagonalJoin(M`phi(g), N`phi(g)) : g in Generators(G) ]>;
    D := ZGModuleInit(W, phi);

    return D;
end intrinsic;


/* Here M is a Rec ZGModule, given by a Z-basis w_1, ..., w_n (represented by the rows of M`hnf.).
   We compute the FGModule F w_1 \oplus ... \oplus F w_n, where F = Z/pZ.
*/
intrinsic ZGtoFGModule(M :: Rec, p :: RngIntElt) -> Rec
{}
    local G, D, Gl, g, n, U, Uinv;

    G := Domain(M`phi);

    n := Dimension(Codomain(M`phi));
    U := Transpose( M`hnf );
    Uinv := U^-1;

    D := rec<FGModule | >;
    D`M := RSpace(GF(p), n);
    Gl := GL(n, GF(p));
    gens := Generators(G);
    if #gens eq 0 then
        gens := [ Id(G) ];
    end if;
    D`phi := hom<G->Gl | [ g -> Gl ! MakeMatrixIntegral(Uinv*M`phi(g)*U) : g in gens ]>;

    return D;
end intrinsic;


    
/* Lifts an element v in the FG-module D/pD to an element in D.
*/
intrinsic LiftFGModuleElt(D :: Rec, v :: ModTupFldElt) -> ModTupRngElt
{}
   local w; 
   
   w := Vector(Integers(), v);
   return Vector(Rationals(), w)*D`hnf;
end intrinsic;




/* ****************************** IsElementary ******************************
 *
 * IsElementary(G::Grp) -> BoolElt: Returns true if G is elementary. By definition
 * this means that G is p-elementary for a prime p dividing #G.
 *
****************************************************************************** */

intrinsic IsElementary(G :: Grp) -> BoolElt
{}
    local p;

    if #G eq 1 then
       return true;
    end if;
    for p in PrimeDivisors(#G) do
        if IspElementary(G, p) then
	    return true;
	end if;
    end for;
    return false;
end intrinsic;


/* ****************************** IspElementary ******************************
 *
 * IspElementary(G :: Grp, p :: RngIntElt) -> BoolElt: Returns true if G is p-elementary. 
 * By definition this means that G is the direct product of a p-group (which may be trivial) and
 * a cyclic group of order prime to p. 
 *
****************************************************************************** */


intrinsic IspElementary(G :: Grp, p :: RngIntElt) -> BoolElt
{}
    local P, f, Q;

    if IsCyclic(G) then
        return true;
    end if;
    if not IsNilpotent(G) then
       return false;
    end if;

/*    
    for p in PrimeDivisors(#G) do
        P := SylowSubgroup(G, p);
	f, Q := HasComplement(G, P);
	if IsCyclic(Q) then
	    return true;
	end if;
    end for;
*/
    P := SylowSubgroup(G, p);
    f, Q := HasComplement(G, P);
    if IsCyclic(Q) then
        return true;
    end if;

    return false;
end intrinsic;



/* ****************************** ElementarySubgroups ******************************
 *
 * ElementarySubgroups(G :: GrpPerm) -> SeqEnum: Returns the list of elementary
 * subgroups of G.
 *
****************************************************************************** */

intrinsic ElementarySubgroups(G :: GrpPerm) -> SeqEnum
{}
    local ES, S, s;

    ES := [];
    S := Subgroups(G);
    for s in S do
       if IsElementary(s`subgroup) then
           Append(~ES, s`subgroup);
       end if;
    end for;
    return ES;
end intrinsic;


/* ****************************** BrauerInduction ******************************
 *
 * BrauerInduction(chi :: AlgChtrElt) -> SeqEnum: Returns a list of triples
 * [U, phi, c], where U is an elementary subgroup of G, phi is an abelian 
 * character of U, and c_{U,phi} is a multiplicity. One has
               chi = sum_{U,phi} c_{U,phi} ind_U^G(phi)
 *
****************************************************************************** */

intrinsic BrauerInduction(chi :: AlgChtrElt) -> SeqEnum
{}

    local G, ES, CTs, Y, s, i, psi, TG, A, j, y, b, c;

    G := Group(Parent(chi));

    if Degree(chi) eq 1 then
        return [ [* G, chi, 1 *] ];
    end if;

    ES := ElementarySubgroups(G);
    CTs := [* CharacterTable(H) : H in ES *];

    Y := [* *];
    s := 0;
    for i in [1..#ES] do
        for psi in CTs[i] do
           if Degree(psi) eq 1 then
               Append(~Y, [* ES[i], psi *] );
               s := s+1;
           end if;
        end for;
    end for;
    
    TG := CharacterTable(G);
    A := ZeroMatrix(IntegerRing(), #TG, s);
    for i := 1 to #TG do
        j := 1;
        for y in Y do
           A[i, j] := InnerProduct(TG[i], Induction(y[2], G));
           j := j+1;
        end for;
    end for;
    b := Vector( [Integers() ! InnerProduct(chi, psi) : psi in TG] );
    c := Solution(Transpose(A), b);
    /* return [ [* c[i], Y[i] *] : i in [1..Ncols(c)] | c[i] ne 0 ]; */
    return [ [* Y[i][1], Y[i][2], c[i] *] : i in [1..Ncols(c)] | c[i] ne 0 ];

end intrinsic;


/* ****************************** Det ******************************
 *
 * Det(U :: GrpPerm, phi :: AlgChtrElt, lambda :: AlgGrpElt) -> AlgMatElt: 
 * Let chi := ind_U^G(phi). phi is assumed to be linear. 
 * Then we compute Det_chi(lambda), i.e. we
 * compute a representation T associated with chi and then compute
 * Determinant(T(lambda)).
 *
****************************************************************************** */

intrinsic Det(U :: GrpPerm, phi :: AlgChtrElt, lambda :: AlgGrpElt) -> AlgMatElt
{}

    local G, E, Gl, Tphi, T, n, M, g;

    assert Degree(phi) eq 1;
    G := Group(Parent(lambda));
    E := Parent(phi(Id(U)));
    Gl := GL(1,E);

    Tphi := hom< U->Gl | [gen->Gl![phi(gen)] : gen in Generators(U)] >;
    T := Induction(Tphi, G);

    n := Index(G, U);

    M := ZeroMatrix(E,n,n);
    for g in G do
        M := M + ScalarMatrix(E, n, Coefficient(lambda, g)) * T(g);
    end for;

    return Determinant(M);
end intrinsic;



/* ****************************** Det ******************************
 *
 * Det(U :: GrpPerm, phi :: AlgChtrElt, S :: Mtrx) -> AlgMatElt:
 * Let chi := ind_U^G(phi). phi is assumed to be linear. 
 * S is a square matrix with entries in QG.
 * We compute Det_chi(S), i.e. we
 * compute a representation T associated with chi and then compute
 * Determinant( (T(lambda_{i,j}))_{i,j} ).
 *
****************************************************************************** */

intrinsic Det(U :: GrpPerm, phi :: AlgChtrElt, S :: Mtrx) -> AlgMatElt
{}

    local G, E, Gl, Tphi, T, n, M, g, BigM, i, j, Row;

    assert Degree(phi) eq 1;
    G := Group(Parent(S[1,1]));
    E := Parent(phi(Id(U)));
    Gl := GL(1,E);

    Tphi := hom< U->Gl | [gen->Gl![phi(gen)] : gen in Generators(U)] >;
    T := Induction(Tphi, G);

    n := Index(G, U);
    m := Nrows(S);
    
    BigM := ZeroMatrix(E, 0, n*m);
    for i:=1 to m do
        Row := ZeroMatrix(E, n, 0);
        for j := 1 to m do
            M := ZeroMatrix(E,n,n);
            for g in G do
                M := M + ScalarMatrix(E, n, Coefficient(S[i,j], g)) * T(g);
            end for;
	    Row := HorizontalJoin(Row, M);
	end for;
	BigM := VerticalJoin(BigM, Row);
    end for;

    return Determinant(BigM);
end intrinsic;


/* ****************************** InitBrauer ******************************
 *
 * InitBrauer(H :: List) -> List: Computes a BrauerStructure for each component
 * of QG.
 *
****************************************************************************** */

intrinsic InitBrauer(H :: List) -> List
{}

    local B, chi, b;

    B := [* *];
    for h in H do
       b := rec<BrauerStructure | >;
       b`chi := h[2];
       b`Uphis := BrauerInduction(h[2]);
       b`iota := h[3];
       Append(~B, b);
    end for;

    return B;
end intrinsic;


/* ****************************** ReducedNorm ******************************
 *
 * ReducedNorm(BrauerList :: List, lambda :: AlgGrpElt) -> List: 
 * ReducedNorm(lambda :: AlgGrpElt) -> List
 * Computes the reduced norm of lambda as an element in Z(QG) \simeq \prod K_\chi.
 * In each component we use the formula
 *
 *        nr(lambda)_chi = \prod_{U,phi} Det_{\ind_U^G(phi)}^{c_{U,phi}} (lambda)
 *
****************************************************************************** */

intrinsic ReducedNorm(BrauerList :: List, lambda :: AlgGrpElt) -> List
{}

    local N;

    N := [* *];
    for b in BrauerList do
        norm := &*[ Det(t[1], t[2], lambda)^t[3] : t in b`Uphis ];
	Append(~N, b`iota(norm));
    end for;
    return N;
end intrinsic;

intrinsic ReducedNorm(lambda :: AlgGrpElt) -> List
{}

    local N;

    N := [* *];
    for b in Parent(lambda)`Brauer do
        norm := &*[ Det(t[1], t[2], lambda)^t[3] : t in b`Uphis ];
        Append(~N, b`iota(norm));
    end for;
    return N;
end intrinsic;


intrinsic NewtonReducedNorm(lambda :: AlgGrpElt) -> List
{}

    local N, h, chi, n, sigma, s, k;

    N := [* *];
    for h in Parent(lambda)`H do
        chi := h[2];
	n := Integers() ! Degree(chi);
	sigma := [LinExtChar(chi, lambda)];
	s := [ sigma[1] ];
	for k:=2 to n do
	    Append(~s, LinExtChar(chi, lambda^k));
	    Append(~sigma, (-1)^(k+1) * (s[k] + &+[ (-1)^l * s[k-l]*sigma[l] : l in [1..k-1] ]) / k);
	end for;
	Append(~N, h[3](sigma[n]) );
   end for;
   return N;
end intrinsic;

intrinsic NewtonReducedNorm(S :: Mtrx) -> List
{}

    local N, h, chi, n, sigma, s, k, M, m, j;

    N := [* *];
    m := Ncols(S);
    for h in Parent(S[1,1])`H do
        M := S;
        chi := h[2];
        n := Integers() ! Degree(chi);
        sigma := [ &+[ LinExtChar(chi, M[j,j]) : j in [1..m] ] ];
        s := [ sigma[1] ];
        for k:=2 to n*m do
            M := M*S;
            Append(~s, &+[LinExtChar(chi, M[j,j]) : j in [1..m] ]);
            Append(~sigma, (-1)^(k+1) * (s[k] + &+[ (-1)^l * s[k-l]*sigma[l] : l in [1..k-1] ]) / k);
        end for;
        Append(~N, h[3](sigma[n*m]) );
   end for;
   return N;
end intrinsic;



/* ****************************** ReducedNorm ******************************
 *
 * ReducedNorm(S :: Mtrx) -> List
 * Computes the reduced norm of the matrix S as an element in Z(QG) \simeq \prod K_\chi.
 * In each component we use the formula
 *
 *        nr(S)_chi = \prod_{U,phi} Det_{\ind_U^G(phi)}^{c_{U,phi}} (S)
 *
****************************************************************************** */

intrinsic ReducedNorm(S :: Mtrx) -> List
{}

    local N;

    N := [* *];
    for b in Parent(S[1,1])`Brauer do
        norm := &*[ Det(t[1], t[2], S)^t[3] : t in b`Uphis ];
        Append(~N, b`iota(norm));
    end for;
    return N;
end intrinsic;



/* ****************************** ResidueClassGroup ******************************************
 *
 * ResidueClassGroup(QG :: AlgGrp, Ideals :: List) -> GrpAb, List, List
 * Given a list of ideals (\frakg_1, ..., \frakg_r) in Ideals
 * we compute the residue class group ( OC/\frakg )^* and a list m.
 * m is a list of maps which
 * in each component is the inverse of the map returned 
 * by the MAGMA function MultiplicativeGroup.
 * The second returned list is Ideals (just to be compatible with the p-version of
 * this function).
 *
***************************************************************************************/

intrinsic ResidueClassGroup(QG :: AlgGrp, Ideals :: List) -> GrpAb, List, List
{}
    local i, ordA, m, QR, A, f;

    ordA := [];
    m := [* *];

    for i:=1 to #Ideals do
        QR := quo< QG`H[i][4] | Ideals[i] >;
        A, f := MultiplicativeGroup(QR);
        Append(~m, f^-1);
        if NumberOfGenerators(A) eq 0 then
            Append(~ordA, [ 1 ]);
        else
            Append(~ordA, [Order(A.j) : j in [1..NumberOfGenerators(A)]]);
        end if;
    end for;
    OCmodG := AbelianGroup( &cat ordA );

    return OCmodG, m, Ideals;
end intrinsic;


/* ****************************** ResidueClassGroup ******************************************
 *
 * ResidueClassGroup(QG :: AlgGrp, Ideals :: List, p :: RngIntElt) -> GrpAb, List, List
 * Given a list of ideals (\frakg_1, ..., \frakg_r) in Ideals
 * we compute the residue class group ( OC_p/\frakg_p )^* and two lists m and pIdeals.
 * m is a list of maps which
 * in each component it is the inverse of the map returned
 * by the MAGMA function MultiplicativeGroup.
 * pIdeals is the list of p-primary parts of the ideals in Ideals. 
 *
***************************************************************************************/

intrinsic ResidueClassGroup(QG :: AlgGrp, Ideals :: List, p ::RngIntElt) -> GrpAb, List, List
{}
    local i, ordA, m, QR, A, f, pIdeals;

    ordA := [];
    m := [* *];
    pIdeals := [* *];

    for i:=1 to #Ideals do
        Id := pPrimaryPart(Ideals[i], p);
        QR := quo< QG`H[i][4] | Id >;
        A, f := MultiplicativeGroup(QR);
        Append(~m, f^-1);
	Append(~pIdeals, Id);
        if NumberOfGenerators(A) eq 0 then
            Append(~ordA, [ 1 ]);
        else
            Append(~ordA, [Order(A.j) : j in [1..NumberOfGenerators(A)]]);
        end if;
    end for;
    OCmodG := AbelianGroup( &cat ordA );

    return OCmodG, m, pIdeals;
end intrinsic;


/* ****************************** pPrimaryPart ******************************************
 *
 * pPrimaryPart(Id :: RngOrdIdl, p ::RngIntElt) -> RngOrdIdl
 * Given an ideal Id and a rational prime p, this function computes the p-primary
 * part of Id.
 *
***************************************************************************************/

intrinsic pPrimaryPart(Id :: RngOrdFracIdl, p ::RngIntElt) -> RngOrdFracIdl
{}
    local f, Idp, x;
  
    f := Factorization(Id);
    Idp := 1*Order(Id);
    for x in f do
        if IsPowerOf(#quo< Order(Id) | x[1] >, p) then
	    Idp := Idp * x[1]^x[2];
	end if;
    end for;
    return Idp;
end intrinsic;


/* ****************************** TorsionSubgroup ******************************************
 *
 * TorsionSubgroup(QG :: AlgGrp, OCmodG :: GrpAb, m :: List, K1AmodF :: SeqEnum) -> GrpAb, Map:
 * Given the residue class group  OCmodG and a list of generators of K1AmodF
 * we compute the relations given by these generators. As quotient we obtain K_0(ZG, Q)_{tors}
 * or K_0(Z_p[G], Q_p), depending on OCmodG and K1AmodF.
 * We return the abelian group K_0(ZG, Q)_{tors} or K_0(Z_p[G], Q_p),
 * and the map f from OCmodG to this quotient.
 *
***************************************************************************************/

intrinsic TorsionSubgroup(QG :: AlgGrp, OCmodG :: GrpAb, m :: List, K1AmodF :: SeqEnum) -> GrpAb, Map
{}

    local R, count, alpha, row, Id, i, clA, D, v;

    R := [];
    count := 1;
    vprintf RelKTh, 1: "We will compute %o relations\n", #K1AmodF;
    for alpha in K1AmodF do
        assert IsInvertible(alpha);
        if count mod 10 eq 0 then
            vprintf RelKTh, 1: "%o relations computed\n", count;
        end if;
        row := [];
        /* nr := ReducedNorm(QG`Brauer, alpha); */
        nr := NewtonReducedNorm(alpha);
        for i:=1 to #nr do
            v := ElementToSequence( (m[i]) (nr[i]) );
            if #v eq 0 then
                v := [0];
            end if;
            row cat:= v;
        end for;
        Append(~R, row);
        count +:= 1;
    end for;
    DTZG := rec< DTGrp | >;
    DTZG`DT, DTZG`f := quo< OCmodG | [OCmodG ! a : a in R] >;
    DTZG`OCmodG := OCmodG;
    DTZG`m := m;

    return DTZG;
end intrinsic;


/* ****************************** K0RelTorsion ******************************************
 *
 * K0RelTorsion(QG :: AlgGrp, ZG :: Rec, F :: Rec, Ideals :: List, p :: RngIntElt) -> Rec
 * Given ZG, a two-sided ideal F \sseq ZG of a maximal order,  the List Ideals = [ \frg_1,...,\frg_r ]
 * and a prime p, we compute K_0(Z_p[G], Q_p)_tors. The main work is done in TorsionSubgroup.
 *
***************************************************************************************/

intrinsic K0RelTorsion(QG :: AlgGrp, ZG :: Rec, F :: Rec, Ideals :: List, p :: RngIntElt) -> Rec
{}

    local K1AmodF, OCmodG, DT;

    K1AmodF, exact, ps, qs, Qs := Oumf(QG, ZG, F, Ideals,p);
    OCmodG, m, pIdeals := ResidueClassGroup(QG, Ideals,p);
    DT:= TorsionSubgroup(QG, OCmodG, m, K1AmodF);
    DT`p := p;
    DT`Gp := pIdeals;
    return DT;
end intrinsic;


/* ****************************** K0Rel******************************************
 *
 * K0Rel(QG :: AlgGrp, ZG :: Rec, F :: Rec, Ideals :: List, p :: RngIntElt) -> Rec
 * Given ZG, a two-sided ideal F \sseq ZG of a maximal order,  the List Ideals = [ \frg_1,...,\frg_r ]
 * and a prime p, we compute K_0(Z_p[G], Q_p). The main work is done in TorsionSubgroup.
 *
***************************************************************************************/

intrinsic K0Rel(QG :: AlgGrp, ZG :: Rec, F :: Rec, Ideals :: List, p :: RngIntElt) -> Rec
{}

    local K1AmodF, OCmodG, DT, K;

    K1AmodF := Oumf(QG, ZG, F, Ideals,p);
    OCmodG, m, pIdeals := ResidueClassGroup(QG, Ideals,p);
    DT:= TorsionSubgroup(QG, OCmodG, m, K1AmodF);

    K := rec< K0RelGrp | >;
    K`DT := DT`DT;
    K`OCmodG := DT`OCmodG;
    K`m := DT`m;
    K`f := DT`f;
    K`p := p;
    K`Gp := pIdeals;
    K`PrimeIdeals := [* [ v[1] : v in Factorization(p*h[4])]  : h in QG`H*];
    K`PrimeElts := [* Uniformizer(K`PrimeIdeals[i], K`Gp[i])   : i in [1..#K`Gp] *];

    return K;
end intrinsic;


/* ****************************** Uniformizer ******************************************
 *
 * Uniformizer(primes :: SeqEnum) -> SeqEnum: For a SeqEnum of prime ideals (P_1,...P_s),
 * this function returns a SeqEnum of uniformizing elements (pi_1, ..., pi_s) such that
 *          pi_i \equiv 1 (\mod I_{P_l})
 * for all l ne i.
 *
***************************************************************************************/

intrinsic Uniformizer(primes :: SeqEnum, I :: RngOrdIdl) -> SeqEnum
{}

    local OK, X, i, U, M, temp;

    OK := Order(primes[1]);
    X := [ OK ! 1 : p in primes ];
    U := [ OK ! 0 : p in primes ];
    M := [ p^Max(1, Valuation(I, p)) : p in primes ];

    for i:=1 to #primes do
        temp := M[i];
	M[i] := primes[i]^2;
        X[i] := UniformizingElement(primes[i]);
        U[i] := ChineseRemainderTheorem(X, M);
        X[i] := OK ! 1;
	M[i] := temp;
    end for;

    return U;
end intrinsic;



/* ****************************** Uniformizer ******************************************
 *
 * Uniformizer(primes :: SeqEnum) -> SeqEnum: For a SeqEnum of prime ideals (P_1,...P_s), 
 * this function returns a SeqEnum of uniformizing elements (pi_1, ..., pi_s) such that
 * val_{P_i}(pi_j) = \delta_{ij}.
 *
***************************************************************************************/

intrinsic Uniformizer(primes :: SeqEnum) -> SeqEnum
{}
  
    local OK, X, i, U;

    OK := Order(primes[1]);
    X := [ OK ! 1 : p in primes ];
    U := [ OK ! 0 : p in primes ];
    M := [ p : p in primes ];
    
    for i:=1 to #primes do
        M[i] := primes[i]^2;
        X[i] := UniformizingElement(primes[i]);
	U[i] := ChineseRemainderTheorem(X, M);
	X[i] := OK ! 1;
        M[i] := primes[i];
    end for;

    return U;
end intrinsic;


/* ****************************** K0RelLog ******************************************
 *
 * K0RelLog(K :: Rec, TOmega :: Rec) -> List, GrpAbElt
 * Here K is a rec K0Grp and TOmega a rec K0RelElt.
 * This function computes the discrete log of TOmega.
 * K_0(\Zp[G], \Qp[G]) is represented as I_C(p) \times DT(\Zp[G]).
 * The List contains the ideal part, the second return value is an element in K`DT.
 *
 * Changes: 24.2.17 Inserted Transpose(S) to make it compatible with Ruben's RelAlgKTheory.m
 *          Also, without that any computation of Ruben returns false, which contradicts the theory.
 *
***************************************************************************************/

intrinsic K0RelLog(K :: Rec, TOmega :: Rec) -> List
{}

    local QG, BasisA, BasisB, S, nr, IdealPart, i, j, e, row, v;

    QG := K`QG;

    BasisA := LocalBasis(TOmega`A, K`p);
    BasisB := LocalBasis(TOmega`B, K`p);

    S :=  QGMatrix(QG, TOmega`Theta, TOmega`A, BasisA, TOmega`B, BasisB);
    S := Transpose(S);   /* Changed on 24.2.17 */
    /* nr:= ReducedNorm(S); */
    nr:= NewtonReducedNorm(S);
    
    IdealPart := [* 1*h[4] : h in QG`H *];
    for i:=1 to #nr do
        for j:=1 to #K`PrimeIdeals[i] do
            e := Valuation(nr[i], K`PrimeIdeals[i][j]);
	    IdealPart[i] := IdealPart[i]*K`PrimeIdeals[i][j]^e;
	    nr[i] := nr[i]*K`PrimeElts[i][j]^(-e);
	end for;
    end for;

    row := [];
    for i:=1 to #nr do
        v := ElementToSequence( (K`m[i])(Numerator(nr[i])) - ((K`m[i])(Denominator(nr[i]))) );
        if #v eq 0 then
            v := [0];
        end if;
        row cat:= v;
    end for;

    return [* IdealPart, K`f(K`OCmodG!row) *];

end intrinsic;


/* ****************************** K0RelLog ******************************************
 *
 * K0RelLog(K :: Rec, TOmega :: Rec) -> List, GrpAbElt
 * Here K is a rec K0Grp and TOmega a rec K0RelElt.
 * This function computes the discrete log of TOmega.
 * K_0(\Zp[G], \Qp[G]) is represented as I_C(p) \times DT(\Zp[G]).
 * The List contains the ideal part, the second return value is an element in K`DT.
 *
***************************************************************************************/

intrinsic K0RelLog(K :: Rec, nr :: List) -> List
{}

    local QG, IdealPart, i, j, e, row, v;

    QG := K`QG;

    IdealPart := [* 1*h[4] : h in QG`H *];
    for i:=1 to #nr do
        for j:=1 to #K`PrimeIdeals[i] do
            e := Valuation(nr[i], K`PrimeIdeals[i][j]);
	    IdealPart[i] := IdealPart[i]*K`PrimeIdeals[i][j]^e;
	    nr[i] := nr[i]*K`PrimeElts[i][j]^(-e);
	end for;
    end for;

    /* Achtung: Diese Funktion ist fehlerhaft!!! Für jedes i ist nr[i]*K`PrimeElts[i][j]^(-e)
       zu berechnen und anschließend mit dem Chin.Restsatz ein lambda_i zu bestimmen, so dass

       		lambda_i \equiv nr[i]*K`PrimeElts[i][j]^(-e) (\mod K`PrimeIdeals[i][j]^v(\frg_p,i)))

       Alternativ könnte man die Primelemente so bestimmen, dass sie kongruent 1 sind
       mod \frg_i / \frP^v_\frP(\frg_i).
    */
	
    row := [];
    for i:=1 to #nr do
        v := ElementToSequence( (K`m[i])(Numerator(nr[i])) - ((K`m[i])(Denominator(nr[i]))) );
        if #v eq 0 then
            v := [0];
        end if;
        row cat:= v;
    end for;

    return [* IdealPart, K`f(K`OCmodG!row) *];

end intrinsic;




/* ****************************** ClassGroupLog ******************************************
 *
 * ClassGroupLog(cl :: Rec, P :: Rec) -> GrpAbElt
 * Here cl is a rec LocFreeClassGrp  and P a rec ZGModule.
 * This function computes the discrete log of P.
 * It is assumed that P is locally free; this is not checked.
 *
***************************************************************************************/

intrinsic ClassGroupLog(cl :: Rec, P :: Rec) -> GrpAbElt
{}

    local QG, rho, A, index, ps, f, TOmega, p, K, omega, t, off, gamma, i, g, Oki, CoD, D, w, eta, y, IdB; 

    QG := cl`QG;
    G := Group(QG);

    /* Compute [P, id, ZG^rk(P)] as an element of the relative group. */
    rankP := Integers() ! ( Nrows(P`hnf) / #G );
    rho := RegularRep(QG);
    A := RegularModule(rho);
    F := A;
    for i:=2 to rankP do
        F := ZGModuleDirectSum(F, A);
    end for;
    TOmega := rec< K0RelElt | >;
    TOmega`A := P;
    TOmega`B := F;
    TOmega`Theta := IdentityMatrix(Rationals(), rankP * #G);

    /* Compute primes p for which P_p ne Z_pG^rankP. */
    index := ZGModuleIndex(F, P);
    ps := [ f[1] : f in Factorization( Numerator(index) )];
    ps cat:= [ f[1] : f in Factorization( Denominator(index) )];

    print "From ClassGroupLog: ps = ", ps;
    
    log := cl`cl ! 0;
    for p in ps do
        print "p = ", p;
        K := RelativeGroup(cl`QG, p);
	print "Compute the log in K0Rel";
        omega := K0RelLog(K, TOmega);
	print "log computed";
	t := ElementToSequence( ((K`f)^-1)(omega[2]) );
        off := 0;
	IdB := [* *];
	for i:=1 to #QG`H do 
            // print "i = ", i;
            /* Compute gamma_i as in Bley/Wilson, Section 5. */
	    g := K`m[i];
            OKi := Order(cl`Ideals[i]);
	    CoD := Codomain(g);
	    D := Domain(g);
	    w := [ t[off + i] : i in [1..Ngens(CoD)] ];
	    off := off + Ngens(CoD);

	    eta := OKi ! (g^-1)(CoD!w);
	    if eta eq (OKi ! 0) then
	        eta := (OKi ! 1);
		assert #CoD eq 1;
	    end if;
	    y := &*[ K`PrimeElts[i][j]^Valuation(omega[1][i], K`PrimeIdeals[i][j]) : j in [1..#K`PrimeIdeals[i]] ];
	    gamma := y * eta;

            /* Set up data to do the weak approximation step. */
            fac := Factorization(cl`Ideals[i]);
            alpha := [* OKi ! 1 : f in fac *];
            for j:=1 to #fac do
	        if pPrimaryPart(fac[j][1], p) eq fac[j][1] then
		    alpha[j] := gamma;
		end if;
	    end for;
	    beta := Approximation(alpha, fac, cl`Ideals[i],  cl`InfinitePlaces[i]);

            /* Finally compute the ideal whose class in the ray classgroup represents [P_p,id, Z_pG^rk(P)]. */
	    Append(~IdB, pPrimaryPart(gamma*OKi, p) * beta);
	end for;

        /* Now read IdB (which comes componentwise) as an element of the locally free class group. */
        cl_elt := [];
        for i:=1 to #IdB do
            NumB, DenB := IdealNumDen(IdB[i]);
            v := ElementToSequence( (cl`m[i])(NumB) - (cl`m[i])(DenB) );
            if #v eq 0 then
                v := [0];
            end if;
            cl_elt cat:= v;
        end for;

        /* Add the elements of cl for the different p. */
        print "Contribution for p = ", p, " ----> ", cl`f (cl`rcgp ! cl_elt);
        log +:= cl`f (cl`rcgp ! cl_elt);
    end for;
    return log;
end intrinsic;


/* ****************************** Approximation ******************************************
 *
 * Approximation(alpha :: SeqEnum, fac :: SeqEnum, IdG :: RngOrdIdl, ArchPlaces :: SeqEnum) -> RngCycElt
 * alpha = [alpha_\frp : \frp \mid IdG] is a sequence of elements of K, K = Quot(OK), OK = Order(IdG), 
 * fac is the factorization of IdG as obtained by Factorization(IdG),
 * ArchPlaces ist the sequence of all real archemedian places of K.
 * This function computes a totally positive beta in K, such that
 *		v_\frp(\alpha_\frp \beta - 1) gt v_\frp(IdG), \forall \frp \mid IdG.
 *
***************************************************************************************/

intrinsic Approximation(alpha :: List, fac :: SeqEnum, IdG :: RngOrdIdl, ArchPlaces :: SeqEnum) -> RngCycElt
{}

    local OK, uniformizer, pi, c, f, j, Q, s1, Qmal, s2, num, den, y, beta;

    OK := Order(IdG);
    uniformizer := Uniformizer( [f[1] : f in fac] );
    pi := &*[ uniformizer[j]^( -Valuation(alpha[j], fac[j][1]) ) : j in [1..#fac] ] / One(OK);
    c := [];
    for j:=1 to #fac do
       Q, s1 := quo< OK | fac[j][1]^fac[j][2] >;
       Qmal, s2 := MultiplicativeGroup(Q);
       num, den := numden(alpha[j]*pi, fac[j][1]);
       Append(~c, OK ! s2( (s2^(-1))(s1(den)) - (s2^(-1))(s1(num)) ) ); 
    end for;

    y := ChineseRemainderTheorem(c, [f[1]^f[2] : f in fac]);
    beta := y * pi;
    if #ArchPlaces gt 0 then
        beta := MakeTotallyPositive(beta, IdG);
    end if;
    return beta;
end intrinsic;

/* The same as above, but alpha is of type SeqEnum instead of List. */
intrinsic Approximation(alpha :: SeqEnum, fac :: SeqEnum, IdG :: RngOrdIdl, ArchPlaces :: SeqEnum) -> RngCycElt
{}

    local OK, uniformizer, pi, c, f, j, Q, s1, Qmal, s2, num, den, y, beta;

    OK := Order(IdG);
    uniformizer := Uniformizer( [f[1] : f in fac] );
    pi := &*[ uniformizer[j]^( -Valuation(alpha[j], fac[j][1]) ) : j in [1..#fac] ] / One(OK);
    c := [];
    for j:=1 to #fac do
       Q, s1 := quo< OK | fac[j][1]^fac[j][2] >;
       Qmal, s2 := MultiplicativeGroup(Q);
       num, den := numden(alpha[j]*pi, fac[j][1]);
       Append(~c, OK ! s2( (s2^(-1))(s1(den)) - (s2^(-1))(s1(num)) ) );
    end for;

    y := ChineseRemainderTheorem(c, [f[1]^f[2] : f in fac]);
    beta := y * pi;
    if #ArchPlaces gt 0 then
        beta := MakeTotallyPositive(beta, IdG);
    end if;
    return beta;
end intrinsic;



/* ****************************** MakeTotallyPositive ******************************************
 *
 * intrinsic MakeTotallyPositive(a :: FldOrdElt, IdG :: RngOrdIdl) -> FldOrdElt
 * Computes a totally positive  element b in K, K = Parent(a), such that v_\frp(a- b) ge v_\frp(IdG).
 * This is essentially [Cohen 2, Alg. 4.2.20].
 *
***************************************************************************************/

intrinsic MakeTotallyPositive(a :: FldOrdElt, IdG :: RngOrdIdl) -> FldOrdElt 
{}
    local F2, r, A, rk, y, count, versuch, b, v, B, rk1, w;

    F2 := GF(2);
    r := Signature( Parent(a) );
    w := RealSign(a);
    if w eq Vector([F2!0 : i in [1..r]]) then
        return a;  /* a is already totally positive. */
    end if;
    A := Matrix(F2, r, 0, []);
    bound := 5;
    rk := 0;
    y := [];
    count := 0; versuch := 1;
    while rk lt r do
        b := Parent(a) ! Random(IdG, bound);
	count +:= 1;
	v := RealSign(1+b);
	B := HorizontalJoin(A, Matrix(r,1,ElementToSequence(v)));
	rk1 := Rank(B);
	if rk lt rk1 then
	    A := B;
	    rk := rk1;
	    Append(~y, 1+b);
	end if;
	if count eq 2^(versuch+1) * r then
	    versuch +:= 1;
	    bound +:= 5;
	end if;
    end while;
    
    w := (A^(-1)) * Matrix(r,1, ElementToSequence(RealSign(a)));
    b := a;
    for i:=1 to r do
        if w[i][1] ne F2!0 then
	    b *:= y[i];
	end if;
    end for;
    return b;
end intrinsic;


/* ****************************** RealSign ******************************************
 *
 * intrinsic RealSign(a :: FldOrdElt) -> ModTupFldElt
 * Computes the vector v in F_2^r, r = number of real places, such that sign( sigma_i(a) ) = (-1)^v[i].
 *
***************************************************************************************/

intrinsic RealSign(a :: FldOrdElt) -> ModTupFldElt
{}
    local F2, r, c, v, i;

    F2 := GF(2);
    r := Signature( Parent(a) );
    c := Conjugates(a);
    R := RealField( Parent(c[1]) );
    v := Vector( r, [F2!0 : i in [1..r]] );
    for i:=1 to r do
        if R!c[i] lt 0 then
	    v[i] := F2!1;
	end if;
    end for;
    return v;
end intrinsic;
    
/* ****************************** IdealNumDen ******************************************
 *
 * IdealNumDen(I :: RngOrdFracIdl) -> RngOrdIdl, RngOrdIdl
 * Computes the numerator and denominatorof the fractional ideal I. Numerator and 
 * denominator are relatively prime.
 *
***************************************************************************************/

intrinsic IdealNumDen(I :: RngOrdFracIdl) -> RngOrdIdl, RngOrdIdl
{}
    local Num, Den, F, f;

    Num := 1*Order(I);
    Den := 1*Order(I);
    F := Factorization(I);
    for f in F do
        if f[2] lt 0 then
            Den *:= f[1]^(-f[2]);
        end if;
        if f[2] gt 0 then
            Num *:= f[1]^f[2];
        end if;
    end for;
    return Num, Den;
end intrinsic;


/* ****************************** numden ******************************************
 *
 * numden(xi :: FldOrdElt, P :: RngOrdIdl) -> RngOrdElt, RngOrdElt
 * Computes a  numerator beta and denominator gamma of the field element xi = beta / gamma,
 * such that beta, gamma are integral and relatively prime to the given prime ideal P. 
 * It is assumed that the valuation of xi at P is 0.
 *
***************************************************************************************/

intrinsic numden(xi :: FldOrdElt, P :: RngOrdIdl) -> RngOrdElt, RngOrdElt
{}

    local OK, F, IdA, IdB, f, b, beta, gamma, p, j;

    assert Valuation(xi, P) eq 0;
    OK := Order(P);

    /* Compute the ideal numerator and denominator of xi*OK. Could be replaced by IdealNumDen. */
    F := Factorization(xi*OK);
    IdA := 1*OK;
    IdB := 1*OK;
    for f in F do
        if f[2] lt 0 then
	    IdB := IdB * f[1]^(-f[2]);
	end if;
	if f[2] gt 0 then
            IdA := IdA * f[1]^f[2];
        end if;
    end for;

    b, beta := IsPrincipal(IdA);
    if b then
        /* Ideal numerator and denominator are already principal. */
        gamma := beta / xi;
	return beta, gamma;
    end if;

    /* Find a prime Q such that IdA*Q is principal. This could be improved a lot. This was
       probably never tested, since in most applications in our context the rings of integers
       have class number 1. 
    */
    p := 2;
    found := false;
    while not found do
        F := Factorization(p*OK);
	j := 1;
	while not found and j le #F do
	    if F[j][1] ne P then
	        b, beta := IsPrincipal(F[j][1]*IdA);
		if b then
		    found := true;
		    gamma := beta / xi;
		else
		    j := j+1;
		end if;
	     end if;
	end while;
	p :=NextPrime(p);
     end while;
     return beta, gamma;
end intrinsic;



/* ****************************** CentralSimpleAlgRep ******************************************
 *
 * CentralSimpleAlgRep(A :: AlgAss) -> List:
 * Given a central simple algebra (as a AlgAss over a finite field) we compute the centre
 * F of A and an isomorphism A -> M = M_n(F). List contains all information for the functions
 * beta: M -> A and betainv : A -> M.
 *
***************************************************************************************/

intrinsic CentralSimpleAlgRep(A :: AlgAss) -> List
{}

    local Fp, F, FBasisA, L, FBasisL, m, B, b, om, T, j, y;

    Fp := BaseField(A);
    F := Centre(A);
    FBasisA := FindFBasis(A,F);
    L := MinimalLeftIdeals(A : Limit:=1)[1];
    FBasisL := FindFBasis(L, F);

    m := #FBasisL;    s := #Basis(F);
    B := [];
    for b in FBasisL do
        for om in Basis(F) do
            Append(~B, ElementToSequence( (A!om)*b) );
        end for;
    end for;
    B :=  Matrix(B) ;
    b :=Vector([Fp!0 : i in [1..m*s]]);
    T := ZeroMatrix(F, m, m*s);


    for j:=1 to m*s do
        b[j] := Fp ! 1;
        y := Solution(B, b);
	for k:=1 to m do
	    T[k, j] := F ! [y[t] : t in [(k-1)*s+1..k*s]];
	end for;
	b[j] := Fp ! 0;
    end for;

    D := [];
    for a in FBasisA do
       Tj := [];
       for b in FBasisL do
           Append(~Tj, a*b);
       end for;
       Tj := Transpose( Matrix(Tj));
       S := ZeroMatrix(F, Nrows(Tj), Ncols(Tj));
       for i:=1 to Nrows(Tj) do
           for j:=1 to Ncols(Tj) do
               S[i,j] := F!Tj[i,j];
           end for;
       end for;
       Append(~D, T*S);
    end for;

    Phi := [];
    for d in D do
        x :=[];
        for i:=1 to Nrows(d) do
            for j:=1 to Ncols(d) do
                Append(~x, d[i,j]);
            end for;
        end for;
        Append(~Phi, x);
    end for;
    Phi := Transpose( Matrix(Phi) );

    /* The following if  is necessary because there is a problem 
       with inversion of 1x1 matrices in MAGMA V2-17-9 
       Changed 28.9.2011   
    */
    PhiInv := MyInvertPhi(Phi);

    M := MatrixAlgebra(F, m);

    m := #FBasisA;    s := #Basis(F);
    B := [];
    for b in FBasisA do
        for om in Basis(F) do
            Append(~B, ElementToSequence( (A!om)*b) );
        end for;
    end for;
    B :=  Matrix(B) ;
    b :=Vector([Fp!0 : i in [1..m*s]]);
    T := ZeroMatrix(F, m, m*s);


    for j:=1 to m*s do
        b[j] := Fp ! 1;
        y := Solution(B, b);
        for k:=1 to m do
            T[k, j] := F ! [y[t] : t in [(k-1)*s+1..k*s]];
        end for;
        b[j] := Fp ! 0;
    end for;


    return M, [* FBasisA, Basis(M), Phi, PhiInv, T *];
end intrinsic; 

intrinsic MyInvertPhi(Phi :: AlgMatElt) -> AlgMatElt
{}
    local i;

    // print "In MyInvertPhi";

    F := Parent(Phi[1,1]);
    Fbasis := Basis(F);
    Fp := BaseField(F);
    p := Characteristic(Fp);
    d := Dimension(F);
    Fq := GF(p^d);
    alpha := Fq.1;
    mipo := MinimalPolynomial(alpha);
    C := [c : c in CartesianProduct( [[z : z in Fp] : i in [1..d]] )];
    beta_found := false; i := 1;
    while not beta_found do
        c := C[i];
	beta := &+[ c[j]*Fbasis[j] : j in [1..d] ];
	if MinimalPolynomial(beta) eq mipo then
	    beta_found := true;
	end if;
	i := i+1;
    end while;

    graph_iota := []; graph_iota_inv := []; Felts := [];
    for e in C do
        x := &+[e[i]*alpha^(i-1) : i in [1..#e]];
        y := &+[e[i]*beta^(i-1) : i in [1..#e]];
        Append(~graph_iota, <x,y>);
        Append(~graph_iota_inv, <y,x>);
        Append(~Felts, y);
    end for;
    iota := map<Fq->Felts | graph_iota >;
    iota_inv := map<Felts->Fq | graph_iota_inv >;
    Y := Matrix(Fq, Nrows(Phi), Ncols(Phi), [iota_inv(Phi[i,j]) : i,j in [1..Nrows(Phi)]]);
    Z := Y^-1; 
    PhiInv := Matrix(F,Nrows(Phi), Ncols(Phi), [iota(Z[i,j]) : i,j in [1..Nrows(Phi)]]);
    // print "Leave my InvertPhi";
    return(PhiInv);
end intrinsic;

        


/* ****************************** beta ******************************************
 *
 * beta(m :: AlgMatElt, Q :: List) -> AlgAssElt:
 * beta : M =M_n(F) -> A is a isomorphism of M to the central simple F-algebra A.
 * Q is a list computed by CentralSimpleAlgebraRep.
 *
***************************************************************************************/

intrinsic beta(m :: AlgMatElt, Q :: List) -> AlgAssElt
{}

    local PhiInv, FBasisA, A, z, a;

    PhiInv := Q[4];
    FBasisA := Q[1];
    A := Parent(FBasisA[1]);
    z := PhiInv * Transpose( Matrix([ ElementToSequence(m) ]) );
    a := &+[(A!z[i,1])*FBasisA[i] : i in [1..#FBasisA]];
    return a;
end intrinsic;


/* ****************************** betainv ******************************************
 *
 * betainv(a :: AlgAssElt, Q :: List) -> AlgMatElt:
 * The inverse of beta. Q is a list computed by CentralSimpleAlgebraRep.
 *
***************************************************************************************/

intrinsic betainv(a :: AlgAssElt, Q :: List) -> AlgMatElt
{}

    local Phi, F, T, FBasisM, M, x, z, m;

    Phi := Q[3];
    F := Parent(Phi[1,1]);
    T := Q[5];
    FBasisM := Q[2];
    M := Parent(FBasisM[1]);
    x := Transpose( Matrix([ [F!y : y in ElementToSequence(a)] ]) );
    z := Phi * T * x;
    m := &+[(M!z[i,1])*FBasisM[i] : i in [1..#FBasisM]];
    return m;
end intrinsic;



/* ****************************** LocalBasis ******************************************
 *
 * LocalBasis(D :: Rec, p :: RngIntElt) -> List:
 * For a p-locally free ZG-module D we compute a list [w_1, ..., w_m] such that
 * w_1, ..., w_m is a ZpG-basis of D_p.
 *
 * Changed on 18.9.2012: We now test whether D is locally free
 *
***************************************************************************************/
    
intrinsic LocalBasis(D :: Rec, p :: RngIntElt) -> List, BoolElt
{}

    local G, Dp, F, FG, J, R, psi, JDp, DpmodJ, idem, comp, i, A, e, gens, B, M, Q, v, w;

    G := Domain(D`phi);
    d := Nrows(D`hnf) / #G;
    if not IsIntegral(d) then
        return [**], false;
    else
        d := Integers() ! d;
    end if;

    Dp := ZGtoFGModule(D, p);
    F := GF(p);
    FG := GroupAlgebra(F, G);
    J := JacobsonRadical(FG);
    R, psi := quo<FG | J>;
    // print "Compute JDp";
    JDp := FGSubModule(Dp, J);
    // print "Compute DpmodJ";
    DpmodJ, tau := FGQuotientModule(Dp, JDp);

    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]);
        gens := [ GroupRingAction(DpmodJ, e, b) : b in Basis(DpmodJ`M) ];
        // print "#gens = ", #gens;
        Append(~C, FGSubModule(DpmodJ, gens));
    end for;

    B := [* *];
    for i:=1 to #comp do
        /* Berechne die Basis in der i-ten Komponente */
        M, Q := CentralSimpleAlgRep(comp[i]);
        v, isfree := MoritaBasis(M, Q, C[i], psi, R);
        if #v ne d or 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 [ LiftFGModuleElt(D, v) : v in w ], true;
end intrinsic;


/* ****************************** MoritaBasis ******************************************
 *
 * MoritaBasis(M :: AlgMat, Q :: List, P :: Rec, psi :: Map, R :: AlgAss) -> List
 * R is a semisimple Fp-algebra, M is a component of R, Q is a list computed by
 * CentralSimpleAlgebraRep which is needed to view M as a submodule of R, P is a free
 * M-module given as a Rec FGModule, psi is a map from R to FpG.
 *
 * We compute a M-basis of P.
 *
 * Changed on 18.9.2012: We also check whether P is free over M
 *
***************************************************************************************/

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

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

    n := Degree(M);
    E := Matrix(n, n, [(psi^-1)(R ! beta(b, Q)) : b in Basis(M)]);
    gens := [ GroupRingAction(P, E[1,1], b) : b in Basis(P`M) ];
    V := FSubModule(P, gens); 
    F := CoefficientRing(M);
    FBasisV := FindFBasis(V, M, Q, psi, R);
    if not IsDivisibleBy(#FBasisV, n) then
        return [**], false;
    end if;

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

    return w, true;
end intrinsic;



/* ****************************** LocalBasis ******************************************
 *
 * LocalBasis(D :: Rec, p :: RngIntElt) -> List:
 * For a p-locally free ZG-module D we compute a list [w_1, ..., w_m] such that
 * w_1, ..., w_m is a ZpG-basis of D_p.
 *
***************************************************************************************/
    
intrinsic OldLocalBasis(D :: Rec, p :: RngIntElt) -> List
{}

    local G, Dp, F, FG, J, R, psi, JDp, DpmodJ, idem, comp, i, A, e, gens, B, M, Q, v, w;

    G := Domain(D`phi);
    Dp := ZGtoFGModule(D, p);
    F := GF(p);
    FG := GroupAlgebra(F, G);
    J := JacobsonRadical(FG);
    R, psi := quo<FG | J>;
    // print "Compute JDp";
    JDp := FGSubModule(Dp, J);
    // print "Compute DpmodJ";
    DpmodJ, tau := FGQuotientModule(Dp, JDp);

    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]);
        gens := [ GroupRingAction(DpmodJ, e, b) : b in Basis(DpmodJ`M) ];
        // print "#gens = ", #gens;
        Append(~C, FGSubModule(DpmodJ, gens));
    end for;

    B := [* *];
    for i:=1 to #comp do
        /* Berechne die Basis in der i-ten Komponente */
        M, Q := CentralSimpleAlgRep(comp[i]);
        v := MoritaBasis(M, Q, C[i], psi, R);
        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 [ LiftFGModuleElt(D, v) : v in w ];
end intrinsic;


/* ****************************** MoritaBasis ******************************************
 *
 * MoritaBasis(M :: AlgMat, Q :: List, P :: Rec, psi :: Map, R :: AlgAss) -> List
 * R is a semisimple Fp-algebra, M is a component of R, Q is a list computed by
 * CentralSimpleAlgebraRep which is needed to view M as a submodule of R, P is a free
 * M-module given as a Rec FGModule, psi is a map from R to FpG.
 *
 * We compute a M-basis of P.
 *
***************************************************************************************/

intrinsic OldMoritaBasis(M :: AlgMat, Q :: List, P :: Rec, psi :: Map, R :: AlgAss) -> List
{}

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

    n := Degree(M);
    E := Matrix(n, n, [(psi^-1)(R ! beta(b, Q)) : b in Basis(M)]);
    gens := [ GroupRingAction(P, E[1,1], b) : b in Basis(P`M) ];
    V := FSubModule(P, gens); 
    F := CoefficientRing(M);
    FBasisV := FindFBasis(V, M, Q, psi, R);

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

    return w;
end intrinsic;

/* ****************************** TestLocalBasis ******************************************
 *
 * TestLocalBasis(QG :: AlgGrp, P :: Rec, basis :: SeqEnum, p :: RngIntElt) -> BoolElt
 * Checks whether the elements of basis form a ZpG-basis of the ZG-module P_p.
 *
***************************************************************************************/

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

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

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



/* ****************************** QGMatrix ******************************************
 *
 * QGMatrix(QG :: AlgGrp, Theta :: Mtrx, A :: Rec, BasisA :: SeqEnum, B :: Rec, BasisB :: SeqEnum) -> Mtrx
 * Theta defines a QG-isomorphism A_Q -> B_Q. We compute a matrix S in Gl_m(QG) such that
 * Theta(BasisA) = BasisB*S.
 *
***************************************************************************************/

intrinsic QGMatrix(QG :: AlgGrp, Theta :: Mtrx, A :: Rec, BasisA :: SeqEnum, B :: Rec, BasisB :: SeqEnum) -> Mtrx
{}

    local m, n, i, g, C, S, j, b, x, G;

    G := Group(QG);
    m := #BasisA; n := Dimension(Codomain(A`phi));
    C := [];
    for i:=1 to m do
        for g in G do
        Append(~C, ElementToSequence( GroupRingAction(B, QG!g, BasisB[i])) ) ;
        end for;
    end for;
    C := Transpose(Matrix(C));

    S := ZeroMatrix(QG, m, m);
    for j:=1 to m do
        b := BasisA[j]*Transpose(Theta);
	x := Solution(Transpose(C), b);
	for i:=1 to m do
	    S[i,j] := QG ! [x[l] : l in [(i-1)*#G+1..i*#G]];
	end for;
    end for;

    return S;
end intrinsic;
    

/************************ CommonSplittingField **************************************
 *
 * CommonSplittingField(QU :: AlgGrp,  QG :: AlgGrp) -> FldCyc:
 * This function returns the minimal cyclotomic field which contains all character values
 * of all absolutely irreducible characters of U and G.
 *
****************************************************************************************/
	
intrinsic CommonSplittingField(QU :: AlgGrp,  QG :: AlgGrp) -> FldCyc
{}

    local S, h, chi, E;

    S := [];
    for h in QG`H do
        chi := h[2];
        S cat:= [chi(g) : g in Group(QG)];
    end for;
    
    for h in QU`H do
        chi := h[2];
        S cat:= [chi(g) : g in Group(QU)];
    end for;
    E :=MinimalCyclotomicField(S);
    if Degree(E) eq 1 then
        return CyclotomicField(1);/* because MinimalCyclotomicField returns the rationals in this case */
    else
        return E;
    end if;
end intrinsic;


/************************ Induction **************************************
 *
 * Induction(KrelU :: Rec, KrelG :: Rec, t :: GrpAbElt) -> GrpAbElt
 * This function computes the induction map K_0(Z_p[U], Q_p)_tors --> K_0(Z_p[G], Q_p)_tors.
 *
****************************************************************************************/

intrinsic Induction(KrelU :: Rec, KrelG :: Rec, t :: GrpAbElt) -> GrpAbElt
{}

    local E,alpha, Tau, HG, HU, G, U, beta, i, chi, v, j, tau, row;

    E := CommonSplittingField(KrelU`QG, KrelG`QG);
    alpha := LiftToOC(KrelU,t);

    Tau := EmbeddingsToE(KrelU`QG, E);
    HG := KrelG`QG`H; 
    HU := KrelU`QG`H;
    G := Group(KrelG`QG);
    U := Group(KrelU`QG);
    beta := [];
    for i:=1 to #HG do
        chi := Restriction(HG[i][2], U);
	v := [];
	for j:=1 to #alpha do
	    v cat:= [tau(alpha[j])^ScalarProduct(U, tau, HU[j][2], HU[j][3], chi) : tau in Tau[j]];
	end for;
        Append(~beta, &*v);
    end for;
    row :=[];
    m := KrelG`m;
    Tau := EmbeddingsToE(KrelG`QG, E);

    for i:=1 to #beta do
	iota := Tau[i][1]; OKi := HG[i][4];
        /* For some groups there are PROBLEMS with coercion and domains of maps in the next line ??? */
        v := ElementToSequence( (m[i]) (OKi ! (iota^-1)(beta[i])) );
        if #v eq 0 then
            v := [0];
        end if;
        row cat:= v;
    end for;

    return KrelG`f( KrelG`OCmodG ! row );
end intrinsic;



/************************ Restriction **************************************
 *
 * Restriction(KrelG :: Rec, KrelU :: Rec, t :: GrpAbElt) -> GrpAbElt
 * This function computes the restriction map K_0(Z_p[G], Q_p)_tors --> K_0(Z_p[U], Q_p)_tors.
 *
****************************************************************************************/

intrinsic Restriction(KrelG :: Rec, KrelU :: Rec, t :: GrpAbElt) -> GrpAbElt
{}

    local E,alpha, Tau, HG, HU, G, U, beta, i, chi, v, j, tau, row;

    E := CommonSplittingField(KrelU`QG, KrelG`QG);
    alpha := LiftToOC(KrelG,t);

    Tau := EmbeddingsToE(KrelG`QG, E);
    HG := KrelG`QG`H;
    HU := KrelU`QG`H;
    G := Group(KrelG`QG);
    U := Group(KrelU`QG);
    beta := [];
    for j:=1 to #HU do
        chi := Induction(HU[j][2], G);
        v := [];
        for i:=1 to #alpha do
            v cat:= [tau(alpha[i])^ScalarProduct(G, tau, HG[i][2], HG[i][3], chi) : tau in Tau[i]];
        end for;
        Append(~beta, &*v);
    end for;
    row :=[];
    m := KrelU`m;
    Tau := EmbeddingsToE(KrelU`QG, E);

    for j:=1 to #beta do
        // iota := HG[j][3]; OKj := HG[j][4];
	// v := ElementToSequence( (m[j]) (OKj ! iota( beta[j] )) );          /* Problem */
	// print "j = ", j;
	iota := Tau[j][1]; OKj := HU[j][4];
       /* For some groups there are PROBLEMS with coercion and domains of maps in the next line ??? */
	v := ElementToSequence( (m[j]) (OKj ! (iota^-1)(beta[j])) );

        if #v eq 0 then
            v := [0];
        end if;
        row cat:= v;
    end for;

    return KrelU`f( KrelU`OCmodG ! row );
end intrinsic;



/************************ Quotient **************************************
 *
 * Quotient(KrelG :: Rec, KrelQ :: Rec, f :: Map, t :: GrpAbElt) -> GrpAbElt
 * This function computes the quotient map K_0(Z_p[G], Q_p)_tors --> K_0(Z_p[G/U], Q_p)_tors,
 * where U is a normal subgroup.
 * Here f is the natural map G -> G/U.
 *
****************************************************************************************/

intrinsic OldQuotient(KrelG :: Rec, KrelQ :: Rec, f :: Map, t :: GrpAbElt) -> GrpAbElt
{}

    local alpha, HG, HU, G, U, i, chi, v, j, row;

    alpha := LiftToOC(KrelG,t);

    HG := KrelG`QG`H;
    HQ := KrelQ`QG`H;
    G := Group(KrelG`QG);
    Q := Group(KrelQ`QG);
    X := KrelG`QG`X; /* The orbits of characters of G */

    row :=[];
    m := KrelQ`m;

    for j:=1 to #HQ do
        /* The next two lines were corrected on 25.9.2017 */
        chi := LiftCharacter(HQ[j][2], f, G);
        i := [i : i in [1..#X] | chi in X[i]][1];
        v := ElementToSequence( (m[j]) (alpha[i]) );

        if #v eq 0 then
            v := [0];
        end if;
        row cat:= v;
    end for;

    return KrelQ`f( KrelQ`OCmodG ! row );
end intrinsic;

intrinsic Quotient(KrelG :: Rec, KrelQ :: Rec, f :: Map, t :: GrpAbElt) -> GrpAbElt
{}

    local alpha, HG, HU, G, U, i, chi, v, j, row;

    alpha := LiftToOC(KrelG,t);

    HG := KrelG`QG`H;
    HQ := KrelQ`QG`H;
    G := Group(KrelG`QG);
    Q := Group(KrelQ`QG);
    X := KrelG`QG`X; /* The orbits of characters of G */

    row :=[];
    m := KrelQ`m;

    for j:=1 to #HQ do
        /* The next two lines were corrected on 25.9.2017 */
        chi := LiftCharacter(HQ[j][2], f, G);
	/* Find the orbit X[i] with chi in X[i] and compute iota such that iota(psi) = chi where psi = X[i][1]. */
        i := [i : i in [1..#X] | chi in X[i]][1];
	psi := X[i][1];
	G_Kpsi,_,q := AutomorphismGroup(CoefficientField(psi));
	iota := [q(g) : g in G_Kpsi | &and[q(g)(psi(h)) eq chi(h) : h in G]][1];

        // v := ElementToSequence( (m[j]) (alpha[i]) );
	v := ElementToSequence( (m[j]) ( iota(alpha[i]) ) );

        if #v eq 0 then
            v := [0];
        end if;
	print "j = ", j, "v = ", v;
        row cat:= v;
    end for;

    return KrelQ`f( KrelQ`OCmodG ! row );
end intrinsic;

/* Just an auxiliary function which computes
            < tau \circ psi, chi >_G 
*/

intrinsic ScalarProduct(U :: GrpPerm, tau :: Map, psi :: AlgChtrElt, iota :: Map, chi :: AlgChtrElt) -> RngIntElt
{}

    local u;

    return Integers() ! (&+[ tau( iota(psi(u)) )*chi(u^-1) : u in U] / #U);

end intrinsic;    

    
/*************************** LiftToOC *******************************************
 *
 * LiftToOC(KrelU :: Rec, t :: GrpAbElt) -> List:
 * t is an element in KrelU`DT = K_0(Z_p[G], Q_p)_tors. This group is a quotient of the
 * multiplicative group (O_C / \frg)^\times. We compute an element x \in O_C which maps to t.
 *
*******************************************************************************************/

intrinsic LiftToOC(KrelU :: Rec, t :: GrpAbElt) -> List
{}
    
    local s,m, pos, L, i, C, OKi, n, v;

    s := ElementToSequence( (KrelU`f^-1)(t) );
    m := KrelU`m;
    pos := 1;
    L := [* *];
    for i:=1 to #m do
        C := Codomain(m[i]);
	OKi := KrelU`QG`H[i][4];
	n := NumberOfGenerators(C);
	v := [s[j] : j in [pos..pos+n-1]];
	Append(~L, OKi ! (m[i]^-1)( C!v ) );
	pos := pos+n;
    end for;

    return L;
end intrinsic;
	

intrinsic EmbeddingsToE(QU :: AlgGrp, E :: FldCyc) -> List
{}

    local Tau, i, L, LX, mipo, roots;

    Tau := [* *];
    for i:=1 to #QU`H do
        L := QU`H[i][1];
        if Degree(L) eq 1 then
            theta := L!1;
	end if;
	if #Generators(L) eq 1 then
	    theta := L.1;
        else
            theta := CyclotomicPrimitiveElement(L);
	end if;
        mipo := MinimalPolynomial(theta);
        /* We construct K, which is "equal" to L,because we want a simple extension, given
           by a primitive element. This allows to define homomorphisms.
        */
        K := NumberField(mipo : DoLinearExtension:=true);
        KX := PolynomialRing(K);
        s := hom<K->L | theta>;
        roots := Roots(KX!mipo); /* This did not work for G = D_{38}, U = C_{19} and resulted
                                    in a segmentation fault with version 2.12-1.
                                    It works in version 2.13-3.
                                 */
        /* If necessary we interchange roots so that Tau[i][1] is the identity. */
        for k:=2 to #roots do
	    if roots[k][1] eq K.1 then
	        roots[k][1] := roots[1][1];
		roots[1][1] := K.1;
	    end if;
	end for;
        Append(~Tau, [ (s^-1)*hom<K->E | E ! s(r[1])> : r in roots ]);
    end for;

    return Tau;
end intrinsic;

/* Returns a primitive element for a cyclotomic field. This is usually "smaller" than
   the random element returned by the function PrimitiveElement.
*/

intrinsic CyclotomicPrimitiveElement(C :: FldCyc) -> FldCycElt
{}

    local xi, y;

    // print "In CyclotomicPrimitiveElement with C = ", C;
    xi := &*[y : y in Generators(C)];
    if IsPrimitive(xi) then
        return xi;
    else
        return PrimitiveElement(C);
    end if;
end intrinsic;





/* ****************************** InitGroupAlgebra ******************************
 *
 * InitGroupAlgebra(G::Grp, p :: RngIntElt) -> Rec: Computes the group K_0(Z_p[G], Q_p)
 * and enough additional data to solve the discrete logarithm problem.
 *
  ************************************************************************* */

intrinsic InitGroupAlgebra(G :: Grp) -> AlgGrp
{}

    local  QG, T, Omega, Zbasis;

    QG := GroupAlgebra(Rationals(), G); 

    e := Exponent(G); n := #G;
    F := CyclotomicField(e);
    IrrMod := IrreducibleModules(G, F);
    IrrRep := [ Representation(V) : V in IrrMod ];
    IrrChar := [ Character(d) : d in IrrRep ];
    Omega := CompOmega(G);      /* Omega = Galois group of Q(zeta_e) / Q */
    QG`Mods, QG`Reps, QG`X := OmegaOrbitsOfReps(IrrMod, Omega);
    QG`H := ComputeCharacterFields(QG);



    // T := CharacterTable(G);

    // Omega := CompOmega(G);      /* Omega = Galois group of Q(zeta_e) / Q */
    // QG`X := OmegaRepresentatives(T, Omega);
    // QG`H := ComputeKi(G, QG`X);

    QG`B := ComputeB(QG, QG`H);
    QG`Psi := ComputePsi(QG);
    Zbasis := QGCentre(QG, QG`X, QG`H);  /* NEVER CHANGE THIS LIST: it is used in the computation of PsiInv */
    QG`OC := Zbasis;
    // QG`Brauer := InitBrauer(QG`H);    /* Auskommentiert am 17.2.17 wegen eines MAGMA Fehlers */
   return QG;
end intrinsic;

			
/* ****************************** K0RelSub ******************************
 *
 * Subtraction of two elements in a relative group
 *
 * Implemented: 26.9.2012 by Werner Bley
 *
  ************************************************************************* */

intrinsic K0RelSub(a :: List, b :: List) -> List
{Subtracts two elements in a relative algebraic K-group}
    local c, i;

    c := [* [* a[1,i] / b[1,i] : i in [1..#a[1]] *], a[2] - b[2] *];
    return c;
end intrinsic;


/* ****************************** K0RelAdd ******************************
 *
 * Addition of two elements in a relative group
 *
 * Implemented: 26.9.2012 by Werner Bley
 *
  ************************************************************************* */

intrinsic K0RelAdd(a :: List, b :: List) -> List
{Adds two elements in a relative algebraic K-group}
    local c, i;

    c := [* [* a[1,i] * b[1,i] : i in [1..#a[1]] *], a[2] + b[2] *];
    return c;
end intrinsic;


/* ****************************** K0RelEquality ******************************
 *
 * Test for equality
 *
 * Implemented: 26.9.2012 by Werner Bley
 *
  ************************************************************************* */

intrinsic K0RelEquality(a :: List, b :: List) -> BoolElt
{Tests equality in a relative algebraic K-group}

    return &and[a[1,i] eq b[1,i] : i in [1..#a[1]]] and a[2] eq b[2];
end intrinsic;

/* ****************************** K0RelEqualityModTorsion ******************************
 *
 * Test for equality modulo torsion
 *
 * Implemented: 26.9.2012 by Werner Bley
 *
  ************************************************************************* */

intrinsic K0RelEqualityModTorsion(a :: List, b :: List) -> BoolElt
{Tests equality modulo torsion in a relative algebraic K-group} 

    return &and[a[1,i] eq b[1,i] : i in [1..#a[1]]];
end intrinsic;


/* ****************************** K0RelIsTrivialModTorsion ******************************
 *
 * Test for triviality mod torsion
 *
 * Implemented: 26.9.2012 by Werner Bley
 *
  ************************************************************************* */

intrinsic K0RelIsTrivialModTorsion(a :: List) -> BoolElt
{Tests equality modulo torsion in a relative algebraic K-group}

    return &and[ P eq 1*Order(P) : P in a[1] ];
end intrinsic;


/* ****************************** K0RelIsTrivial******************************
 *
 * Test for triviality
 *
 * Implemented: 26.9.2012 by Werner Bley
 *
  ************************************************************************* */

intrinsic K0RelIsTrivial(a :: List) -> BoolElt
{Tests equality in a relative algebraic K-group}

    return &and[ P eq 1*Order(P) : P in a[1] ] and a[2] eq Zero(Parent(a[2]));
end intrinsic;


