/* ******************************************************************************
   ******************************* PART I ****************************************
   ************************* General Functions  **********************************
   ******************************************************************************* */


// The function "HurwitzMove"  performs the i-th braid move 

BraidMove:=function(V,i)
	c:=V[i]*V[i+1]*V[i]^-1;
return Insert(Remove(V,i+1),i,c); end function;

/* The function "BraidOrbit" determines the orbit of a given generating 
   triple "V" under the braid group. */

BraidOrbit:=function(V)
orb:={ }; sortorb:={ }; Trash:={V}; 
repeat
    ExtractRep(~Trash,~gens); Include(~orb, gens);
		for i in [1..#V-1] do newgens:=BraidMove(gens,i);
			if newgens notin orb then 
				Include(~Trash, newgens);
			end if; 
		end for;
until IsEmpty(Trash);

for gens in orb do  test:=true;
for k in [1..#V-1] do 
    if Order(gens[k]) gt Order(gens[k+1]) then 
		test:=false; break k; 
	end if; 
end for;
    if test then 
		Include(~sortorb, gens); 
	end if;
end for;
return sortorb; end function;

/* The function "GTrps" determines in a fist step the set of all generating triples 
   of signature [3,3,3] for G=He(3) or (Z_3)^2. In the second step, it partitions this 
   set into braid orbits. A representative for each orbit is returned. */ 

GTrps:=function(G)
n:=#G;
Triples:={};
El:=Set(G);
El:=Exclude(El,Id(G));
    for g1 in El do 
        for g2 in El do 
            g3:=(g1*g2)^(-1);
                if g3 ne Id(G) and n eq #sub<G|g1,g2> then 
                    Triples:=Include(Triples,[g1,g2,g3]);
                end if;
        end for;
    end for;
Repres:={};
	while not IsEmpty(Triples) do 
		V:=Rep(Triples); 
		Include(~Repres,V); 
		orb:=BraidOrbit(V);
			for W in orb do 
				Exclude(~Triples, W);
			end for; 
	end while; 
return Repres;
end function;

/* The function "CharCan" returns the character of the canonical representation 
   of the action which corresponds to the given generating triple "V" for G. 
   It also needs the character table "CT" of the group G as an input. */

CharCan:=function(V,CT)
F:=CyclotomicField(3);
ze:=F.1^2;
g1:=V[1]; g2:=V[2];
chi:=CT[1];
for c in CT do 
	if Degree(c) eq 1 then 
		if ze eq c(g1) and ze eq c(g2) then 
			chi:=c; break c;
		end if;
	end if;
end for;
return chi;
end function;


/* The function "ConjTup" performs the simultaneous conjugation of a multiset of 
   generating triples by an automorphism of the group. */

ConjTup:=function(phi,Tup);
NewTup:={* *};
for V in Tup do
	Vnew:=[];
		for g in V do 
			Append(~Vnew,phi(g)); 
		end for;
	Include(~NewTup,Vnew);
end for;
return NewTup;
end function;

// The function "SetToSeq" turns a set or multiset into a sequence

SetToSeq:=function(TSet)
L:=[];
for V in TSet  do
L:=Append(L,V);
end for;
return L;
end function;

/* ******************************************************************************
   ******************************* PART II ***************************************
   ******************Classification of rigid threefolds ************************** 
   ******************* with canonical singularities  ***************************** 
   ******************************************************************************* */

/* The function "RigidThreefolds" determines the set of triples of generating triples 
   [V1,V2,V3], which are compatible in the sense that they correspond to a rigid threefold. 
   Recall that rigidity means that there are no invariant quadratic differentials
   or equivalently chi_{E_i}*chi_{E_j} is not trivial. 
   Any rigid threefold obtained by a diagonal action of "G" is represented by such a triple. 
   However many triples still correspond to isomorphic threefolds, because we have not yet
   considered the full action of "Br3^3xAut(G)xSym(3)".  */ 

RigidThreefolds:=function(G)
TriplesOfGTriples:={};
CT:=CharacterTable(G);
GenerTriples:=GTrps(G);
	for V1 in GenerTriples do
	for V2 in GenerTriples do 
	for V3 in GenerTriples do 
		XE1:=CharCan(V1,CT); 
		XE2:=CharCan(V2,CT); 
		XE3:=CharCan(V3,CT);
		Chi:=XE1*XE2 + XE1*XE3 + XE2*XE3;
			if InnerProduct(Chi,PrincipalCharacter(G)) eq 0 then 
				TriplesOfGTriples:=Include(TriplesOfGTriples,{* V1,V2,V3 *});
			end if;
	end for; end for; end for; 
return TriplesOfGTriples;
end function;

// ************** MAIN ROUTINE FOR THE CLASSIFICATION OF THREEFOLDS *********************

/* The function "FullActionThreefolds" determines the orbits of the action of 
   "Br3^3xAut(G)xSym(3)" on the set of triples of generating triples computed with
   "RigidThreefolds". */

FullActionThreefolds:=function(G)
Trips:=RigidThreefolds(G); Repres:={};
Aut:=AutomorphismGroup(G); 
f,P:= PermutationRepresentation(Aut);

while not IsEmpty(Trips) do 
	Trp:=Rep(Trips);
	TrpS:=SetToSeq(Trp); 
	V1:=TrpS[1];V2:=TrpS[2];V3:=TrpS[3];
	Include(~Repres,Trp);
	orb1:=BraidOrbit(V1); orb2:=BraidOrbit(V2); orb3:=BraidOrbit(V3);
		for W1 in orb1 do 
			for W2 in orb2 do 
				for W3 in orb3 do 
					for p in P do
						phi:=p @@ f; 
						NewTrp:={* W1,W2,W3 *};
						Exclude(~Trips, ConjTup(phi,NewTrp)); 
					end for; 
				end for; 
			end for; 
		end for; 
end while; 
return Repres; end function;


/* The function "ClassifyThreefolds" computes number of singularities "Ngor" and "Nter", 
   the geometric genus "pg" and the Betti numbers "b3" and "b2" of the threefolds determined 
   with "FullActionThreefolds2". It creates a file "RigidThreefolds.txt", which contains 
   this data and generating triples for each threefold. */ 

ClassifyThreefolds:=function(G)
CT:=CharacterTable(G); p:=PrincipalCharacter(G); Ngor:=0; Nter:=0;
F:="RigidThreefolds.txt";
fprintf F, "Rigid product quotient threefolds with group G=%o: \n\n\n", GroupName(G);
count:=1; 

Trips:=FullActionThreefolds(G);
for Trp in Trips do
	TrpS:=SetToSeq(Trp);
	V1:=TrpS[1];V2:=TrpS[2];V3:=TrpS[3];
	XE1:=CharCan(V1,CT); XE2:=CharCan(V2,CT); XE3:=CharCan(V3,CT);
	XE1bar:=ComplexConjugate(XE1); XE2bar:=ComplexConjugate(XE2); XE3bar:=ComplexConjugate(XE3);
	chi11:=XE1bar*XE2 +  XE1*XE2bar +  XE1bar*XE3 +  XE1*XE3bar + XE2bar*XE3 + XE2*XE3bar + 3*p;
	chi21:=XE1bar*XE2*XE3 +  XE1*XE2bar*XE3  + XE1*XE2*XE3bar + 2*(XE1 + XE2 + XE3); 
	pg:=InnerProduct(XE1*XE2*XE3,p);
		if pg eq 1 then 
			b2:=InnerProduct(chi11,p);
			b3:=2;
			Ngor:=3*b2; Nter:=0;
		else 
			Nter:=9;
			h11:=InnerProduct(chi11,p); b2:=h11;
			h12:=InnerProduct(chi21,p); b3:=2*h12;
			Ngor:=3*(h11-h12-2);
		end if;
		fprintf F, "%o)",count;
		fprintf F, " Singularities: Ngor=%o,  Nter=%o \n ", Ngor, Nter;
		fprintf F, "\n Cohomology: pg=%o, b3=%o, b2=%o \n", pg, b3, b2;	
		fprintf F, "\n Generating triples: \n \n V1=%o \n V2=%o \n V3=%o\n", V1, V2, V3;	
		fprintf F, "\n\n\n";
		count:=count+1;
end for;
return "Classification_accomplished"; end function;

/* *************** Run the following function for the groups:  *****************
   		   SmallGroup(9,2)=(Z_3)^2 and SmallGroup(27,3)=He(3)
   ********************************************************************** */


// G1:=SmallGroup(9,2);
// ClassifyThreefolds(G1);

// G2:=SmallGroup(27,3); 
// ClassifyThreefolds(G2);


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


/* ******************************************************************************
   ******************************* PART III ************************************* 
   ******************Classification of smooth rigid fourfolds *******************  
   ****************************************************************************** 
   ****************************************************************************** */

// The function "stab" determines the stabilizer set of a generating triple 

stab:= function(V,G)
M:= {} ;
	for g in V do
		for n in [1 .. Order(g)] do
			M := M join Conjugates(G,g^n) ;
		end for;
	end for;
return M;
end function ;

/* In analogy to the function "RigidThreefolds" we determine the set of quadruples of 
   of generating triples [V1,V2,V3,V4], which are compatible in the sense that they 
   correspond to a rigid fourfold obtained by a free action. 
   Again many of them correspond to isomorphic fourfolds, because we still need to 
   consider the full action of "Br3^4xAut(G)xSym(4)". */ 

RigidSmoothFourFolds:=function(G)
QuadOfGTriples:={};
CT:=CharacterTable(G);
GenerTriples:=GTrps(G);
	for V1 in GenerTriples do
	for V2 in GenerTriples do 
	for V3 in GenerTriples do 
	for V4 in GenerTriples do 
	if #(stab(V1,G) meet stab(V2,G) meet stab(V3,G) meet stab(V4,G)) eq 1 then 
		XE1:=CharCan(V1,CT); XE2:=CharCan(V2,CT); 
		XE3:=CharCan(V3,CT); XE4:=CharCan(V4,CT);		
		Chi:=XE1*XE2 + XE1*XE3 + XE1*XE4 + XE2*XE3 + XE2*XE4 + XE3*XE4;
			if InnerProduct(Chi,PrincipalCharacter(G)) eq 0 then 
				QuadOfGTriples:=Include(QuadOfGTriples,{* V1,V2,V3,V4 *});
			end if;
	end if;
	end for; end for; end for; end for; 
return QuadOfGTriples;
end function;

/* The function "FullActionFourfolds" determines the orbits of the action of 
   "Br3^4xAut(G)xSym(4)" on the set of quadruples of generating triples computed 
   with "RigidSmoothFourFolds". */

FullActionFourfolds:=function(G)

Quads:=RigidSmoothFourFolds(G); Repres:={};
Aut:=AutomorphismGroup(G); 
f,P:= PermutationRepresentation(Aut);

while not IsEmpty(Quads) do 
	Quad:=Rep(Quads);
	QuadSeq:=SetToSeq(Quad); 
	V1:=QuadSeq[1];V2:=QuadSeq[2];V3:=QuadSeq[3];V4:=QuadSeq[4];
	Include(~Repres,Quad);
	orb1:=BraidOrbit(V1); orb2:=BraidOrbit(V2);orb3:=BraidOrbit(V3); orb4:=BraidOrbit(V4);
		for W1 in orb1 do 
			for W2 in orb2 do 
				for W3 in orb3 do 
					for W4 in orb4 do 
						for p in P do
							phi:=p @@ f; 
							NewQuad:={* W1,W2,W3,W4 *};
							Exclude(~Quads, ConjTup(phi,NewQuad)); 
		end for; end for; end for; end for; end for;
end while; 
return Repres; end function;


// G1:=SmallGroup(9,2);
// FullActionFourfolds(G1);

// G2:=SmallGroup(27,3); 
// FullActionFourfolds(G2);

/* ****************** The output is the following: *******************************
   A single orbit for each exceptional group represented by:

> G1:=SmallGroup(9,2);
> FullActionFourfolds(G1);
{   
    {*  
        [ G1.1^2, G1.1^2 * G1.2, G1.1^2 * G1.2^2 ],
        [ G1.1^2, G1.1 * G1.2^2, G1.2 ],
        [ G1.1^2 * G1.2^2, G1.1, G1.2 ],
        [ G1.2, G1.1^2 * G1.2, G1.1 * G1.2 ]
    *}
}
>
> G2:=SmallGroup(27,3);
> FullActionFourfolds(G2);
{   
    {*  
        [ G2.1^2 * G2.2, G2.1 * G2.3^2, G2.2^2 ],
        [ G2.1^2 * G2.2^2 * G2.3^2, G2.2, G2.1 * G2.3 ],
        [ G2.1^2 * G2.2^2 * G2.3, G2.2^2 * G2.3^2, G2.1 * G2.2^2 * G2.3^2 ],
        [ G2.1 * G2.2 * G2.3, G2.1 * G2.2^2 * G2.3, G2.1 ]
    *}
}

*/