// First we compute the set of Types from Lemma 2.2. 
// The function "Theta" computes for a given type "T" the number 
// "Theta(T)" (cf. Definition 3.1). 
// The function "TestCond" is a boolean function, which decides 
// if a type T satisfies the conditions given in Definition 3.1. 
// Note that by lemma 2.2 the length of the signatures is bounden by 8.


Theta:= function(T)
b:= -2;
for m in T do 
b:= b + 1 - 1/m ;
end for ;
return b ;
end function ;


TestCond:=function(T)
test:=true;
a:=Theta(T) ;
if a gt 0 then 
	for n in T do 
		if 4/(n*a) notin IntegerRing()
		then 
		test:=false;
		break n ;
		end if;
	end for; 
else test:=false;
end if;
return test;
end function;


ShortTupel:=function(T)
S:=[];
for i in [1..#T] do 
if T[i] ne 1 then 
S:=Append(S,T[i]);
end if;
end for;
return S;
end function;


ListTupels:= [] ;
for m1 in [1..10]  do for m2 in [m1..10] do 
for m3 in [m2..10] do for m4 in [m3..10] do
for m5 in [m4..10] do for m6 in [m5..30] do
for m7 in [m6..30] do for m8 in [m7..30] do
T:=[Integers()| m1,m2,m3,m4,m5,m6,m7,m8] ;
a:=Theta(T) ;
if 
TestCond(T) 
then
ListTupels:=Append(ListTupels,ShortTupel(T)) ;
end if;
end for; end for; end for ; end for;
end for ; end for; end for ;
end for ;


// The function "Breuer" checks if a group of Order "GOrd" can act 
// on a compact Riemann surface C of genus "g". If "GOrd" is greater 
// than the maximum allowed group order given in Breuers table "false" 
// is returned, else "true" is returned. This works if 2 <= g <= 48, 
// if g >= 49 the function returns "true". 


Breuer:=function(GOrd,g)
test:=true;
Breuer_Table:=
[1,48,168,120,192,150,504,336,320,432,240,120,360,1092,
504,720,1344,168,720,228,480,1008,192,216,720,750,
624,1296,672,264,720,372,768,1320,544,672,1728,444,912,
936,960,410,1512,516,1320,2160,384,408];
if g le 48 then 
	if GOrd gt Breuer_Table[g] then 
	test:=false;
	end if;
end if;
return test;
end function;


// The function  "TripleTypes_Gord"  determines for a given group order "GOrd" the 
// set of triples of the form "(GOrd,T_1,T_2)", such that  "GOrd = alpha(T_1)*alpha(T_2)/2".


TripleTypes_Gord:=function(GOrd)
Set:={ };
L2:=ListTupels;
for T1 in ListTupels do 
for T2 in L2 do 
alph_1:=4/Theta(T1);
alph_2:=4/Theta(T2);
n:= (alph_1*alph_2)/2; 
if n eq GOrd
and Breuer(GOrd,IntegerRing()!(alph_1 + 1)) and Breuer(GOrd,IntegerRing()!(alph_2 + 1)) 
then                         
Include(~Set,[ [GOrd] ,T1,T2 ]);
end if;
end for ;
L2 := Reverse(Prune(Reverse(L2))) ;  
end for; 
return Set;
end function;


// This is the set of exceptional group orders.  


BadOrders:={1920, 1152, 768, 512, 384, 256};


// The function "ElementsOfOrd" determines all elements 
// of order "n" in a given group "G". 


ElementsOfOrd:=func<G, n | {g: g in G| Order(g) eq n}>;


// The function "TuplesOfGivenOrder" creates a sequence of the same 
// length as the input sequence "type", whose entries are subsets of the group 
// "G" in the input, and precisely the subsets of elements of order
// the corresponding entry of "type".


TuplesOfGivenOrders:=function(G,type)
SEQ:=[];
for i in [1..#type] do
	el:=ElementsOfOrd(G,type[i]);
	if IsEmpty(el) then return [];
		else Append(~SEQ,el);
	end if; 
end for;
return SEQ;
end function;


// The function "TupleToSeq" transforms a tuple into a sequence.


TupleToSeq:=function(cands)
T:=[];
for t in cands do 
Append(~T,t);
end for;
return T;
end function;


// The function "ExSphGens" is a boolean function which decides 
// if a group "G" has a system of generators of type "T". 


ExSphGens:=function(G,T)
test:=false;
t:=T[#T];
SetCands:=TuplesOfGivenOrders(G,Remove(T,#T));
if not IsEmpty(SetCands) then
	cands:= CartesianProduct(SetCands);
		for cand in cands do
			m:=[];
			for i in {1..(#T-1)} do m[i]:=cand[i];
			end for;
		if Order(&*m) eq t then
		if #sub<G|TupleToSeq(cand)> eq #G then
			test:=true; break cand;
		end if; 
		end if; 
		end for; 
end if;
return test;
end function;


// The function "stab" computes the stabilizer set of a 
// spherical system of generators "A" of a group "G" according 
// to Definition 3.3. 


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


// To implement the group actions described in section 4 we need the 
// automorphism group of a given group "G" as an explicit set.
// This is done with the function "AutGr". 
// The function "HurwitzMove" creates a single Hurwitz move.
 

AutGr:=function(G)
Aut:=AutomorphismGroup(G); A:={ Aut!1 };
repeat
for g1 in Generators(Aut) do
for g2 in A do
Include (~A,g1*g2);
end for; end for;
until #A eq #Aut; return A; end function;


HurwitzMove:=function(seq,idx)
return Insert(Remove(seq,idx),idx+1,seq[idx]^seq[idx+1]);
end function;


// The function "HurwitzOrbit" creates the Orbit of a given spherical system of 
// generators "seq" under the Hurwitz action. To save memory it returns only 
// the ordered spherical systems of generators (i.e. ord(g_i)=m_i). 


HurwitzOrbit:=function(seq)
orb:={ }; Purgatory:={ seq };
repeat
ExtractRep(~Purgatory,~gens); Include(~orb, gens);
for k in [1..#seq-1] do hurgens:=HurwitzMove(gens,k);
if hurgens notin orb then Include(~Purgatory, hurgens);
end if; end for;
until IsEmpty(Purgatory);
orbcut:={ };
for gens in orb do test:=true;
for k in [1..#seq-1] do
if Order(gens[k]) gt Order(gens[k+1]) then test:=false; break k;
end if;
end for;
if test then Include(~orbcut, gens);
end if;
end for; return orbcut; end function;


// The function "SphGens" creates all the spherical systems of "G" of type "T".  


SphGens:=function(G, T)
Vect:={};
t:=T[#T];
SetCands:=TuplesOfGivenOrders(G,Remove(T,#T));
for cand in CartesianProduct(SetCands) do
    m:=[]; for i in {1..(#T-1)} do m[i]:=cand[i];end for;
	pr:=&*m;
	if Order(pr) eq t then
    if #sub<G|TupleToSeq(cand)> eq #G  then
       Include(~Vect, TupleToSeq(cand) cat [pr^-1]);
	end if; 
	end if; 
end for;
return Vect;
end function;


// The function "B1" determines for each B_r x Aut(G) orbit 
// exactly one representative. 


B1:=function(G,type) ;
W:=SphGens(G,type) ;
Aut:=AutGr(G);
Repres:={ }; 
while not IsEmpty(W) do
Tup:=Rep(W);
Include(~Repres,Tup);
orb1:=HurwitzOrbit(Tup);
for g1 in orb1 do for phi in Aut do
Exclude(~W, phi(g1));
end for; 
end for;
end while;
return Repres;
end function;


// The function "B2" determines for each orbit of the Hurwitz action 
// exactly one representative.


B2:=function(G,type) ;
W:=SphGens(G,type) ;
Repres:={ }; 
while not IsEmpty(W) do
Tup:=Rep(W);
Include(~Repres,Tup);
orb1:=HurwitzOrbit(Tup);
for g1 in orb1 do
Exclude(~W, g1);
end for; 
end while;
return Repres;
end function;


// The function  "FindComp" computes for each pair in  B1(G,T) x B2(G,S) the 
// intersection of the stabilizer sets.  
// If the intersection is trivial, then the pair will be returned.  
// According to Lemma 4.4 we have at least one representative for each B_r x B_s x Aut(G)-orbit. 


FindComp:=function(G,sig1,sig2)
Comps:={ };
for K1 in B1(G,sig1) do 
for K2 in B2(G,sig2) do 
if  IsEmpty(stab(K1,G) meet stab(K2,G)) then
Include(~Comps, [K1,K2]) ;
end if ;
end for ;
end for ;
return [* Comps, sig1, sig2 *] ;
end function ;


// In some cases the function "FindComp" returns more pairs then B_r x B_s x Aut(G)-orbits.
// We have to identify equivalent pairs. This is done as follows: 
// The function "PartOrb"  separates the pairs which are in different B_r x B_s x Aut(G)-orbits 
// using the method from Lemma 4.4. 
// With the function "Orbi" we identify the remaining pairs under the action of B_r x B_s x Aut(G). 
// The function "AlleOrbits" then returns for each B_r x B_s x Aut(G)-orbit exactly one representative. 


Orbi:=function(W,G)
Aut:=AutGr(G);
Repres := {  };
while not IsEmpty(W) do
Tup:=Rep(W);
Include(~Repres,Tup);  
orb1:=HurwitzOrbit(Tup[1]); orb2:=HurwitzOrbit(Tup[2]);
for g1 in orb1 do for g2 in orb2 do for phi in Aut do
DTup:=[phi(g1),phi(g2)];
Exclude(~W, DTup);
end for; end for; end for;
end while;
return Repres;
end function;


PartOrb:=function(U,gens,G )
Set:={};
for T in gens do 
if U eq T[1] then 
Set:=Include(Set,T);
end if;
end for;
if #Set ne 1 then 
Orbsource:=Orbi(Set,G);
else 
Orbsource:= Set;
end if;
return Orbsource;
end function;

AlleOrbits:=function(comps,G)
Orbits:=[ ];
Set:={};
for T in comps do 
Set:=Include(Set,T[1]);
end for;
for S in Set do 
Pa:=PartOrb(S,comps,G);
Orbits:=Append(Orbits,Pa);
end for;
return Orbits;
end function;


// The function "flip" exchanges the tuples "T" and "S". In the case where 
// "T" and "S" are different it computes both "M1:=FindComp(G,S,T)" and "M2:=FindComp(G,T,S)". 
// It compares the cardinality of "M1" and M2" and returns the smaller set. 
// Doing this we can avoid to compute the full  B_r x B_s x Aut(G)-orbit in most cases.
// If not, we can at least speed up the computation, because we have to identify less pairs 
// under the action of B_r x B_s x Aut(G). 


flip:=function(comps,G)
kleincomps:=comps;
U:=comps[2];
V:=comps[3];
if V ne U then
comps2:=FindComp(G,V,U);
if #comps2[1] le #comps[1] then 
kleincomps:=comps2;
end if;
end if;
return kleincomps;
end function;


// With the function "Faelle" we do a case differentiation as described in section 4 after Lemma 4.4. 


Faelle:=function(Comps);
Text:="We have to identify the pairs under the action of  B_r x B_s x Aut(G)."; 
t:=1;
if #Comps eq 0 then 
Text:="This group has no disjoint pair of spherical systems of generators."; 
t:=0;
else 
Set:={};
for T in Comps do 
Set:=Include(Set,T[1]);
end for;
if #Set eq #Comps then 
Text:="These are all orbits of B_r x B_s x Aut(G).";
t:=0;
end if;
end if;
return [* Text,t *];
end function;


// With the function "Perfect_test" we search through the database of perfect 
// groups if there is a perfect group "G" of order "GOrd" admitting a pair "(A1,A2)" 
// of spherical systems of generators of type "(T1,T2)".  


Perfect_test:=function(GOrd,T1,T2)
Text:="There are no perfect groups of this order admitting a disjoint pair of spherical systems of generators
with T_i in N.";
DB:=PerfectGroupDatabase();
nr:=NumberOfGroups(DB,GOrd);
if nr ge 1 then 
	for i in [1..nr] do
f:=PermutationRepresentation(DB, GOrd, i);
G:=Image(f);
	if ExSphGens(G,T1) and ExSphGens(G,T2) then 
	Text:="We have to investigate this case further";
	Text:=[* G, T1, T2, Text *];		
	break i;
	end if;
	end for;
end if;
return Text;
end function;


// The function "Orbifold" determines the polygonal group of type "T".


Orbifold:=function(seq)
F:=FreeGroup(#seq); Rel:={F![1..#seq]};
for i in [1..#seq] do
Include(~Rel,F.i^(seq[i])); end for;
return quo<F|Rel>;
end function;


// The function "AbelGroups" computes all possible groups, which have the correct 
// abelianization w.r.t the given types. With this step we can exclude the 
// groups which don't have the right abelianization. 


SetofTypes:=function(triples)
Set:={};
for T in triples do
Set:=Set join {T[2],T[3]};
end for;
return Set;
end function;


PossQuotients:=function(SetTuples)
Quot:={};
for T in SetTuples do
	group:=AbelianQuotient(Orbifold(T));
	I:=IdentifyGroup(group);
	G:=SmallGroup(I[1],I[2]);
	Li:=LowIndexSubgroups(G,I[1]);
		for H in Li do
		Q:=quo<G|H>;
		Include(~Quot,IdentifyGroup(Q));
		end for;
end for;
return Quot;
end function;


AbelGroups:=function(GOrd, triples)
Quot:=PossQuotients(SetofTypes(triples));
		max:=0;
		for Id in Quot do
		if Id[1] ge max then
		max:= Id[1];
		end if;
		end for;
checked:=[* *];
if not IsEmpty(SetofTypes(triples)) then  
P:= SmallGroupProcess(GOrd);
     repeat
     H := Current(P);  
		AbG:=AbelianQuotient(H);	 
		if  #AbG le max then 
		if  IdentifyGroup(AbelianQuotient(H)) in Quot then
			Append(~checked,H); 
		end if;
		end if;
     Advance(~P);
until IsEmpty(P);
end if;
return checked;
end function;


// The function "Main_loop1" is the main routine to classify all surfaces, 
// such that "|G| ge 2000" or "|G| in {1920, 1152, 768, 512, 384, 256}". 
//
// With  the command "Main_loop1(n)" the script classifies all surfaces 
// where "|G|<=n" and "|G| ge 2000" or "|G|<=n" and "|G| in {1920, 1152, 768, 512, 384, 256}".
// 
// The function "Main_loop2" is the main routine to classify all surfaces, 
// such that "|G| <=2000" and "|G| notin {1920, 1152, 768, 512, 384, 256}". 
//
// With the command "Main_loop2(n1,n2,1)" the script classifies all surfaces 
// where "n1<=|G|<=n2" and "T_i \ne [2^8]". 
//
// With the command "Main_loop2(n1,n2,0)" the script classifies all surfaces 
// where "n1<=|G|<=n2" and "T_1 = [2^8]". 
//
// The function "Main_loop2" calls the function "Classify". This function performs 
// step2-step4  according to section 5. The input is a triple of the form 
// "(GOrd,T1,T2)" (from step 1) (cf. section 5). 


Classify:=function(T)
n:=T[1][1];
if n notin BadOrders then  
			for G in SmallGroups(IntegerRing()!n :Warning:= false) do 
				if ExSphGens(G,T[2]) then 
				if ExSphGens(G,T[3]) then 
					if 	8 in {#T[2],#T[3]} and  #G in {48,96,168}  then				
						Text:="This is an exceptional case.";
						Tuple:=[* IdentifyGroup(G),T[2],T[3], Text *];					
					PrintFile("loop2.txt",Tuple);
					else
					OrbHalf:=FindComp(G,T[2],T[3]);
						if Faelle(OrbHalf[1])[2] eq 1 then
							OrbHalf:=flip(OrbHalf,G);
						end if;
						Er:=[* IdentifyGroup(G),OrbHalf[2],OrbHalf[3],
						OrbHalf[1], Faelle(OrbHalf[1])[1] *] ;
						PrintFile("loop2.txt",Er);
							if Faelle(OrbHalf[1])[2] eq 1 then 
								FullActionOrbits:=AlleOrbits(OrbHalf[1],G);
								PrintFile("loop2.txt",FullActionOrbits);
							end if;
					end if;
				end if;
				end if;
			end for;
end if; 
return "done";
end function;


Main_loop2:=function(n1,n2,switch);
for n in [n1..n2] do 
for T in TripleTypes_Gord(n) do 
	if switch eq 1 then 
		if #T[2] ne 8 and  #T[3] ne 8 then 
			Classify(T);
		end if;
	else
		if 8 in {#T[2],#T[3]} then
		Classify(T);
		end if;
	end if;
end for;
end for;
return "main loop2 done";
end function;


Main_loop1:=function(n);
for GOrd in [1..n] do
	if GOrd ge 2001 then 
		Triples:=TripleTypes_Gord(GOrd);
		if #Triples ne 0 then 
			for T in Triples do 
				Tuple:=[* GOrd, T[2], T[3] *];				
				Text:="The following group order is treated:";
					PrintFile("exceptional.txt",Text);
					PrintFile("exceptional.txt",GOrd);
					Text:=Perfect_test(GOrd,T[2],T[3]);  
					PrintFile("exceptional.txt",Text);
					PrintFile("exceptional.txt","*"); PrintFile("exceptional.txt","*");
					if 4/Theta(T[2]) ne 168 and 4/Theta(T[3]) ne 168 then 
					Text:="We have to search for non-perfect groups in the following case:"; 	
						PrintFile("exceptional.txt",Text);
						PrintFile("exceptional.txt",Tuple);
					PrintFile("exceptional.txt","*"); PrintFile("exceptional.txt","*");	
					end if;
			end for;
		end if;
	else 
		if GOrd in BadOrders then 
			Text:="The following exceptional group order is treated:";			
			PrintFile("exceptional.txt",Text);
			PrintFile("exceptional.txt",GOrd);
			Triples:=TripleTypes_Gord(GOrd);
				Grps:=AbelGroups(GOrd,Triples);			
				Set:={};
				for T in Triples do 
				for G in Grps do
					if ExSphGens(G,T[2]) and ExSphGens(G,T[3]) then				
					if #FindComp(G,T[2],T[3])[1] ge 1  then 
					id:=IdentifyGroup(G);
					Include(~Set,id);
					Tuple:=[* id, T[2], T[3] *];		
							Text:="We have to investigate the following triple:";
							PrintFile("exceptional.txt",Text);
							PrintFile("exceptional.txt",Tuple);
					end if; end if;
				end for;
				end for;
				if IsEmpty(Set) then 
Text:="There are no groups of this order admitting a disjoint pair of spherical systems of generators
with T_i in N.";
					PrintFile("exceptional.txt",Text); PrintFile("exceptional.txt","*"); PrintFile("exceptional.txt","*");
				end if;
		end if;
	end if;
end for;
return "main loop1 done";
end function;


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

// This is the source code to perform the computations in Proposition 5.7.
// To run the script, it is necessary to load the main program from above.  

G:=SmallGroup(48,48);

SphGens2:=function(G,R)
Gens:={ };
for x1 in R do for x2 in R do for x3 in R do
for x4 in R do for x5 in R do for x6 in R do
for x7 in R do pr:=x1*x2*x3*x4*x5*x6*x7;
if pr in R then
if sub<G|x1,x2,x3,x4,x5,x6,x7> eq G then
Include(~Gens, [x1,x2,x3,x4,x5,x6,x7,pr^-1]);
end if; end if; end for; end for; end for; end for; end for; 
end for; end for; return Gens; end function;

B4:=SphGens2(G,ElementsOfOrd(G,2) diff Class(G,G.1));
B5:=SphGens2(G,ElementsOfOrd(G,2) diff Class(G,G.1*G.2));

OrbitRep:=function(W)
RMenge := { }; 
while not IsEmpty(W) do
Tup:=Rep(W); Include(~RMenge,Tup);   
orb2:=HurwitzOrbit(Tup);
for g2 in orb2 do Exclude(~W, g2); end for;
end while; return RMenge; end function;

R4:=OrbitRep(B4); R5:=OrbitRep(B5); 


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

// This is the source code to perform the computations in Remark 5.8.
// To run the script, it is necessary to load the main program from above.  

G:=SmallGroup(128,36);
test:=false;
A1:=[G.1*G.2*G.4*G.6,G.1*G.4*G.5*G.6,G.2*G.3*G.4*G.7];
B1:=A1;
A2:=[G.1*G.2*G.3*G.6*G.7,G.2*G.5*G.7,G.1*G.3*G.4*G.7]; 
B2:=[G.2*G.5*G.6*G.7,G.1*G.3*G.4*G.6*G.7,G.1*G.2*G.3*G.6*G.7];
Aut:=AutGr(G);
orb1:=HurwitzOrbit(A2); orb2:=HurwitzOrbit(A1);
for g1 in orb1 do for g2 in orb2 do for phi in Aut do
if [phi(g1),phi(g2)] eq [B1,B2] then 
test:=true;
break g1;
end if;
end for; end for;
end for;
test;