// TuplesOfGivenOrder creates a sequence of the same length as the input
// sequence type, whose entries are subsets of the group 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-1] do
  EL:={g: g in G| Order(g) eq type[i]};
  if IsEmpty(EL) then return [{}];
  else Append(~SEQ,EL);
end if; end for;
return SEQ;
end function;

// This transforms a tuple into a sequence

TupleToSeq:=func<tuple|[x: x in tuple]>;

// Now we create all sets of spherical generators of a group of the
// prescribed signature.

VectGens:=function(G, type)
Vect:={}; SetCands:=TuplesOfGivenOrders(G,type);
for cand in CartesianProduct(SetCands) do
  if Order(&*cand) eq type[#type] then
    if #sub<G|TupleToSeq(cand)> eq #G then 
      Include(~Vect, Append(TupleToSeq(cand),(&*cand)^-1));
end if; end if; end for; ;
return Vect;
end function;

// HurwitzOrbit, starting from a sequence seq of elements of a group,
// creates all sequences of elements which are equivalent to the given one
// for the equivalence relation generated by the Hurwitz moves
// and returns (to spare memory) only the ones whose entries have 
// orders disposed as the ones in seq

HurwitzMove:= func<seq,j|Insert(Remove(seq,j),j+1,seq[j]^seq[j+1])>;

HurwitzOrbit:=function(seq)
orb:={ }; shortorb:={  }; Trash:={ seq };
repeat ExtractRep(~Trash,~gens); Include(~orb, gens);
  for k in [1..#seq-1] do 
    newgens:=HurwitzMove(gens,k);
    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..#seq] do 
    if not Order(gens[k]) eq Order(seq[k]) then 
      test:=false; break k;
  end if; end for;
  if test then Include(~shortorb, gens); end if;
end for;
return shortorb;
end function;

// OrbitsVectGens creates a sequence 
// whose entries are subsets of the group in the input,
// and precisely the subsets  are the orbits under the
// Hurwitz moves

OrbitsVectGens:=function(G, type)        
Orbits:=[];   Vects:=VectGens(G, type);
while not IsEmpty(Vects) do
  v:=Rep(Vects);  orb:=HurwitzOrbit(v);
  Append(~Orbits, orb); Vects:=Vects diff orb;
end while;
return Orbits;
end function;
 
// Disjoint checks if two spherical systems of generators are disjoint 

Stab:= function(seq, G)
M:= {} ;
for i in [1..#seq] do g:=seq[i]; 
  for n in [1 .. (Order(g)-1) ] do
    M := M join Conjugates(G,g^n) ;
end for; end for;
return M;
end function ;

Disjoint:=func<G,seq1, seq2 | IsEmpty( Stab(seq1,G) meet Stab(seq2,G))>;

// Homology computes the first homology group
// of the surface associated to the disjoint pair
// of spherical systems of generators (seq1, seq2) of G

Poly:=function(seq, gr) 
F:=FreeGroup(#seq);  Rel:={}; Q:=Id(F);
for i in {1..#seq} do 
  Q:=Q*F.(i); Include(~Rel,F.(i)^(Order(seq[i])));   
end for;
Include(~Rel,Q); P:=quo<F|Rel>; 	
return P, hom<P->gr|seq>;
end function;

Homology:=function(G, seq1,  seq2)
T1,f1:=Poly(seq1,G);   T2,f2:=Poly(seq2,G);
T1xT2:=DirectProduct(T1,T2);
 GxG,inG:=DirectProduct(G,G);
if Category(G) eq GrpPC then 
  n:=NPCgens(G); 
  else n:=NumberOfGenerators(G); 
end if;
Diag:= hom<G->GxG| [inG[1](G.j)*inG[2](G.j): j in [1..n]]>(G);
f:=hom<T1xT2->GxG| inG[1](seq1) cat inG[2](seq2)>;
Pi1:=Rewrite(T1xT2,Diag@@f); Pi1:=Simplify(Pi1);
return Pi1, AbelianQuotient(Pi1);
end function;

// Surfaces is the main function.
// It takes as input a group G and two signatures.
// It calls the previous function and moreover identifies
// distinct pairs of Hurwitz's orbits of spherical systems of
// generators under the action of Aut(G). 
// Note that we consider only the action of Out(G), since 
// the Inn(G)-action was already taken care of in HurwitzOrbit.

Surfaces:=function(G, type1,type2) 
R:={}; Aut:=AutomorphismGroup(G);
F,q:=FPGroup(Aut);     O1,p:=OuterFPGroup(Aut);
Out,k:=PermutationGroup(O1);
V1:=OrbitsVectGens(G, type1); 
if type1 eq type2 then V2:=V1;
  else V2:=OrbitsVectGens(G,type2); 
end if;
W:=Set(CartesianProduct({1..#V1},{1..#V2}));
for pair  in W do 
  v1:= Rep(V1[pair[1]]); v2:= Rep(V2[pair[2]]);
  if not Disjoint(G, v1,v2) then Exclude(~W, pair);
 end if;
end for;
while not IsEmpty(W) do 
  pair:=Rep(W);    
  v1:= Rep(V1[pair[1]]); v2:= Rep(V2[pair[2]]);Include(~R,[v1,v2]);
  for phi in Out do 
    if exists(x){y: y in W | (q(phi@@(p*k)))(v1) in V1[y[1]] and 
   	 (q(phi@@(p*k)))(v2) in V2[y[2]]} then Exclude(~W, x); 
      if type1 eq type2 then Exclude(~W, <x[2],x[1]>); 
    end if; end if;
    if IsEmpty(W) then break phi; end if;
end for; end while;
printf "Number of components: %o\n", #R;
for r in R do  Pi1, H1:= Homology(G,r[1],r[2]);
	printf "Spherical generators:\n%o\n%o\n", r[1],r[2];
	printf "Homology:\n%o\n", H1;
	//printf "FundamentalGroup:\n%o\n", Pi1;
end for;
return {};
end function;