# The generators used here are Donald's.
#
# Although they are in lower case, they are not inverses.
#
# Relative to the generators in ~steger/donald/SU21/C11/C11-0-FewElts
# the generators used here are:
#
#   j   <-->  J
#   b   <-->  A    (and b^-1  <-->  a=A^-1)
#   u   <-->  K1
#   v   <-->  K8

FreeGp:=FreeGroup("j","u","v","b");
xj:=FreeGp.1;
xu:=FreeGp.2;
xv:=FreeGp.3;
xb:=FreeGp.4;

Rels:=[
xu^4,
xv^8,
Comm(xu,xj),
Comm(xv,xj),
xj^-3*xv^2,
xu*xv*xu*xv^-1*xu*xv^-1,

(xb*xj)^2*(xv*xu^2)^-1,
Comm(xb,xv*xu^2),
xb^3,
(xb*xv*xu^3)^3

# These are the omitted relations
#,
#(xb*xu*xj)^4,
#(xb*xu^2*xb*xu*xj)^3,
#(xb*xu^2*xv*xu*xj^-1*xb*xu*xj)^3,
#(xb*xu^3*xv*xj^-1*xb*xu*xj)^3

];

GammaBarFP:=FreeGp/Rels;
j:=GammaBarFP.1;
u:=GammaBarFP.2;
v:=GammaBarFP.3;
b:=GammaBarFP.4;

ix3:=Group(u,v,b*j,j*b);
ix3FCA:=FactorCosetAction(GammaBarFP,ix3);
Index(GammaBarFP,ix3) = 3;
AbelianInvariants(ix3)=[2,2];

ix864:=GroupWithGenerators([v*u*b*j*u^-1, u^-1*j^-1*b*j^2, u^2*v*b*u*j^-2]);
ix864FCA:=FactorCosetAction(GammaBarFP,ix864);
Index(GammaBarFP,ix864) = 864;
AbelianInvariants(ix864)=[0,0];

#####
#
# Calculation of abelianizations for ix864, and its subgroups of
# index~21, 7, and 3, called
# 
#   ix18144, ix6048, and ix2592
#   
# Calculation also of the actions on these abelianizations by the
# normalizers of these groups.

ix864MAQ:=MaximalAbelianQuotient(ix864);    
ix864Ab:=Image(ix864MAQ);
IsAbelian(ix864Ab);

ix864FCA:=FactorCosetAction(GammaBarFP,ix864);;
ix864Norm:=Normalizer(GammaBarFP,ix864);
StructureDescription(ix864Norm/ix864) = "C3";

IsOne(ix864Ab.1);
toVec864:=function(x)
  local er,vec,l,j,ix;
  er:=ExtRepOfObj(x);
  l:=Size(er)/2;
  vec:=0*[1..2];
  for j in 2*[0..l-1] do
    ix:=er[j+1]-1;
    if ix > 0 then
      vec[ix]:=vec[ix]+er[j+2];
    fi;
  od;
  return vec;
end;

gens864:=[ix864.1*ix864.2^-1,ix864.3];
mtxInv864:=TransposedMat(
  List(gens864,x->toVec864(x^ix864MAQ))
)^-1;

j*v^2 in Normalizer(GammaBarFP,ix864);
ZMtx864:=mtxInv864*TransposedMat(
  List(gens864,x->toVec864((x^(j*v^2))^ix864MAQ))
);
StructureDescription(Group(ZMtx864)) = "C3";
IsZero(ZMtx864^2+ZMtx864+One(ZMtx864));

ix288:=GroupWithGenerators([j*v^2,j^2*b*u^-1*b,j^2*b^-1*u*j]);
ix288=Normalizer(GammaBarFP,ix864);
ix288FCA:=FactorCosetAction(GammaBarFP,ix288);
Index(GammaBarFP,ix288) = 288;
	
ix18144:=GroupWithGenerators(
  v^2*b*j^-1*v^-2*b^-1*j, v^2*b^-1*j^-1*v^-2*b*j, v*b*j^-1*v^-2*b^-1*v*j, 
  v*b^-1*j^-1*v^-2*b*v*j, v^-2*b*j*v^2*b^-1*j^-1, v^-2*b^-1*j*v^2*b*j^-1, 
  v^-1*b^-1*j*v^2*b*v^-1*j^-1, b*j^-1*v^-2*b^-1*v^2*j, b^-1*j*v^2*b*v^-2*j^-1, 
  j^2*u*v*b*v^-2*b*u, j*u*v*b*v^-2*b*u*j, j*v*b^-1*v*b*u^-1*b^-1*u*b^-1, 
  j*v^-1*b^-1*j*v^2*b*v^-1*j^-2, j*b*v*b*v^-2*b*v*b
);
  
ix18144FCA:=FactorCosetAction(GammaBarFP,ix18144);;
IsNormal(GammaBarFP,ix18144);
Index(ix864,ix18144) = 21;
StructureDescription(GammaBarFP/ix18144) = "C3 x PSU(3,3)";
StructureDescription(ix3/ix18144) = "PSU(3,3)";
StructureDescription(ix864/ix18144) = "C7 : C3";
AbelianInvariants(ix18144) = [0,0,0,0,0,0,0,0,0,0,0,0,0,0];

ix18144MAQ:=MaximalAbelianQuotient(ix18144);    
ix18144Ab:=Image(ix18144MAQ);
IsAbelian(ix18144Ab);

toVec18144:=function(x)
  local er,vec,l,j;
  er:=ExtRepOfObj(x);
  l:=Size(er)/2;
  vec:=0*[1..14];
  for j in 2*[0..l-1] do
    vec[er[j+1]]:=vec[er[j+1]]+er[j+2];
  od;
  return vec;
end;

gens18144:=GeneratorsOfGroup(ix18144);
mtxInv18144:=TransposedMat(
  List(gens18144,x->toVec18144(x^ix18144MAQ))
)^-1;

jMtx18144:=mtxInv18144*TransposedMat(
  List(gens18144,x->toVec18144((x^(j^-1))^ix18144MAQ))
);
uMtx18144:=mtxInv18144*TransposedMat(
  List(gens18144,x->toVec18144((x^(u^-1))^ix18144MAQ))
);
vMtx18144:=mtxInv18144*TransposedMat(
  List(gens18144,x->toVec18144((x^(v^-1))^ix18144MAQ))
);
bMtx18144:=mtxInv18144*TransposedMat(
  List(gens18144,x->toVec18144((x^(b^-1))^ix18144MAQ))
);

ix18144PU:=
  GroupWithGenerators([uMtx18144,vMtx18144,jMtx18144*bMtx18144]); 
Size(ix18144PU) = 6048;
StructureDescription(ix18144PU) = "PSU(3,3)";

ix18144PUCT:=CharacterTable(ix18144PU);
ix18144PUCC:=ConjugacyClasses(ix18144PUCT);;
List(ix18144PUCC,c->Trace(Representative(c)))
  = 2*Irr(ix18144PUCT)[3]{[1..14]};
Irr(ix18144PUCT)[3];
  
ix18144PUC3:=
  GroupWithGenerators([jMtx18144,uMtx18144,vMtx18144,bMtx18144]);
StructureDescription(ix18144PUC3) = "C3 x PSU(3,3)";

ix18144C3:=GroupWithGenerators([jMtx18144*vMtx18144^2]);
ix18144C3 = Center(ix18144PUC3);
StructureDescription(ix18144C3) = "C3";
IsZero((jMtx18144*vMtx18144^2)^2+jMtx18144*vMtx18144^2+One(jMtx18144));

hom18144:=GroupHomomorphismByImages(GammaBarFP,ix18144PUC3,
  [j,u,v,b],[jMtx18144,uMtx18144,vMtx18144,bMtx18144]);

ix6048:=GroupWithGenerators([
  j^-1*u^-1*b*j^2*u*v*b^-1*v*u*b, j^2*b*u^-1*b*j^-1*u^-1*b*j^2, 
  j*b^-1*u^-1*v^-1*b^-1*v^-1*j*v*b*j^-1*b*u^-1*j^-1,
  j*b^-1*u^-1*v^-1*b^-1*v^-1*j^-1*u*b*u^-1*v*b,
  j*b^-1*u^-1*v^-1*b^-1*v^-2*b*v^2*b*u^-1*j^-1,
  j*b^-1*u^-1*v^-1*b^-1*v^-1*b^2*v^-1*b*j^2*b^-1,
  u*b^-1*u*v^-1*b*u^-1*b^-1*v^3*b
]);
Index(GammaBarFP,ix6048) = 6048;
Index(ix6048,ix18144) = 3; 
AbelianInvariants(ix6048) = [0,0,0,0,0,0];

ix6048FCA:=FactorCosetAction(GammaBarFP,ix6048);;
ix6048Norm:=Normalizer(GammaBarFP,ix6048);
StructureDescription(ix6048Norm/ix6048) = "C6 x C3";

ix6048MAQ:=MaximalAbelianQuotient(ix6048);
ix6048Ab:=Image(ix6048MAQ);

toVec6048:=function(x)
  local er,vec,l,j;
  er:=ExtRepOfObj(x);
  l:=Size(er)/2;
  vec:=0*[1..6];
  for j in 2*[0..l-1] do
    vec[er[j+1]]:=vec[er[j+1]]+er[j+2];
  od;
  return vec;
end;

gens6048:=GeneratorsOfGroup(ix6048){[1,2,3,5,6,7]};
mtxInv6048:=TransposedMat(
  List(gens6048,x->toVec6048(x^ix6048MAQ))
)^-1;

Z1Mtx6048:=mtxInv6048*TransposedMat(
  List(gens6048,x->toVec6048((x^((b*v*b*j)^-1))^ix6048MAQ))
);
Z2Mtx6048:=mtxInv6048*TransposedMat(
  List(gens6048,x->toVec6048((x^((v*b*v^-1)^-1))^ix6048MAQ))
);

StructureDescription(Group(Z1Mtx6048,Z2Mtx6048)) = "C6 x C3";
Eigenvalues(Field(E(3)),Z1Mtx6048^2) = [E(3),E(3)^2];

ix2592:=GroupWithGenerators( [ b*j^-2*b^-1*v*b,
  v^2*b*j^-1*v^-2*b^-1*j, v*b*j^-1*v^-2*b^-1*v*j,
  v*b^-1*j^-1*v^-2*b*v*j, v^-2*b*j*v^2*b^-1*j^-1,
  v^-2*b^-1*j*v^2*b*j^-1, v^-1*b^-1*j*v^2*b*v^-1*j^-1,
  b*j^-1*v^-2*b^-1*v^2*j, j*u*v*b*v^-2*b*u*j ]);

ix2592FCA:=FactorCosetAction(GammaBarFP,ix2592);;
Index(GammaBarFP,ix2592) = 2592;
IsSubgroup(ix2592,ix18144);
IsNormal(ix2592,ix18144);
Index(ix2592,ix18144) = 7;
IsSubgroup(ix864,ix2592);
Index(ix864,ix2592) = 3;
IsNormal(ix864,ix2592);

ix2592FCA:=FactorCosetAction(GammaBarFP,ix2592);;
ix2592Norm:=Normalizer(GammaBarFP,ix2592);
StructureDescription(ix2592Norm/ix2592) = "C3 x C3";

AbelianInvariants(ix2592) = [ 0, 0, 7 ];
ix2592MAQ:=MaximalAbelianQuotient(ix2592);
ix2592Ab:=Image(ix2592MAQ);
IsOne(ix2592Ab.1^7);

toVec2592:=function(x)
  local er,vec,l,j,ix;
  er:=ExtRepOfObj(x);
  l:=Size(er)/2;
  vec:=0*[1..2];
  for j in 2*[0..l-1] do
    ix:=er[j+1]-1;
    if ix > 0 then
      vec[ix]:=vec[ix]+er[j+2];
    fi;
  od;
  return vec;
end;

gens2592:=GeneratorsOfGroup(ix2592){[1,8]};
mtxInv2592:=TransposedMat(
  List(gens2592,x->toVec2592(x^ix2592MAQ))
)^-1;

StructureDescription(Normalizer(GammaBarFP,ix2592)/ix2592) = "C3 x C3";

Z1Mtx2592:=mtxInv2592*TransposedMat(
  List(gens2592,x->toVec2592((x^((j*v^2)^-1))^ix2592MAQ))
);
j*v^2 in Normalizer(GammaBarFP,ix2592);
not j*v^2 in ix3;
IsZero(Z1Mtx2592^2+Z1Mtx2592+One(Z1Mtx2592));

Z2Mtx2592:=mtxInv2592*TransposedMat(
  List(gens2592,x->toVec2592((x^((j^2*b*u^-1*b)^-1))^ix2592MAQ))
);
j^2*b*u^-1*b in Normalizer(GammaBarFP,ix2592);
j^2*b*u^-1*b in ix3;
IsOne(Z2Mtx2592);
	
StructureDescription(Group(Z1Mtx2592,Z2Mtx2592)) = "C3";

Index(ix864Ab,
  Group(List(GeneratorsOfGroup(ix18144),x->x^ix864MAQ))) = 3;
Index(ix864Ab,
  Group(List(GeneratorsOfGroup(ix6048),x->x^ix864MAQ))) = 1;
Index(ix864Ab,
  Group(List(GeneratorsOfGroup(ix2592),x->x^ix864MAQ))) = 3;
Index(ix2592Ab,
  Group(List(GeneratorsOfGroup(ix18144),x->x^ix2592MAQ))) = 7;
Index(ix6048Ab,
  Group(List(GeneratorsOfGroup(ix18144),x->x^ix6048MAQ))) = 3;

#####
  
J:=j; A:=b;

K1:=u;
K2:=u^2;
K3:=u^3;

K4:=v*u^2*v*j^9;
K5:=u*v*u^2*v*j^9;
K6:=u^2*v*u^2*v*j^9;
K7:=u^3*v*u^2*v*j^9; 

K8:=v;
K9:=u*v;
K10:=u^2*v;
K11:=u^3*v;

K12:=u^2*v*u;
K13:=u*v*u;
K14:=u^3*v*u;
K15:=v*u;

K16:=u*v*u^2;
K17:=u^2*v*u^2;
K18:=u^3*v*u^2;
K19:=v*u^2;

K20:=u^2*v*u^3;
K21:=u*v*u^3;
K22:=u^3*v*u^3;
K23:=v*u^3;

KList:=[One(K1),K1,K2,K3,K4,K5,K6,K7,K8,K9,K10,K11,
        K12,K13,K14,K15,K16,K17,K18,K19,K20,K21,K22,K23];

NoTors:=x->MovedPoints(x^ix864FCA)=[1..864];
NoTors288:=x->MovedPoints(x^ix288FCA)=[1..288];

ListX([1..24],[0..11],
  function(kix,jix)
    if kix=1 and jix=0 then
      return false;
    else
      return true;
    fi;
  end,
  function(kix,jix)
    return not NoTors288(KList[kix]*j^jix);
  end,
  function(kix,jix)
    return KList[kix]*j^jix;
  end
) = [J^4,J^8];

List([J^4,J^8],
  x->288-Size(MovedPoints(x^ix288FCA))) = [288,288];

ForAll(ListX([1..24],[0..11],
  function(kix,jix)
    if kix=1 and jix=0 then return true; fi;
        return NoTors(KList[kix]*j^jix);
  end),
  x->x
);

NoTors(K20*J^5*A^-1*K12*J);
NoTors(K17*J^7*A^-1*K20*J^8*A^-1*K1*J^4);

TorsList:=[
  K20*J^5*A^-1*K17,
  K20*J^5*A^-1*K5*J^4,
  K20*J^5*A^-1*K21*J^2,
  K20*J^6*A^-1*K23*J^5,
  K20*J^6*A^-1*K4*J^7,
  K20*J^6*A^-1*K16*J^3,
  K18*J^7*A^-1*K12*J^2,
  K18*J^7*A^-1*K2*J^4,
  K18*J^7*A^-1*K9*J^3,
  K17*J^7*A^-1*K2*J^3,
  K17*J^7*A^-1*K12*J,
  K17*J^7*A^-1*K9*J^2,
  K18*J^6*A^-1*K8*J^3,
  K18*J^6*A^-1*K1*J^4,
  K18*J^6*A^-1*K13*J^2,
  K17*J^8*A^-1*K17*J,
  K17*J^8*A^-1*K5*J^5,
  K17*J^8*A^-1*K21*J^3,
  K17*J^7*A^-1*K20*J^8*A^-1*K13*J^2,
  K17*J^7*A^-1*K20*J^8*A^-1*K8*J^3,
  K17*J^7*A^-1*K20*J^8*A^-1*K1*J^4,
  
  K20*J^5*A^-1*K12*J,
  K18*J^7*A^-1*K21*J^3,
  
  K20*J^5*A^-1*K21*J^2,
  K20*J^5*A^-1*K16*J^2,
  K20*J^5*A^-1*K9*J^2,
  K20*J^5*A^-1*K13*J^2,
  K20*J^6*A^-1*K16*J^3,
  K20*J^6*A^-1*K21*J^3,
  K20*J^6*A^-1*K13*J^3,
  K20*J^6*A^-1*K9*J^3,
  
  K20*J^5*A^-1*K21*J^3,
  K20*J^5*A^-1*K16*J^3,
  K20*J^5*A^-1*K9*J^3,
  K20*J^5*A^-1*K13*J^3];;
  
ForAll(TorsList,NoTors);
List(TorsList,x->Order(x^ix864FCA));

Filtered(TorsList,x->not NoTors288(x)) = 
  [K20*J^5*A^-1*K12*J,K18*J^7*A^-1*K21*J^3];
List([K20*J^5*A^-1*K12*J,K18*J^7*A^-1*K21*J^3],
  x->288-Size(MovedPoints(x^ix288FCA))) = [6,6];
List([K20*J^5*A^-1*K12*J,K18*J^7*A^-1*K21*J^3],
  x->Difference([1..288],MovedPoints(x^ix288FCA)));

(K20*J^5*A^-1*K12*J)^(b*v) in ix288;
(K18*J^7*A^-1*K21*J^3)^(b*v) in ix288;
J^4 in ix288;
J^8 in ix288;

# This shows that ix288 is generated by its finite order elements,
# hence that the corresponding surface is simply connected.
Group(Concatenation(ListX(
  GeneratorsOfGroup(ix288),GeneratorsOfGroup(ix288),
  [(K20*J^5*A^-1*K12*J)^(b*v),J^4],
  function(x1,x2,y)
   return [y^(x1*x2),y^(x1*x2^-1),y^(x1^-1*x2),y^(x1^-1*x2^-1),
     y^(x2*x1),y^(x2*x1^-1),y^(x2^-1*x1),y^(x2^-1*x1^-1)];
   end
))) = ix288;

Index(ix288,Group(Concatenation(ListX(
  GeneratorsOfGroup(ix288),GeneratorsOfGroup(ix288),GeneratorsOfGroup(ix288),
  [(K20*J^5*A^-1*K12*J)^(b*v)],
  function(x1,x2,x3,y)
   return [
     y^(x1*x2*x3),y^(x1*x2^-1*x3),y^(x1^-1*x2*x3),y^(x1^-1*x2^-1*x3),
     y^(x2*x1*x3),y^(x2*x1^-1*x3),y^(x2^-1*x1*x3),y^(x2^-1*x1^-1*x3),
     y^(x1*x3*x2),y^(x1*x3*x2^-1),y^(x1^-1*x3*x2),y^(x1^-1*x3*x2^-1),
     y^(x2*x3*x1),y^(x2*x3*x1^-1),y^(x2^-1*x3*x1),y^(x2^-1*x3*x1^-1),
     y^(x3*x1*x2),y^(x3*x1*x2^-1),y^(x3*x1^-1*x2),y^(x3*x1^-1*x2^-1),
     y^(x3*x2*x1),y^(x3*x2*x1^-1),y^(x3*x2^-1*x1),y^(x3*x2^-1*x1^-1),
     y^(x1*x2*x3^-1),y^(x1*x2^-1*x3^-1),y^(x1^-1*x2*x3^-1),y^(x1^-1*x2^-1*x3^-1),
     y^(x2*x1*x3^-1),y^(x2*x1^-1*x3^-1),y^(x2^-1*x1*x3^-1),y^(x2^-1*x1^-1*x3^-1),
     y^(x1*x3^-1*x2),y^(x1*x3^-1*x2^-1),y^(x1^-1*x3^-1*x2),y^(x1^-1*x3^-1*x2^-1),
     y^(x2*x3^-1*x1),y^(x2*x3^-1*x1^-1),y^(x2^-1*x3^-1*x1),y^(x2^-1*x3^-1*x1^-1),
     y^(x3^-1*x1*x2),y^(x3^-1*x1*x2^-1),y^(x3^-1*x1^-1*x2),y^(x3^-1*x1^-1*x2^-1),
     y^(x3^-1*x2*x1),y^(x3^-1*x2*x1^-1),y^(x3^-1*x2^-1*x1),y^(x3^-1*x2^-1*x1^-1)
    ];
   end
))));

#####
ix1:=Group(u^-1,v^-1,j^-1,b^-1);
ix4:=Group(u^-1,v^-1,j^-1,b*u^-1*b^-1);
ix12a:=Group(u^-1,v^-1,j^-1,b*v^-1*b^-1);
ix24a:=Group(u^-1,v^-1,j^-1,b*v^2*b^-1);
ix6:=Group(u^-1,v^-1,b*j);
ix3:=Group(u^-1,v^-1,b*j,j*b);
ix12b:=Group(u^-1,v^-1,b*u^-1*b^-1);
ix9:=Group(u^-1,j^-1,v*u^2*v^-1,b^-1*u^-1*b^-1);
ix18:=Group(u^-1,j^-1,b^-1*u^-1*b^-1,v*b*v^-1*b^-1);
ix16:=Group(v^-1,b^-1,u*v^-1*u^-1);
ix24b:=Group(v^-1,b^-1,u*b^-1*u^-1,u^-1*b^-1*u,u^2*v^-1*u^-2);
ix8:=Group(u*v^-1,j^-1,u^2,u*b*u^-1*b^-1);
ix24c:=Group(u*v^-1,j^-1,u^2,u*b*v^-1*b^-1,b*u^2*b^-1);
ix24d:=Group(u*v^-1,u^2,u*b*u^-1*b^-1);
ix24e:=Group(u*v^-1,u^2,u*b*v^-1*b^-1,b*u^2*b^-1,b*j*u^-1*b);

b^-1*v*b^-1*v^-1*b^-1 in ix24a;

Index(GammaBarFP,ix1) = 1;
Index(GammaBarFP,ix4) = 4;
Index(GammaBarFP,ix12a) = 12;
Index(GammaBarFP,ix24a) = 24;
Index(GammaBarFP,ix6) = 6;
Index(GammaBarFP,ix3) = 3;
Index(GammaBarFP,ix12b) = 12;
Index(GammaBarFP,ix9) = 9;
Index(GammaBarFP,ix18) = 18;
Index(GammaBarFP,ix16) = 16;
Index(GammaBarFP,ix24b) = 24;
Index(GammaBarFP,ix8) = 8;
Index(GammaBarFP,ix24c) = 24;
Index(GammaBarFP,ix24d) = 24;
Index(GammaBarFP,ix24e) = 24;

AbelianInvariants(ix1) = [3]; 
AbelianInvariants(ix4) = [2,3]; 
AbelianInvariants(ix12a) = [2,2,3]; 
AbelianInvariants(ix24a) = [3,4]; 
AbelianInvariants(ix6) = [4]; 
AbelianInvariants(ix3) = [2,2]; 
AbelianInvariants(ix12b) = [2]; 
AbelianInvariants(ix9) = [2,3]; 
AbelianInvariants(ix18) = [3]; 
AbelianInvariants(ix16) = [2,3]; 
AbelianInvariants(ix24b) = [2,2,3]; 
AbelianInvariants(ix8) = [3,3]; 
AbelianInvariants(ix24c) = [2,2,3,3]; 
AbelianInvariants(ix24d) = [3,3]; 
AbelianInvariants(ix24e) = [3,3];