#####
#
# This is GAP code
# 
# The file is ~steger/donald/SU21/a2p3/a2p3-3-proj.gap
#
# See file ~steger/donald/SU21/a2p3/D.txt for the definition of the
# division algebra $\D$ and for all the notation.
#
#####

Read("/home/steger/donald/SU21/a2p3/a2p3-proj.gap");
Read("/home/steger/donald/SU21/NormElem.gap");

#####
#
# Free group
#
######

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

RelnList:=[
(Af*Bf^-1)^3,
(Bf*Cf*Af*Cf*Af^-2)^3,
(Cf^-1*Bf^-1*Cf^-1*Af^-1)^3,
Af*Bf^-1*Cf^-1*Bf*Cf*Af*Bf*Cf^-1*Bf^-1*Af*Cf*Af,
Bf*Cf^-1*Af^-1*Cf^-1*Bf^-1*Af*Cf^-1*Af^-1*Cf^-1*Bf^-1*Cf^-1*Af^-1*Bf*Cf,
Cf*Bf^2*Cf*Af*Cf^-1*Bf^-1*Cf^-1*Bf^-2*Cf^-1*Bf*Cf*Af,
Af^-1*Bf^-1*Cf^-1*Bf^-1*Cf^-1*Bf^-2*Cf^-1*Af^-2*Cf^-1*Bf^2*Cf*Bf*Cf*Bf*Af^-1,
Bf^-1*Cf^-1*Bf^-3*Af^-1*Cf^-1*Bf^-2*Cf^-1*Bf^-1*Cf^-1*Af^-1*Cf^-1*Af^-1*Bf*Af^-1*Cf^-1,
Af^-2*Cf^-1*Bf^-1*Cf*Bf^2*Cf*Af*Bf^2*Cf*Bf^2*Cf*Af*Bf*Cf^-1,
Bf*Af^-1*Cf*Bf^2*Cf*Af^3*Bf^-1*Cf^-1*Bf^-1*Cf^-1*Bf^-2*Cf*Bf*Cf,
Bf^-1*Cf^-1*Bf^-1*Cf^-1*Af*Bf^-1*Cf^-1*Bf^-1*Cf^-1*Bf*Cf*Bf*Cf*Bf*Af^-1*Cf*Bf^-1*Af^-1*Cf^-1*Bf^-1,
Cf*Bf^2*Cf*Af*Bf*Cf^-1*Af^-1*Cf^-1*Bf^-1*Af^3*Bf^-1*Cf^-1*Bf^-1*Cf^-1*Bf^-3,
Af*Bf*Cf*Af*Bf^-1*Cf^-1*Af^-1*Cf^-1*Bf^-1*Af^-1*Cf^-1*Af^-1*Cf^-1*Bf*Cf*Af*Bf*Cf^-1*Af^-1*Cf^-1*Bf^-1*Af,
Af^-1*Bf*Cf*Af^-1*Bf*Cf*Af^-1*Cf^-1*Bf^-3*Cf^-1*Af^-2*Cf^-1*Af^-1*Bf*Af^-1*Cf^-1*Bf^-1*Cf^-1*Af^-1*Bf*Cf,
Af*Cf*Bf^-1*Cf*Af*Bf*Cf*Af*Cf*Bf^-1*Cf^-1*Af^-1*Cf*Bf^-1*Cf^-1*Bf^-1*Cf^-1*Bf*Cf*Af*Bf*Af^-1*Bf*Cf
];

GammaBarFP:=FreeGp/RelnList;;
AFP:=GammaBarFP.1;
BFP:=GammaBarFP.2;
CFP:=GammaBarFP.3;

index3aFP:=Group(AFP,CFP,BFP*AFP*BFP^-1,BFP^3,BFP*CFP*BFP^-1);
index3bFP:=Group(AFP,CFP,BFP^2,BFP*CFP*BFP^-1);

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

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

# Check the indexes of the subgroups

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

# Check normality

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

# Free group versions of these subgroups
FPtoFGSubgp:=g->Group(List(GeneratorsOfGroup(g),UnderlyingElement));

index3aFG:=FPtoFGSubgp(index3aFP);
index3bFG:=FPtoFGSubgp(index3bFP);

#####
#
# Generators of $(a=2,p=3,\{3\})$
#
#####

#From: "Donald Cartwright" <donald.cartwright@gmail.com>
#To: steger@uniss.it
#Subject: (a=2,p=3) results [Re: (a=1,p=5) ... hmmm]
#Date: Mon, 24 Dec 2007 16:41:23 +1100

ADM:=One(mm)*[
[1/10*(-2*Sqrt(-2)-1)*Zm^2+1/10*(-Sqrt(-2)-7)*Zm+1/10*Sqrt(-2),
1/20*(-Sqrt(-2)+2)*Zm^2+1/20*(3*Sqrt(-2)-4)*Zm-3/10,
1/5*Sqrt(-2)*Zm^2+1/5*Sqrt(-2)*Zm+3/5*Sqrt(-2)],
[1/10*(-3*Sqrt(-2)+6)*Zm^2-3/10*Sqrt(-2)*Zm+1/5*(3*Sqrt(-2)+3),
-1/5*Zm^2+1/5*(-2*Sqrt(-2)+1)*Zm+1/5*(Sqrt(-2)+1),
1/10*(-Sqrt(-2)-2)*Zm^2+1/10*(-3*Sqrt(-2)-2)*Zm+1/5*(-Sqrt(-2)-3)],
[1/20*(9*Sqrt(-2)+6)*Zm^2+1/20*(9*Sqrt(-2)+24)*Zm+1/10*(6*Sqrt(-2)-3),
1/10*(-3*Sqrt(-2)-6)*Zm^2-3/10*Sqrt(-2)*Zm+1/5*(3*Sqrt(-2)-3),
1/10*(2*Sqrt(-2)+3)*Zm^2+1/2*(Sqrt(-2)+1)*Zm+1/10*(7*Sqrt(-2)+8)]];

BDM:=One(mm)*[
[3/10*Zm^2+1/10*(-Sqrt(-2)+1)*Zm+3/10*Sqrt(-2),
1/20*Zm^2+1/20*(-Sqrt(-2)-5)*Zm+1/20*(3*Sqrt(-2)+6),
1/30*(4*Sqrt(-2)-1)*Zm^2+1/30*(5*Sqrt(-2)+1)*Zm+1/10*(3*Sqrt(-2)+2)],
[1/10*(-Sqrt(-2)+4)*Zm^2+1/10*(-Sqrt(-2)-2)*Zm+1/5*(Sqrt(-2)+5),
1/10*(-2*Sqrt(-2)-1)*Zm^2+1/10*(-Sqrt(-2)+1)*Zm+1/10*(Sqrt(-2)-4),
-1/10*Sqrt(-2)*Zm^2+1/5*(-Sqrt(-2)+1)*Zm+2/5],
[1/20*(6*Sqrt(-2)-3)*Zm^2+1/20*(15*Sqrt(-2)+3)*Zm+1/20*(21*Sqrt(-2)+18),
1/10*(-3*Sqrt(-2)-3)*Zm^2+1/10*(-4*Sqrt(-2)+1)*Zm+1/10*(-Sqrt(-2)+4),
1/5*(Sqrt(-2)-1)*Zm^2+1/5*(Sqrt(-2)-1)*Zm+1/5*(3*Sqrt(-2)-3)]];

CDM:=One(mm)*[
[1/2*Sqrt(-2)*Zm^2+1/2*(Sqrt(-2)-2)*Zm,
0,
0],
[0,
1/2*(-Sqrt(-2)+1)*Zm^2+1/2*(-2*Sqrt(-2)+1)*Zm+1/2*(-3*Sqrt(-2)+2),
0],
[0,
0,
-1/2*Zm^2+1/2*(Sqrt(-2)+1)*Zm-1/2*s]];

ADMVec:=DMtoDMVec(ADM);;
BDMVec:=DMtoDMVec(BDM);;
CDMVec:=DMtoDMVec(CDM);;

ADMVec in DMVec;
BDMVec in DMVec;
CDMVec in DMVec;

IsOne(iotaDM(ADM)*ADM);
IsOne(iotaDM(BDM)*BDM);
IsOne(iotaDM(CDM)*CDM);

# Check that the relations hold (projectively) for all our matrix
# generators

GammaBarDM:=Group(ADM,BDM,CDM);;
DMHom:=GroupHomomorphismByImages(FreeGp,GammaBarDM,
  [Af,Bf,Cf],[ADM,BDM,CDM]);;
ForAll(RelnList,x->IsPOne(x^DMHom));
	
#####
#
# Construct matrices for the group generators acting on~$\D$ by
# conjugation
#
######

DMtoPDM:=mtx->TransposedMat(
  List(basisDM,
    bmtx->Coefficients(basisDMVec,DMtoDMVec(mtx*bmtx*mtx^-1))
));

APDM:=DMtoPDM(ADM);;
BPDM:=DMtoPDM(BDM);;
CPDM:=DMtoPDM(CDM);;

IsZero(LieBracket(iotaPDM,APDM));
IsZero(LieBracket(iotaPDM,BPDM));
IsZero(LieBracket(iotaPDM,CPDM));

GammaBarPDM:=Group(APDM,BPDM,CPDM);;
FGtoPDM:=GroupHomomorphismByImages(FreeGp,GammaBarPDM,
  [Af,Bf,Cf],[APDM,BPDM,CPDM]);;
ForAll(RelnList,x->IsOne(x^FGtoPDM));

######
#
# Find the matrices for the group generators acting by conjugation on
# the self-adjoint, traceless part of~$\D$.
#
#####

PDMtoPSA:=mtx->TransposedMat(
  List(DMVecSelfAdjBasis,SAVec->
    Coefficients(DMVecSelfAdjBasis,mtx*SAVec)
));

APSA0:=PDMtoPSA(APDM);;
BPSA0:=PDMtoPSA(BPDM);;
CPSA0:=PDMtoPSA(CPDM);;

#####
#
# Adjust the basis of the self-adjoint part of~$\D$ so that
# the group generators correspond to integral matrices.
#
#####

GammaBarPSA0:=Group(APSA0,BPSA0,CPSA0);;

ILatPSA0:=InvariantLattice(GammaBarPSA0)^-1;
APSA:=APSA0^ILatPSA0;
BPSA:=BPSA0^ILatPSA0;
CPSA:=CPSA0^ILatPSA0;

GammaBarPSA:=Group(APSA,BPSA,CPSA);;
FGtoPSA:=GroupHomomorphismByImages(FreeGp,GammaBarPSA,
  [Af,Bf,Cf],[APSA,BPSA,CPSA]);;
ForAll(RelnList,x->IsOne(x^FGtoPSA));

######
#
# Do calculations $\mod 2$
#
######

toMod2:=x->x*One(GF(2));
toMod2M:=mtx->fOnM(toMod2,mtx);

APSAMod2:=toMod2M(APSA);;
BPSAMod2:=toMod2M(BPSA);;
CPSAMod2:=toMod2M(CPSA);;

GammaBarPSAMod2:=Group(APSAMod2,BPSAMod2,CPSAMod2);;
FGtoPSAMod2:=GroupHomomorphismByImages(FreeGp,GammaBarPSAMod2,
  [Af,Bf,Cf],[APSAMod2,BPSAMod2,CPSAMod2]);;
ForAll(RelnList,x->IsOne(x^FGtoPSAMod2));
	
Size(GammaBarPSAMod2) = 192;
#StructureDescription(GammaBarPSAMod2) = "(((C4 x C4) : C3) : C2) : C2";

index3aPSAMod2:=Image(FGtoPSAMod2,index3aFG);
index3bPSAMod2:=Image(FGtoPSAMod2,index3bFG);

index3aPSAMod2 = GammaBarPSAMod2;
Index(GammaBarPSAMod2,index3bPSAMod2) = 3;

######
#
# Do calculations $\mod 3$
#
######

toMod3:=x->x*One(GF(3));
toMod3M:=mtx->fOnM(toMod3,mtx);

APSAMod3:=toMod3M(APSA);;
BPSAMod3:=toMod3M(BPSA);;
CPSAMod3:=toMod3M(CPSA);;

GammaBarPSAMod3:=Group(APSAMod3,BPSAMod3,CPSAMod3);;
FGtoPSAMod3:=GroupHomomorphismByImages(FreeGp,GammaBarPSAMod3,
  [Af,Bf,Cf],[APSAMod3,BPSAMod3,CPSAMod3]);;
ForAll(RelnList,x->IsOne(x^FGtoPSAMod3));
	
Size(GammaBarPSAMod3) = 3*13*3^6;
ChSGammaBarPSAMod3:=ChiefSeries(GammaBarPSAMod3);
Length(ChSGammaBarPSAMod3) = 5;
List(ChSGammaBarPSAMod3,Size)=[3*13*3^6,13*3^6,3^6,3^3,1];
Mod3H1:=ChSGammaBarPSAMod3[3];
Mod3H0:=ChSGammaBarPSAMod3[4];
Mod3H0 = Center(Mod3H1);
StructureDescription(Mod3H0) = "C3 x C3 x C3";
StructureDescription(Mod3H1/Mod3H0) = "C3 x C3 x C3";
#
# whence follows
# StructureDescription(Mod3H1) = "(C3 x C3 x C3) . (C3 x C3 x C3)"
#
Mod3Q1Hom:=NaturalHomomorphismByNormalSubgroup(GammaBarPSAMod3,Mod3H1);;
Mod3Q1:=Image(Mod3Q1Hom);
StructureDescription(Mod3Q1) = "C13 : C3";

index3aPSAMod3:=Image(FGtoPSAMod3,index3aFG);
index3bPSAMod3:=Image(FGtoPSAMod3,index3bFG);

Index(GammaBarPSAMod3,index3aPSAMod3) = 3;
index3bPSAMod3 = GammaBarPSAMod3;

FGtoMod3Q1:=CompositionMapping(Mod3Q1Hom,FGtoPSAMod3);
ForAll(RelnList,x->IsOne(x^FGtoMod3Q1));

index3aMod3Q1:=Image(FGtoMod3Q1,index3aFG);
index3bMod3Q1:=Image(FGtoMod3Q1,index3bFG);

Index(Mod3Q1,index3aMod3Q1) = 3;
index3bMod3Q1 = Mod3Q1;

######
#
# Do calculations $\mod 5$
#
######

toMod5:=x->x*One(GF(5));
toMod5M:=mtx->fOnM(toMod5,mtx);

APSAMod5:=toMod5M(APSA);;
BPSAMod5:=toMod5M(BPSA);;
CPSAMod5:=toMod5M(CPSA);;

GammaBarPSAMod5:=Group(APSAMod5,BPSAMod5,CPSAMod5);;
FGtoPSAMod5:=GroupHomomorphismByImages(FreeGp,GammaBarPSAMod5,
  [Af,Bf,Cf],[APSAMod5,BPSAMod5,CPSAMod5]);;
ForAll(RelnList,x->IsOne(x^FGtoPSAMod5));

ChSGammaBarPSAMod5:=ChiefSeries(GammaBarPSAMod5);
Length(ChSGammaBarPSAMod5) = 3;
Mod5H1:=ChSGammaBarPSAMod5[2];
IsSimple(Mod5H1);
#IsomorphismGroups(PSU(3,5),Mod5H1);
#IsomorphismGroups(PGU(3,5),GammaBarPSAMod5);

index3aPSAMod5:=Image(FGtoPSAMod5,index3aFG);
index3bPSAMod5:=Image(FGtoPSAMod5,index3bFG);

Index(GammaBarPSAMod5,index3aPSAMod5) = 3;
index3bPSAMod5 = GammaBarPSAMod5;

######
#
# Do calculations $\mod 7$
#
######

toMod7:=x->x*One(GF(7));
toMod7M:=mtx->fOnM(toMod7,mtx);

APSAMod7:=toMod7M(APSA);;
BPSAMod7:=toMod7M(BPSA);;
CPSAMod7:=toMod7M(CPSA);;

GammaBarPSAMod7:=Group(APSAMod7,BPSAMod7,CPSAMod7);;
FGtoPSAMod7:=GroupHomomorphismByImages(FreeGp,GammaBarPSAMod7,
  [Af,Bf,Cf],[APSAMod7,BPSAMod7,CPSAMod7]);;
ForAll(RelnList,x->IsOne(x^FGtoPSAMod7));
	
#IsomorphismGroups(PGU(3,7),GammaBarPSAMod7);
#Size(GammaBarPSAMod7) = 5663616;

index3aPSAMod7:=Image(FGtoPSAMod7,index3aFG);
index3bPSAMod7:=Image(FGtoPSAMod7,index3bFG);

index3aPSAMod7 = GammaBarPSAMod7;
index3bPSAMod7 = GammaBarPSAMod7;

######
#
# Do calculations $\mod 11$
#
######

toMod11:=x->x*One(GF(11));
toMod11M:=mtx->fOnM(toMod11,mtx);

APSAMod11:=toMod11M(APSA);;
BPSAMod11:=toMod11M(BPSA);;
CPSAMod11:=toMod11M(CPSA);;

GammaBarPSAMod11:=Group(APSAMod11,BPSAMod11,CPSAMod11);;
FGtoPSAMod11:=GroupHomomorphismByImages(FreeGp,GammaBarPSAMod11,
  [Af,Bf,Cf],[APSAMod11,BPSAMod11,CPSAMod11]);;
ForAll(RelnList,x->IsOne(x^FGtoPSAMod11));
	
#IsSimple(GammaBarPSAMod11);
#Size(GammaBarPSAMod11) = 212427600;
#StructureDescription(GammaBarPSAMod11) = "PSL(3,11)";

index3aPSAMod11:=Image(FGtoPSAMod11,index3aFG);
index3bPSAMod11:=Image(FGtoPSAMod11,index3bFG);

index3aPSAMod11 = GammaBarPSAMod11;
index3bPSAMod11 = GammaBarPSAMod11;

#####
#
# Consider abelianizations
#
#####

# (a=2,p=3,\{2\})
AbelianInvariants(GammaBarFP) = [2,2,3];
AbelianInvariants(GammaBarPSAMod2) = [2,2];
AbelianInvariants(GammaBarPSAMod3) = [3];
AbelianInvariants(GammaBarPSAMod5) = [3];

# (a=2,p=3,\{2\},D_3)
AbelianInvariants(index3aFP) = [2,2,13];
AbelianInvariants(index3aPSAMod2) = [2,2];
AbelianInvariants(index3aPSAMod3) = [13];
AbelianInvariants(index3aPSAMod5) = [];

# (a=2,p=3,\{2I\})
AbelianInvariants(index3bFP) = [2,2,2,2,3];
AbelianInvariants(index3bPSAMod2) = [2,2,2,2];
AbelianInvariants(index3bPSAMod3) = [3];
AbelianInvariants(index3bPSAMod5) = [3];