/*

This program contains Magma functions for calculating various homotopical
 functors and related constructions. For the program to run properly Version
 2.3 of Magma needs to have been installed, and the file homotopyhelp.m (which
 contains detailed explanations of the functions) needs to be present in the
 same directory as this program.

In writing the program the following conventions have been followed: an
 enumerated group is represented by a string (such as eG) whose first letter is
 a lower case "e", whereas a finitely presented group is represented by a
 string whose first letter is upper case; a homomorphism from a group G to a
 group H is denoted by GhomH.

I would like to thank John Cannon for his very generous help in simplifying,
 drastically shortening, and speeding up an earier version of this program.

 */ 



/* ******************************************************************** */
print "Type   'help();'   for information on the Homotopy Package.
";

forward ProductRelations, Remainder;

procedure help();
System("more homotopyhelp.m");
end procedure;
/* ********************************************************************* */



/* ********************************************************************* */    
EnumeratedGroup := function(G, p);

    if IsPrime(p) then
        eG, GhomeG := pQuotient(G, p, 63);
        eGhomG := hom< eG -> G | [ (eG.i)@@GhomeG :
                                   i in [1..NumberOfPCGenerators(eG)] ] >;
    else
        CSpc := CosetSpace(G, sub<G|>: CosetLimit := 2000000);
        GhomeG, eG := CosetAction(CSpc);
        eGhomG := hom< eG -> G | [ (eG.i)@@GhomeG : 
                                   i in [1..#Generators(G)] ] >;
    end if;
    return eG, GhomeG, eGhomG;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
EnumeratedKernel := function( delta );

    K := {};
    G := Domain(delta);
    trivialelt := delta(Id(G));
    for x in G do
        if x in K then continue; end if;
        if delta(x) eq trivialelt then 
            K := NormalClosure( G, sub<G | K, x>);
        end if;
    end for;
    return K;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
PresentedGroup := function( eG );

    G, GhomeG := FPGroup( eG);
    NGens := Type(eG) eq GrpPC select NumberOfPCGenerators(eG) else Ngens(eG);
    eGhomG := hom< eG -> G | [ (eG.i)@@GhomeG : i in [1..NGens] ] >;
    return G, GhomeG, eGhomG;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
SelfAction := function(G, p);

    eG, GhomeG, eGhomG := EnumeratedGroup(G, p);
    GactG := map< car<G,G>->G | x:->eGhomG(GhomeG(x[2])^GhomeG(x[1])) >;
    NiceGeneratorsG := &join {Generators(UpperCentralSeries(eG)[i]) : 
                              i in [2..#UpperCentralSeries(eG)] };
    T := Transversal(eG, sub< eG | NiceGeneratorsG > );
    return GactG, eGhomG(T) join eGhomG(NiceGeneratorsG);

end function;
/* ******************************************************************** */



/* ******************************************************************** */
Action := function(GhomeQ, HhomeQ, p);

    G := Domain(GhomeQ);
    H := Domain(HhomeQ);
    eQ := Codomain(GhomeQ);
    eG, GhomeG, eGhomG := EnumeratedGroup(G, p);
    // If H=G then set eH=eG etc., else enumerate H.
    if G cmpeq H then
        eH:=eG; HhomeH:=GhomeG; eHhomH:=eGhomG;
    else
        eH, HhomeH, eHhomH := EnumeratedGroup(H, p);
    end if;
    imG := GhomeQ(G);
    m := IsPrime(p) select #PCGenerators(imG) else #Generators(imG);
    S := [ eGhomG(x) : x in eG, i in [1..m] | GhomeQ(eGhomG(x)) eq imG.i ];
    imGhomG := hom< imG->G | S>;
    HactG := map< car<H,G>->G | x :-> imGhomG(GhomeQ(x[2])^HhomeQ(x[1])) >;
    NiceGeneratorsG := &join {imGhomG(Generators(UpperCentralSeries(eQ)[i] 
                       meet imG)): i in [2..#UpperCentralSeries(eQ)] };
    T:=Transversal(eG, sub< eG | GhomeG(NiceGeneratorsG) >);
    return HactG, eGhomG(T) join NiceGeneratorsG;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
Tensor := function(G, H, GactH, HactG, NiceGeneratorsG, NiceGeneratorsH, p)

    nG := #Generators(G);
    nH := #Generators(H);
    n := nG + nH;
    eG, GhomeG, eGhomG := EnumeratedGroup(G, p);

    // If H=G then set eH=eG etc., else enumerate H
    if H cmpeq G then
        eH:=eG; HhomeH:=GhomeG; eHhomH:=eGhomG;
    else
        eH, HhomeH, eHhomH := EnumeratedGroup(H, p);
    end if;

    /* GH:=G*H/J will be a quotient of the free product G*H in which the
       commutator subgroup [G,H] is isomorphic to the tensor product of
       G with H. */

    F := FreeGroup(n);
    GhomF := hom<G->F | [F.i : i in [1..nG]] >;
    HhomF := hom<H->F | [F.i : i in [nG+1..n]] >;

    // The function ProductRelations, defined below, yields generators for J
    rels := ProductRelations( GhomF, HhomF, GhomeG, eGhomG, NiceGeneratorsG, 
                              HhomeH, eHhomH, NiceGeneratorsH, GactH, HactG );

    // The function Remainder completes the function Tensor.
    return Remainder( GhomF, HhomF, GhomeG, eGhomG, NiceGeneratorsG, HhomeH, 
                      eHhomH, NiceGeneratorsH, GactH, HactG, rels, p );

end function;
/* ********************************************************************* */



/* ********************************************************************* */
Remainder := function( GhomF, HhomF, GhomeG, eGhomG, NiceGeneratorsG, HhomeH, 
                       eHhomH, NiceGeneratorsH, GactH, HactG, rels, p );

    // This function completes the function Tensor.
    
    F := Codomain( GhomF );
    G := Domain( GhomeG );
    H := Domain( HhomeH );
    eG := Codomain( GhomeG );
    eH := Codomain( HhomeH );
    nG := #Generators(G);
    nH := #Generators(H);
    n := nG + nH;
    GH, FhomGH := quo<F | rels>;
    GHhomG := hom< GH -> G | [G.i : i in [1..nG]] cat [Id(G) : i in [1..nH]]>;
    eGH, GHhomeGH, eGHhomGH := EnumeratedGroup(GH, p);

    eGtensorH :=
             ncl< eGH | [ GHhomeGH(FhomGH(GhomF(G.i))) : i in [1..nG] ]>
                                  meet
             ncl< eGH | [ GHhomeGH(FhomGH(HhomF(H.i))) : i in [1..nH] ]>;

    GxHmapeGtensorH := map< car<G,H>->eGtensorH | x:->
            ( GHhomeGH(FhomGH(GhomF(x[1]))), GHhomeGH(FhomGH(HhomF(x[2]))) )>;

    eTCentre:=
      Centre(eGH) meet sub<eGH|[GHhomeGH(FhomGH(GhomF(G.i))) : i in [1..nG]]>;

    eTCentrehomG := hom<eTCentre->G | x:-> GHhomG(eGHhomGH(x))>;

    GacteGtensorH := map<car<G,eGtensorH>->eGtensorH |
                                  x:-> x[2]^GHhomeGH(FhomGH(GhomF(x[1])))^-1>;

    /* The remainder of this procedure constructs the homomorphism delta.
    This is done by means of the semidirect product of G with H, which we
    denote by GsdH. */

    S := { GhomF(x)^HhomF(y) = GhomF(HactG(<y,x>)) :
                                     y in Generators(H), x in NiceGeneratorsG }
                      join
        { GhomF(r) : r in Relations(G) } join { HhomF(r) : r in Relations(H) };
    GsdH, FhomGsdH := quo< F | S >;

    // The following is in fact not quite a homomorphism! 
    GHmapGsdH := hom< GH->GsdH | [GsdH.i : i in [1..n]] >;
    eGsdH, GsdHhomeGsdH, eGsdHhomGsdH := EnumeratedGroup(GsdH, p);
    eGhomeGsdH := hom<eG->eGsdH| x:->GsdHhomeGsdH(FhomGsdH(GhomF(eGhomG(x))))>;
    imG := eGhomeGsdH(eG);
    m := IsPrime(p) select #PCGenerators(imG) else #Generators(imG);
    S := [ eGhomG(x) : x in eG, i in [1..m] | eGhomeGsdH(x) eq imG.i ];
    imageGhomG := hom< imG -> G | S >;
    delta := hom< eGtensorH->G |
                    x :-> imageGhomG(GsdHhomeGsdH(GHmapGsdH(eGHhomGH(x)))) >;

    return eGtensorH, delta, eTCentre, eTCentrehomG, GxHmapeGtensorH,
                                                               GacteGtensorH;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
ProductRelations := function( GhomF, HhomF, GhomeG, eGhomG, NiceGeneratorsG, 
                              HhomeH, eHhomH, NiceGeneratorsH, GactH, HactG );

    // This function produces the relations for G*H/J in the function Tensor.

    F := Codomain( GhomF );
    G := Domain( GhomF );
    H := Domain( HhomF );
    rels := { GhomF(LHS(x)*RHS(x)^-1) : x in Relations(G) }
                join { HhomF(LHS(x)*RHS(x)^-1) : x in Relations(H) };
    rels := SetToSequence(rels);
        for z in Generators(G) do
                for x in NiceGeneratorsG do
                for y in NiceGeneratorsH do
                x1:=eGhomG(GhomeG(x));
                y1:=eHhomH(HhomeH(y));
                c:=(GhomF(x1), HhomF(y1)) ^ GhomF(z);
                v:=GhomeG(z^-1*x1*z); v:=eGhomG(v);
                w:=HhomeH(GactH(<z,y1>)); w:=eHhomH(w);
                d:=(HhomF(w),GhomF(v));
                rels:= Append(rels,c*d);
                end for;
                end for;
        end for;
        for z in Generators(H) do
                for x in NiceGeneratorsG do
                for y in NiceGeneratorsH do
                x1:=eGhomG(GhomeG(x));
                y1:=eHhomH(HhomeH(y));
                c:=HhomF(z^-1) * (GhomF(x1),HhomF(y1)) * HhomF(z);
                v:=HhomeH(z^-1*y1*z); v:=eHhomH(v);
                w:=GhomeG(HactG(<z,x1>)); w:=eGhomG(w);
                d:=(HhomF(v),GhomF(w));
                rels:= Append(rels,c*d);
                end for;
                end for;
        end for;
    rels:=SequenceToSet(rels);
    return rels;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
TensorSquare := function(G, p);

    GactG, NiceGensG := SelfAction(G, p);
    eGtensorG, delta := Tensor(G, G, GactG, GactG, NiceGensG, NiceGensG, p);
    return eGtensorG, delta;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
Jgroup := function(G, p);

    _, delta := TensorSquare(G,p);
    return EnumeratedKernel(delta);

end function;
/* ******************************************************************** */



/* ******************************************************************** */
TensorCentre := function(G, p);

    GactG, NiceGensG := SelfAction(G, p);
    _, _, TC, phi := Tensor(G, G, GactG, GactG, NiceGensG, NiceGensG, p);
    return TC, phi;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
TensorCube := function(G, p);

    // TS=TensorSquare and TC=TensorCube 
    eG, GhomeG, eGhomG := EnumeratedGroup(G, p);
    GactG, NiceGensG := SelfAction(G, p);
    eTS, delta, _, _, _, GacteTS :=
         Tensor(G, G, GactG, GactG, NiceGensG, NiceGensG, p );
    TS, TShomeTS, eTShomTS := PresentedGroup( eTS );
    TShomeG := hom< TS->eG | x:->GhomeG(delta(TShomeTS(x))) >;
    TSactG, NiceGensG := Action(GhomeG, TShomeG, p);
    GactTS:=map<car<G,TS>->TS|x:->eTShomTS(GacteTS(<x[1],TShomeTS(x[2])>))>;

    // We now construct a generating set for TS closed under the action of G.
    A := SetToIndexedSet(Generators(eTS));
    i := 1;
    while i le #A do
        y := A[i];
        for x in Generators(eG) do
            Include(~A, GacteTS(<eGhomG(x),y>));
        end for;
        i := i + 1;
    end while;
    NiceGensTS := eTShomTS(A);

    Tcube, delta := Tensor(TS, G,TSactG, GactTS, NiceGensTS, NiceGensG, p);
    return Tcube, delta;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
Jtilde := function(G, p);

    GactG, NiceGensG := SelfAction(G, p);
    eGtensorG, delta, _, _, phi :=
                Tensor(G, G, GactG, GactG, NiceGensG, NiceGensG, p );
    JG := sub< eGtensorG |
           { x : x in Centre(eGtensorG) | delta(x) eq delta(Id(eGtensorG)) }>;

    D := sub< eGtensorG | { phi(x,x) : x in Generators(G) },
                          { phi(x,y)*phi(y,x) : x, y in Generators(G) } >;

    eGsymmetricG, phi := quo< eGtensorG | D >;
    return phi(JG);

end function;
/* ******************************************************************** */



/* ******************************************************************** */
Exterior := function(GhomQ, HhomQ, p);

    G := Domain(GhomQ);
    H := Domain(HhomQ);
    Q := Codomain(GhomQ);
    nG := #Generators(G);
    nH := #Generators(H);
    n := nG+nH;
    eQ, QhomeQ, eQhomQ := EnumeratedGroup(Q, p);
    // If G=Q then set eG=eQ etc., else enumerate G.
    if G cmpeq Q then
        eG:=eQ; GhomeG:=QhomeQ; eGhomG:=eQhomQ;
    else
        eG, GhomeG, eGhomG := EnumeratedGroup(G, p);
    end if;
    // If H=Q then set eH=eQ etc., else enumerate H.
    if H cmpeq Q then
        eH:=eQ; HhomeH:=QhomeQ; eHhomH:=eQhomQ;
    else
        eH, HhomeH, eHhomH := EnumeratedGroup(H, p);
    end if;

    GhomeQ:=hom<G->eQ|x:->QhomeQ(GhomQ(x))>;
    HhomeQ:=hom<H->eQ|x:->QhomeQ(HhomQ(x))>;
    GIH := sub<eQ|GhomeQ(eGhomG(eG)) meet HhomeQ(eHhomH(eH) ) >;

    if (G cmpeq Q) and (H cmpeq Q) then
        HactG, NiceGeneratorsG := SelfAction(G, p);
        GactH := HactG; NiceGeneratorsH := NiceGeneratorsG;
    else
        HactG, NiceGeneratorsG := Action(GhomeQ, HhomeQ, p);
        GactH, NiceGeneratorsH := Action(HhomeQ, GhomeQ, p);
    end if;    

    /* GH:=G*H/J will be a quotient of the free product G*H in which the 
       commutator subgroup [G,H] is isomorphic to the exterior product of G 
       with H. */
    F:=FreeGroup(n);
    GhomF:=hom<G->F | [F.i : i in [1..nG]] >;
    HhomF:=hom<H->F | [F.i : i in [nG+1..n]] >;
    rels := ProductRelations( GhomF, HhomF, GhomeG, eGhomG, NiceGeneratorsG,
                              HhomeH, eHhomH, NiceGeneratorsH, GactH, HactG);

    /* We next add relations (x exterior x') to rels where the images of x and 
       x' are equal in GIH. */
    SG := []; SH := [];
    for x in Generators(GIH) do
        for y in eG do
            if GhomeQ(eGhomG(y)) eq x then x1:=y; break; end if;
        end for;
        Append(~SG,x1);
        for z in eH do
            if HhomeQ(eHhomH(z)) eq x then x2:=z; break; end if;    
        end for;
        Append(~SH,x2);
    end for;
    for i in [1..#Generators(GIH)] do
        x1:=SG[i]; x2:=SH[i];
        rels := rels join {( GhomF(eGhomG(x1)), HhomF(eHhomH(x2)) )};
        for j in [1..#Generators(GIH)] do
            y1:=SG[j]; y2:=SH[j];
            rels := rels join 
              {( GhomF(eGhomG(x1)),HhomF(eHhomH(y2)) )*
               ( GhomF(eGhomG(y1)),HhomF(eHhomH(x2)) )};
        end for;
    end for;

    return Remainder( GhomF, HhomF, GhomeG, eGhomG, NiceGeneratorsG, HhomeH,
                      eHhomH, NiceGeneratorsH, GactH, HactG, rels, p);

end function;
/* ******************************************************************* */



/* ******************************************************************** */
ExteriorSquare := function(G, p);

    GactG, NiceGensG := SelfAction(G, p);
    GhomG := hom<G->G | x:->x >;
    eGexteriorG, delta:= Exterior(GhomG, GhomG, p);
    return eGexteriorG, delta;

end function;
/* ******************************************************************** */
     


/* ******************************************************************** */
Htwo := function(G, p);

    _, delta := ExteriorSquare(G,p);
    return EnumeratedKernel(delta);

end function;
/* ******************************************************************** */



/* ******************************************************************** */
RelativeSchurMultiplier := function(G, N, p);

    // Here N is a normal subgroup of G.
    N1 := Rewrite(G, N);
    GhomG := hom< G->G | x:->x >;
    N1homG := hom< N1->G | x:->x >;
    eGexteriorN, delta, _, _, _, _ := Exterior(GhomG, N1homG, p);
    alpha := hom< Centre(eGexteriorN)->G | x:->delta(x) >;
    return EnumeratedKernel(alpha);

end function;
/* ******************************************************************** */



/* ******************************************************************** */
RelativeEpiCentre := function(G, N, p);

    // Here N is a normal subgroup of G. 
    N1:=Rewrite(G,N);
    GhomG:=hom<G->G|x:->x>;
    N1homG:=hom<N1->G|x:->x>;
    _, _, REC, phi, _, _ := Exterior(GhomG, N1homG, p);
    return REC, phi;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
BaerInvariant := function(G, c, q, p);

    n := #Generators(G);
    F := FreeGroup(n);
    GhomF := hom< G->F | [F.i : i in [1..n]] >;
    Grels := { GhomF(LHS(r)*RHS(r)^-1) : r in Relations(G) };

    G1rels := Grels;
    for i in [1..c] do
        Drels := { r^q : r in G1rels }      join
                 { (r,x) : r in G1rels, x in Generators(F)};
        G1rels := Drels;
    end for;
    Drels := (Drels join { r^q : r in Grels }) diff {Id(F)};
    D, FhomD := quo< F | Drels >;
    eD, DhomeD, eDhomD := EnumeratedGroup(D, p);
    FhomeD := hom< F->eD | x:->DhomeD(FhomD(x)) >;
    N := NormalClosure(eD, sub<eD|FhomeD(Grels)> );
    eG, eDhomeG := quo< eD | N >;
    LCS := LowerCentralSeries(eD);
    d := #LCS;
    eBaer := Kernel(eDhomeG) meet LCS[Minimum(c+1,d)];
    C := UpperCentralSeries(eD)[Minimum(c+1,d)];
    return AbelianGroup(eBaer), eDhomeG(C);

end function;
/* ******************************************************************** */
  


/* ******************************************************************** */
CoveringGroup := function(G, p);

    n:=#Generators(G);
    F:=FreeGroup(n);
    GhomF:=hom<G->F | [F.i : i in [1..n]] >;
    FhomG:=hom<F->G | [G.i : i in [1..n]] >;
    eG, GhomeG, eGhomG := EnumeratedGroup(G,p); 
    FhomeG:=hom<F->eG | x:->GhomeG(FhomG(x)) >;
    Grels:={ GhomF(LHS(x)*RHS(x)^-1) : x in Relations(G) };
    DGrels:={};
    for x in Generators(F) do
        m:=Order(FhomeG(x));
        DGrels:=DGrels join {x^m};
        for r in Grels do
            DGrels:=DGrels join {(r,x)};
        end for;
    end for;
    DGrels:=DGrels diff {Id(F)};
    DG, FhomDG := quo<F | DGrels>;
    DGhomF:=hom<DG->F | [F.i : i in [1..n]]>;
    eDG, DGhomeDG, eDGhomDG := EnumeratedGroup(DG,p); 
    delta:=hom<eDG->eG| x:-> GhomeG(FhomG(DGhomF(eDGhomDG(x))))>;
    eR:=NormalClosure(eDG,sub<eDG|DGhomeDG(FhomDG(Grels))>);
    Htwo := CommutatorSubgroup(eDG,eDG) meet eR;
    DG:=Simplify(DG);
    return DG, delta, Htwo;

end function;
/* ******************************************************************** */




/* ******************************************************************** */
IsCapable := function(G, p);

    DG, delta := CoveringGroup(G, p);
    N:=delta(Centre(Domain(delta)));
    return (#N eq 1);

end function;
/* ******************************************************************** */



/* ******************************************************************** */
IsTerminal := function(G, p);

    DG, delta, SM := CoveringGroup(G, p);
    Dom:=Domain(delta);
    Cod:=Codomain(delta);
    if not(IsNilpotent(Cod)) then
        print " ";
        print " ";
        print "The group is not nilpotent!";
    else
        c:=NilpotencyClass(Cod);
        S:=PrimeDivisors(Order(Cod));
        A:=LowerCentralSeries(Dom)[c+1];
        B:=CommutatorSubgroup(Dom,Dom);
        alpha:=hom<B->Cod | x:->delta(x)>;
        e:=1;
        for i in [1..#S] do
            e:=e*S[i];
        end for;
        SMhomSM:=hom<SM->SM | x:->x^e>;
        C:=SMhomSM(SM);
        return (#(C meet A) eq #A);
    end if;

end function;
/* ******************************************************************** */



/* ******************************************************************** */
ThirdHomologyModP := function(G, p);

if  not(IsPrime(p)) then print "The integer p must be a prime."; d:=0; 

else

    n:=#Generators(G);
    F:=FreeGroup(n);
    GhomF:=hom<G->F | [F.i : i in [1..n]] >;
    FhomG:=hom<F->G | [G.i : i in [1..n]] >;
    Grels := {GhomF(LHS(x)*RHS(x)^-1) : x in Relations(G)};
    Grels := SetToIndexedSet(Grels);
    eG, GhomeG, eGhomG := EnumeratedGroup(G, p);
    GeneratorsOfR:=car<eG,Grels>;
    Phi:=map<GeneratorsOfR->F|x:->x[2]^GhomF(eGhomG(x[1]^-1))>;
    CommutatorRels := 
        {(Phi(<x,Grels[i]>),Grels[j]) : x in eG, i, j in [1..#Grels] | i ge j}
                      join
        {r^p : r in Grels}; 

    // We set D:=F/[R,R]R^p, and S=R/[R,R]R^p. 
    D, FhomD := quo<F|CommutatorRels>;
    eD, DhomeD, eDhomD := EnumeratedGroup(D, p);
    
    eS:=NormalClosure(eD,sub<eD|DhomeD(FhomD(Grels))>);
    S, ShomeS, eShomS := PresentedGroup(eS);
    ShomeD:=hom<S->eD|x:->ShomeS(x)>;
    GhomeD:=hom<G->eD|x:->DhomeD(FhomD(GhomF(x)))>;
    GactS:=map<car<G,S>->S|x:->
                         eShomS( GhomeD(x[1])*ShomeD(x[2])*GhomeD(x[1]^-1) )>;

    NiceGeneratorsS:=Generators(S); 
    SactG := map<car<S,G>->G|x:->x[2]>;
    X, NiceGeneratorsG := SelfAction(G, p);

    mS := #Generators(S);
    mG := #Generators(G);
    m := mS+mG;
    
    GensG := PCGenerators(eG);
    Bas := {};
    V := sub<eS| Id(eS)>;
    for x in PCGenerators(eS) do
        x1 := eShomS(x);
        for y in PCGenerators(eG) do
            y1:=eGhomG(y);
            NewV:=sub<eS| V, x*ShomeS(GactS(y1,x1))>;
            if #NewV gt #V then V:=NewV; Bas:=Bas join {<x1,y1>}; end if;
        end for;
    end for;
 
    /* SG:=S*G/J will be a quotient of the free product S*G in which the 
       commutator subgroup [S,G] is isomorphic to the tensor product of S 
       with G. */
      
    E:=FreeGroup(m);
    ShomE:=hom<S->E | [E.i : i in [1..mS]] >;
    GhomE:=hom<G->E | [E.i : i in [mS+1..m]] >;
    rels := ProductRelations(ShomE, GhomE, ShomeS, eShomS, NiceGeneratorsS,
                             GhomeG, eGhomG, NiceGeneratorsG, SactG, GactS); 

    rels := rels join {(ShomE(w[1]),GhomE(w[2])) : w in Bas};
    SG, EhomSG := quo<E|rels>;
    eSG, _, _ := EnumeratedGroup(SG, p);

    // X will be the subgroup [R,G]/[R,R]R^p of S=R/[R,R]R^p.
    X := {ShomeS(GactS(eGhomG(g),s^-1)*s) : g in eG, s in Generators(S)};
    X:=NormalClosure(eS,sub<eS|X>);

    a:= #eSG * (p^#Bas) /((#eS)*(#eG)*(#X));
    for i in [0..100] do 
        if p^i eq a then d:=i; break; end if; 
    end for;

end if;

    return d;

end function;
/* ******************************************************************** */