/* The main functions of the script are
1) the procedure OrbitDivisors, and
2) the functions "q1Extension" and "q2Extension",  which are defined at the bottom of the file.
*/

//-------------------------//-------------------------//-------------------------
//-------------------------//-------------------------//-------------------------

// PART 1: Orbit Divisors

// This function counts the number of fixpoints of g1^-1*g2 in H < Aut(C)

CountingIntersections:=function(g1,g2,seq, H)
int:=0;
for j in [1..#seq] do 
	c:=0; K:=sub<H|seq[j]>;
	for g in H do 
		if g1^-1*g2 in  {g*k*g^-1: k in K}  then c:=c+1; 
		end if;
	end for; 
	int:=int+(c/#K); 
end for;
return int;
end function;

/*
The following procedure is the main function of the script.
Given: 
i) the group G,
ii) the generating vector gv for G0 determing the covering C->C/G^0=C' with g(C')=q,
ii) the automorphism group H: G0<H<Aut(C), the embedding e:G0 -> H, and
iii) the ramification part w of a generating vector for the covering C->C/H, 
it constructs all the orbit divisors on CxC/G induced by elements of H and prints their intersection products.
*/


OrbitDivisors:=procedure(G,genus,G0,H,e,w)

genus_1:=genus- 1; // g(C)-1



t:=Rep({x: x in G | x notin G0}); // \tau'
O2:={x: x in G| x notin G0 and Order(x) eq 2};

R:={e(t*x): x in O2}; //the ramification curves represented by elements in e(G0) < H

RamCurves:=R;

BranchCurves:=[];
OtherCurves:=[];

P:=Set(H);

while not IsEmpty(R) do
	f:=Rep(R);  //it represents the curve Delta_f, graph of the automorphism f of C
	Gamma:={e(t*h*t^-1)*f*e(h^-1) : h in G0} join {e(t^2*h)*f^-1*e(t*h^-1*t^-1) : h in G0}; // it is the G-orbit  of Delta_f
	Append(~BranchCurves, Gamma);
	R:=R diff Gamma; P:=P diff Gamma;
end while;  //Each element in BranchCurves represents a branch curve

#BranchCurves, "Branch curves";

while not IsEmpty(P) do
	f:=Rep(P); 
	Gamma:={e(t*h*t^-1)*f*e(h^-1) : h in G0} join {e(t^2*h)*f^-1*e(t*h^-1*t^-1) : h in G0}; //the  Curves \gamma(Delta_f), \gamma in G
	Append(~OtherCurves, Gamma);
	P:=P diff Gamma;
end while; //Each element in OtherCurves represents an orbit divisor, which is not a branch curve

#OtherCurves, "Other curves"; 

Curves:= BranchCurves cat OtherCurves; 
//each element in Curves is an orbit of the G-action on the set {Delta_f |f in H}, and represent an orbit divisor.

#Curves, "Orbit divisors"; 


Curves_1:={};
SelfInt:=[]; IntK:=[];

for a in [#BranchCurves+1..#Curves] do 
	Gamma:=Curves[a]; n:=#Gamma; 	
	self:= -2*n* genus_1; GK:=4*n*genus_1;
	for g1 in Gamma do
		for g2 in Gamma diff {g1} do 
			int:=CountingIntersections(g1,g2,w,H);
			self:=self+int;
		end for; 
		for r in RamCurves   do 
			intK:=CountingIntersections(g1,r,w, H);
			GK:=GK-intK; 
	end for; end for;
	Csquare:=self/#G; CK:=GK/#G;
	Append(~SelfInt, Csquare); Append(~IntK, CK);
	if Csquare eq -1 and CK eq -1 then
		Include(~Curves_1, a);
	end if;
end for;

 #Curves_1, "(-1)-curves"; // it returns the number of (-1)-curves.

end procedure;

//-------------------------//-------------------------//-------------------------
//-------------------------//-------------------------//-------------------------

// PART 2: Extensions



/* As explained in [CF18, Section 6.2], (see also [BC04, BCG08]), generating vectors for a group H, 
of type (0; m_1,..., m_r), which differ by a Huwitz move or an automorphism of H yield isomorphic covers
C->C/H = P^1, hence we shall consider only one representative for each orbit */
HurwitzMove0:= func<seq,idx|Insert(Remove(seq,idx),idx+1,seq[idx]^seq[idx+1])>;

// HurwitzOrbit0, starting from a sequence 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 of type (0; m_1,..., m_r),
// and return (to spare memory) only the ones whose entries have never decreasing order.

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

// We construct the different orbits.

Orbit:=function(H,GV)      
RepOrb:={}; 
AutH:= AutomorphismGroup(H);
f,A:=PermutationRepresentation(AutH);
while not IsEmpty(GV) do 
	w:=Rep(GV); Include(~RepOrb,w);
	HOrb:=HurwitzOrbit0(w);
	for v in HOrb do  for a in A do
		Exclude(~GV,(a@@f)(v));
		if IsEmpty(GV) then break v; end if;
	end for;end for; 
end while;
return RepOrb;
end function;



// The Extension H of G^0 we are looking for is a semidirect product: H= G^0 \rtimes Z/2.
// Two such semidirect products given by f,g:Z/2 -> Aut(G^0)  are isomorph, 
// if f(1) and g(1) are in the same conjugacy class of Aut(G^0). 
// The next function takes care of this


AutConjugCl:=function(Aut, order)
Set:={};  Rep:=[];
D:=Divisors(order);
f,A:=PermutationRepresentation(Aut);
list:=[x: x in A | Order(x) in D]; 
for el in list do
	if el notin Set then         
		l:={};
		for a in A do Include(~Set, el^a); Include(~l, (el^a)@@f);
		end for;
		Append(~Rep, l);
 end if; end for;
return Rep;
end function;

// This transforms a  tuple  into a sequence

TupleToSeq:=function(tuple)
seq:=[ ];
for el in tuple  do 	Append(~seq,el); end for;
return seq;
end function;



/* The next functions provides the group H, the embedding e:G0-> H, 
and a generating vector w for the covering C->C/H,  starting form a generating vector for the covering C->C/G^0.

The reason to considering Hurwitz orbits is explained in the files with the input data.
*/

q1Extension:=function(G0,gv)

AutG0:= AutomorphismGroup(G0);
R2:= AutConjugCl(AutG0,2); 
C2:=SmallGroup(2,1);
checked:={};

for i in [1..#R2] do  
	autom2:=Rep(R2[i]);
	f:=hom<C2->AutG0|autom2>;
	H,e:= SemidirectProduct(G0,C2,f); // we construct a semidirect procuct H= G^0 \rtimes Z/2.
	if IdentifyGroup(H)[2] in checked then;
	else
		q1:=e(gv[3]); q2:=e(gv[4]);  Include(~checked, IdentifyGroup(H)[2]); 
		if  IsConjugate(H,q1,q2) then // the ramification elements in the original generating vector must be conjugated in H 
			P:={x: x in H | x notin e(G0) and Order(x) eq 2};
			cands:=CartesianProduct([{q1,q2},P,P,P,P]);  
			GV:={};
			for c in cands do
				if Order(&*c) eq 1 then
				if #sub<H|TupleToSeq(c)> eq #H then 
					Include(~GV,TupleToSeq(c));
			end if; end if; end for; // we check if H has a g.v. of the expected type
			if #GV ne 0 then 
				HO:=Orbit(H,GV); break i; // we have found the extension H
			end if;
end if; end if; end for;

return H,e,HO;
end function;
	



q2Extension:=function(G0)

AutG0:= AutomorphismGroup(G0);
R2:= AutConjugCl(AutG0,2);
C2:=SmallGroup(2,1);

checked:={};
for i in [1..#R2] do 
	autom2:=Rep(R2[i]);
	f:=hom<C2->AutG0|autom2>;
	H,e:= SemidirectProduct(G0,C2,f); // we construct a semidirect procuct H= G^0 \rtimes Z/2.
	if IdentifyGroup(H)[2] in checked then;
	else
		Include(~checked, IdentifyGroup(H)[2]); 
		P:={x: x in H | x notin e(G0) and Order(x) eq 2};
		cands:=CartesianProduct([P,P,P,P,P,P]);
		GV:={};
		for c in cands do
			if Order(&*c) eq 1 then
			if #sub<H|TupleToSeq(c)> eq #H then 
				Include(~GV,TupleToSeq(c));
		end if; end if; end for;
		if #GV ne 0 then 
			HO:=Orbit(H,GV); break i; // we have found the extension H
		end if;
end if; end for;

return H,e,HO; 
end function;