/* ===========================================================================================
In this file we provide a function which determines for a given finite group G and for a given 
type T the set of all generating vectors V for G of type T. 
Moreover, we provide an implementation of the action of the Braid group on the set of 
these generating vectors and, in the case that the types are [1; m] or [1; m, m], also the 
action of the full mapping class group. Recall that different generating vectors may determine 
threefolds with the same invariants. This happens if (but not only if) they differ by the 
actions from above (cf. Computational Remark 3.0.13). The code is used in my thesis and in the 
paper "Mixed Threefolds Isogenous to a Product". Any citations e.g. Remark 3.0.13 refer to my 
thesis. 
=========================================================================================== */

/* The script "ElsOfOrd" returns the set of elements of the group G of order d. */ 

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

/* Let G be a finite group and T=[h;m1,...,mr] be a type. The script "TuplesOfGivenOrder"
creates the set of all elements of the form (h1,...,hr,a1,b1,...,ah,bh), where hi is an 
element of G of order mi and ai, bi are arbitrary elements of G. */  

TuplesOfGivenOrders:=function(G,T)
h:=T[1]; SEQ1:=[ ];
	for i in [1..#T-1] do 
		El:=ElsOfOrd(G,T[i+1]);
			if IsEmpty(El) then return {}; 
				else Append(~SEQ1,El);
			end if; 
	end for;
SEQ2:=[ ];
	for i in [1..2*h] do 
		Append(~SEQ2,Set(G)); 
	end for;
return CartesianProduct(SEQ1 cat SEQ2); end function;

/* Let V=(h1,...,hr,a1,b1,...,ah,bh) be a tuple in the output of "TuplesOfGivenOrders(G,T)". 
The script "ProductOne" decides if the product h1*...*hr*[a1,b1]*...*[ah,bh] is trivial 
or not i.e. condition ii) in the definition of a generating vector (Definition 1.1.6). 
True or false is returned accordingly. */ 

ProductOne:=function(V,T,G);
test:=false; p:=Id(G); h:=T[1];
	for i in [1..#T-1] do p:=p*V[i]; end for;
	for i in [1..h] do
		a:=V[2*(i-1) + #T]; b:=V[2*i-1+#T]; p:=p*a*b*a^-1*b^-1;
	end for;
if p eq Id(G) then test:=true; end if;
return test; end function;

/* The next script allows us to test if a group G admits a generating vector V of type T. 
If a generating vector exists, then the function returns true and V, otherwise it returns 
false and the empty set. */  

ExistVectGens:=function(G,T)
n:=#G;  OneGV:={};
cands:= TuplesOfGivenOrders(G,T);
	for cand in cands do V:=[g : g in cand];
		if ProductOne(V,T,G) then 
			if #sub<G|V> eq  n  then 
				Include(~OneGV,V);
				break cand;
			end if; 
		end if; 
	end for; 
return not IsEmpty(OneGV), OneGV; end function;

/* The script "VectGens" determines the set of all generating vectors for the group G 
of type T. */ 

VectGens:=function(G,T)
n:=#G; Vect:={};
cands:= TuplesOfGivenOrders(G,T);
	for cand in cands do V:=[g : g in cand];
		if ProductOne(V,T,G) then 
			if #sub<G|V> eq  n  then 
				Include(~Vect,V);
			end if; 
		end if; 
	end for; 
return Vect; end function;

/* Given a generating vector V=(h1,...,hr,a1,b1,...,ah,bh) and an integer 1<= i <= r-1, the 
script "HurwitzMove"  performs the i-th Braid move i.e. it returns the generating vector 
(h1,...,h(i-1),hi*h(i+1)*hi^-1,hi,h(i+2)...,hr,a1,b1,...,ah,bh). */

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

/* Let V be a generating vector of type T for the group G. The script "BraidOrbit" determines the 
orbit of V under the action of the braid group Br. To spare memory only the elements in the orbit 
whose entries have increasing order are returned. */

BraidOrbit:=function(V,T)
orb:={ }; sortorb:={ }; Trash:={V}; 
repeat
    ExtractRep(~Trash,~gens); Include(~orb, gens);
		for i in [1..#T-2] 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..#T-2] 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;

/* Let V be a generating vector of type [1; m] or of type [1; m, m]. The script "MapGroupOrbit"
determines the orbit of V under the action of the mapping class group.
If the parameter "switch" is equal to 1, the case [1; m] is considered and otherwise the case 
[1; m, m]. Note that the action of the mapping class group is generated by 2 respectively 4 moves, 
wich are performed by the subscript "Moves". */ 

Moves:=function(V,switch)
if switch eq 1 then 
	t1:=[V[1],V[2]*V[1],V[3]]; t2:=[V[1]*V[2]^-1,V[2],V[3]];
	M:={t1,t2};
else 
	t1:=[V[1],V[2]*V[1],V[3],V[4]]; t2:=[V[1]*V[2]^-1,V[2],V[3],V[4]];
	t3:=[V[2]^-1*V[3]*V[1],V[2],V[3],((V[4]^V[1])^V[2])^(V[1]^-1)];
	t4:=[V[1]^-1,V[2]^-1,(V[4]^V[1])^V[2],(V[3]^V[2])^V[1]];
	M:={t1,t2,t3,t4};
end if;
return M; end function;

MapGroupOrbit:=function(V,switch)
orb:={};  Trash:={V};
repeat
    ExtractRep(~Trash,~gens); Include(~orb, gens);
		for newgens in Moves(gens,switch) do
			if newgens notin orb then Include(~Trash, newgens); end if; 
		end for;
until IsEmpty(Trash);
return orb; end function;

/* Let V be a generating vector for G of type T. The script "OrbitSwitch" decides which group 
action is applied: if the type of V is of the form [1,m] or [1,m,m], then the action 
of the full mapping class group is used, otherwise only the braid group is taken into account. */

OrbitSwitch:=function(V,T)
h:=T[1];
if #T eq 2 and h eq 1 then 
	orb:=MapGroupOrbit(V,1);                         
		elif #T eq 3 and h eq 1 and T[2] eq T[3] then 
			orb:=MapGroupOrbit(V,2);         
				else  
			orb:=BraidOrbit(V,T); 
end if;
return orb; end function;

/* Let G be a finite group and T=[h; m1,..., mr] be a type. 
i)  if T is of the form [h;-] (i.e. no ramification part), then the script "GenVecRep" returns 
a generating vector for G of type T, in case such a vector exists and otherwise the empty set.
ii) if T is not of the form [h;-], then the script "GenVecRep" returns a representative for each 
orbit of the action of the mapping class group, respectively the braid group, on the set of all 
generating vectors for G of type T. */

GenVecRep:=function(G,T)
h:=T[1];
if #T eq 1 then 
	test, Repres:=ExistVectGens(G,T); 
		else;
			gens:=VectGens(G,T); Repres:={};
				while not IsEmpty(gens) do 
					V:=Rep(gens); Include(~Repres,V); orb:=OrbitSwitch(V,T);
						for W in orb do 
							Exclude(~gens, W);
						end for; 
				end while; 
end if;
return Repres; end function;