FreeGp:=FreeGroup("A","B","C");
Af:=FreeGp.1;
Bf:=FreeGp.2;
Cf:=FreeGp.3;

Rels:=[
(Bf^-1*Af)^3,

(Cf*Bf^-1*Cf)^3,

(Af^-3*Bf)^3,

Bf^-1*Cf*Bf^-1*Af^-1*Bf^-1*Cf*Af*Cf^2*Bf^-1*Cf^-1*Bf*Cf^-1,

Bf^-1*Af^-1*Cf^-1*Bf*Af*Bf*Af^-2*Bf*Af^-2*Bf^-1*Cf*Bf^-1,

Af^-1*Bf^2*Af^-1*Cf*Bf^-1*Cf*Bf^-1*Cf*Af*Bf*Cf^-1*Af*Bf^-1*Af^-1,

Af^2*Bf^-1*Af*Bf*Af^-1*Bf*Cf^-1*Bf*Af^2*Bf^-1*Af*Cf^2*Bf^-1*Cf,

Af^-1*Bf^-1*Cf*Af^-1*Bf^2*Af^-1*Cf*Bf^-1*Cf*Af*Bf*Cf*Bf^-1*Cf*Bf^-1*Cf*Bf*Af^-1,

Cf*Bf*Cf^-1*Af*Bf^-1*Af^-3*Bf^-1*Af^-1*Bf^-1*Cf*Bf*Af^-2*Bf^-1*Cf*Bf^-1*Af^-1,

Cf^-1*Bf*Cf^-1*Af*Bf*Af*Bf*Af^-1*Cf^-1*Bf^2*Cf^-1*Af*Bf^-1*Af^-2*Cf^2*Bf^-1,

Af*Bf^-1*Af^-1*Bf^-2*Af^-1*Cf^-1*Bf*Af^2*Bf^-1*Cf*Bf^-1*Cf*Bf^-1*Cf*Af*Cf*Bf^-1*Cf,

Cf*Bf^-1*Cf*Af*Bf*Af*Bf^-2*Af^2*Bf*Af^-1*Cf^-1*Bf*Cf^-1*Af^-1*Bf*Af^-2*Bf^-1,

Cf*Bf*Af^-2*Bf^-1*Cf*Af^-2*Bf^-1*Cf*Af*Bf^2*Cf*Bf^-1*Cf*Bf^-1*Cf*Bf^-1*Af^-1*Bf^-1,

Cf*Bf^-1*Cf*Af*Bf^-1*Af^-1*Bf^-1*Af^-1*Cf^-1*Bf*Cf^-2*Af*Bf*Cf^-1*Bf*Af*Bf*Af^-1*Bf^2,

Bf*Cf^-1*Bf*Cf^-2*Af^2*Bf*Af^-1*Cf*Bf^-1*Cf^-1*Bf^-1*Af*Bf^-1*Af^-1*Bf^-1*Cf*Af*Bf*Af^-1,

Bf*Af^-2*Cf*Bf*Cf^-1*Af*Bf^-1*Af^-3*Cf^-1*Bf*Cf*Bf^-1*Cf^2*Af*Bf*Af^-1*Bf,

Cf^-1*Af*Bf^-1*Af^-1*Bf*Af^-1*Cf*Af*Bf*Af*Bf^-2*Af*Cf^-1*Bf*Cf^-1*Bf*Cf^-2*Af*Bf^-1*Af^-1,

Bf*Af*Bf*Af^-1*Bf^2*Af^-2*Cf*Bf^-1*Cf^2*Af^-1*Cf*Bf^-1*Cf*Bf^-1*Cf*Af*Cf*Bf^-1*Cf,

Cf^-1*Af*Bf*Af*Bf*Af^-1*Cf^-1*Bf*Cf^-1*Af^-1*Bf^-1*Af^-1*Bf^-1*Cf*Af*Cf*Bf^-1*Af^-1*Bf^-1*Cf*Af*Bf,

Af^-1*Cf^2*Bf^-1*Cf*Bf^-1*Cf*Af*Bf*Af^-1*Cf^-1*Bf*Cf*Bf^-1*Cf^-1*Bf*Cf^-1*Af*Bf*Af^-1*Bf^2*Af^-1,

Cf*Bf^-1*Cf*Bf^-1*Cf*Af*Cf^-1*Bf*Cf^-1*Bf^-1*Cf*Af*Cf^-1*Af*Bf^-1*Af^-1*Bf^-1*Af^-1*Cf^-1*Bf*Cf^-2*Af,

Bf*Cf^-1*Bf*Af*Bf^2*Cf^-1*Bf*Af*Bf*Af^-1*Bf^2*Af*Bf*Af^-2*Bf^-1*Cf*Af*Bf*Af^-1*Cf^-1,

Af^-1*Bf*Af^-2*Bf^-1*Cf*Bf*Cf^-1*Af*Bf^-1*Af^-2*Bf^-1*Af^2*Bf*Af^-1*Cf^-1*Bf*Cf^-1*Af^2*Bf^-1,

Bf^-1*Af^-1*Cf^-1*Bf*Cf^-1*Af*Bf^-1*Af^-1*Bf*Af^-1*Cf*Af*Bf^2*Af^2*Bf^-1*Af*Cf*Bf^-1*Cf^2*Af^-1,

Bf^-1*Af^-1*Bf*Af^-2*Cf*Bf^-1*Cf*Bf*Cf^-3*Bf*Cf^-1*Bf^-1*Cf*Af*Cf^-2*Bf*Cf^-1*Bf^-1*Af^2,

Bf*Af^-2*Cf^-1*Bf*Cf^-2*Af^4*Bf*Af^-1*Cf*Bf^-1*Cf^-1*Af*Bf^-1*Af^-1*Bf^-2*Af^-1*Cf^-1*Bf,

Bf^-2*Af^-1*Bf^-1*Cf*Af*Bf*Cf^-1*Af^-1*Cf^-1*Bf*Cf^-1*Bf*Cf^-1*Bf^-1*Af^-1*Cf^-1*Bf^-1*Af^-1*Cf^-1*Bf*Cf^-1*Bf*Af^2,

Cf*Bf*Cf^-1*Af*Bf^-1*Af^-2*Bf*Cf^-1*Bf*Cf^-2*Af^-1*Cf^-1*Bf*Cf^-1*Bf*Af*Bf*Cf*Bf^-1*Cf*Bf^-1*Cf*Af,

Af^-1*Cf*Bf*Cf^-1*Af*Bf^-1*Af^-2*Bf^-1*Af^-1*Cf^-1*Bf*Af*Cf*Af*Bf^2*Cf^-1*Bf^2*Af^-2*Cf*Bf^-1*Cf^2*Af^-1,

Bf*Cf^-1*Bf^-1*Cf*Af*Bf^-2*Af^-1*Cf^-1*Bf*Af*Cf*Af*Bf*Cf*Af*Bf*Cf*Bf^-1*Cf*Af*Bf*Cf^-1*Af^-1*Cf^-1*Bf*Af*Bf*Af*Cf^-1,

Bf*Af^-2*Bf^-1*Cf*Bf^-1*Af^-1*Cf*Bf^-1*Cf*Bf^-1*Cf*Af*Cf^-2*Bf*Cf^-1*Af^-1*Cf*Bf^-1*Cf*Bf^-1*Cf*Af*Cf^-1*Bf*Cf^-2*Bf^-1*Af^-1*Cf^-1*Bf,

Bf*Af^-1*Cf^-1*Bf*Cf^2*Bf^-1*Cf*Bf*Cf^-2*Bf^-1*Af^-1*Cf^-1*Bf*Af^2*Bf^-1*Af^-1*Bf^-1*Af^-1*Bf*Af^-2*Bf^-1*Cf*Af^-2*Cf*Bf^-1*Cf*Bf*Cf^-2,

Bf*Af*Bf*Af^-2*Bf^-1*Cf*Af*Bf*Cf*Bf^-1*Cf*Bf*Af^2*Bf^-2*Af*Bf^-1*Af^-1*Bf^-1*Cf^-1*Bf*Cf^-1*Bf^-1*Af*Bf*Af^-2*Bf^-1*Cf*Bf^-1*Af^-1*Cf^-1*Bf*Cf^-3*Bf*Cf^-1*Bf*Cf^-1,

Bf*Af*Bf*Af^-2*Bf^-1*Cf*Af*Bf*Cf^2*Bf^-1*Cf*Af^-1*Bf*Af^-2*Cf*Bf^-1*Cf*Bf*Cf^-2*Bf*Cf^2*Bf^-1*Cf^-1*Bf*Cf^-1*Af^2*Bf^-1*Af^2*Bf^-1*Cf*Af*Bf*Af*Bf^-1*Af^-1*Bf*Af^-1*Bf^-1*Af^-1*Bf^-2*Af^-2*Bf^2*Af^-1*Bf*Cf^-1*Af*Bf^-2*Af^-1*Bf*Af^-2*Bf^-1*Cf^-1*Af^-1*Cf^-1*Bf*Af*Bf*Af*Cf^-2*Bf*Cf^-1*Af^-2*Cf^-1*Bf*Cf^-1*Bf^-1*Cf
];

# Finitely presented version of $\bar\Gamma$

GammaBarFP:=FreeGp/Rels;
A:=GammaBarFP.1;
B:=GammaBarFP.2;
C:=GammaBarFP.3;

# Important subgroup of $\bar\Gamma$

# These are Donald's
#index3aFP:=Group(A, C, B*A*B^-1, B^3);
#index3bFP:=Group(C, A^-1*B^-1, B^-1*A^-1, A^2*B^-1);

# These are Tim's
index3aFP:=Group(A,C,B*A*B^-1);
index3bFP:=Group(C,B*A,A*B,B*A^-2);

index3aFCA:=FactorCosetAction(GammaBarFP,index3aFP);
index3bFCA:=FactorCosetAction(GammaBarFP,index3bFP);

# List of finite order elements

# These are Donald's
#T1:=B*A^-1;
#T2:=A*B^-1*A^-1*B^-1;
#T3:=C*B^-1*A*B^-1*A^-1*B^-1*C*A*B*B*A^-1*B*B*A^-1*B*A*B*A^-1*B;

# These are Tim's

T1:=A*B^-1;
T2:=B*A*B*A^-1;
T3:=C*B^-1*C;

FOList:=[T1,T2,T3];

# Check the indexes of the subgroup

Index(GammaBarFP,index3aFP) = 3;
Index(GammaBarFP,index3bFP) = 3;

# Check normality of the subgroup

IsNormal(GammaBarFP,index3aFP);
not IsNormal(GammaBarFP,index3bFP);

# Check that normalizers are as expected

Normalizer(GammaBarFP,index3aFP) = GammaBarFP;
Normalizer(GammaBarFP,index3bFP) = index3bFP;

# Calculate abelianizations of the subgroups

AbelianInvariants(GammaBarFP) = [3,4];
AbelianInvariants(index3aFP) = [4,31];
AbelianInvariants(index3bFP) = [2,3,4,4];

# Check which elements of finite order belong to which conjugates of
# which finite index subgroups 

List(FOList,k->k^index3aFCA);
List(FOList,k->k^index3bFCA);

# Check that the first three index~3 groups are torsion free

ForAll(FOList,k->MovedPoints(k^index3aFCA) = [1..3]);
ForAll(FOList,k->MovedPoints(k^index3bFCA) = [1..3]);

# Function to calculate fundamental group

FundGp:=function(G,FOList)
  local e1,e2,e3,e4,FCA,GFO,G0;
  e1:=GeneratorsOfGroup(GammaBarFP);
  e2:=Concatenation(e1,List(e1,Inverse));
  Add(e2,One(G));
  e3:=ListX(e2,e2,\*);
  FCA:=FactorCosetAction(GammaBarFP,G);
  GFO:=Concatenation(
    List(FOList,fo->
      List(Difference([1..Index(GammaBarFP,G)],MovedPoints(fo^FCA)),
        j -> fo^First(e3,e->j^(e^FCA) = 1)
      )
    )
  );
  Print("Size(GFO)=",Size(GFO),"\n");
  e1:=GeneratorsOfGroup(G);
  e2:=Concatenation(e1,List(e1,Inverse));
  Add(e2,One(G));
  e3:=ListX(e2,e2,\*);
  e4:=ListX(e2,e3,\*);
  G0:=Group(ListX(e4,GFO,
    function(elt,fo) return fo^elt; end
  ));
  Print(ForAll(GFO,fo->fo in G0),"\n");
  Print(IsNormal(G,G0),"\n");
  return FactorGroup(G,G0);
end;

StructureDescription(FundGp(GammaBarFP,FOList)) = "C4";

# This calculates the automorphism group of the corresponding surface
AutGp:=G->FactorGroup(Normalizer(GammaBarFP,G),G);

StructureDescription(AutGp(index3aFP)) = "C3";
StructureDescription(AutGp(index3bFP)) = "1";