/* *******************************************************************************
   To compute unmixed product-quotient surfaces with abelian group, at most canonical 
   singularities, geometric genus "pg", and with canonical map having at most "nrb" 
   base points, run "ListSurf(pg,nrb)". 
   The surfaces are stored in the file "CanBase_pg_pg_nrb_nrb*.txt". 
********************************************************************************** */ 
 
 /* *******************************************************************************
   ******************************* PART I ****************************************
   ***************************** Numerology **************************************
   ******************************************************************************* */

/* Let "T:=[m1,..., mr]" be a type. The script "Theta" determines the rational number 
   Theta(T) appearing in Hurwitz' formula */

Theta:= function(T)
	t:= -2;
	for i in [1..#T] do 
		m:=T[i]; t:= t + 1 - 1/m;
	end for;
return t; end function; 

/* The script "OrderedSeq" transforms a (multi-) set or a sequence into an ordered 
   sequence. */

OrderedSeq:=function(mset)
seq:=[ ];
	while #mset ne 0 do 
		Append(~seq, Minimum(mset)); Exclude(~mset, Minimum(mset)); 
	end while;
return seq; end function;

/* if T=[m1,...,mr] is the type of a generating vector of an abelian group. Then  
   it holds: 

         lcm(m1,...m_{i-1},m_{i+1},...,mr)= lcm(m1,...,mr) for all 1<=i<=r. 
   
   The function "TestAbel" returns "true" or "false", accordingly. */

TestAbel:=function(T)
r:=#T; test:=true; m:=Lcm(T);
	for i in [1..r] do
		S:=[T[l]: l in [1..i-1] cat [i+1..r]];
			if Lcm(S) ne m then
				test:=false; break i;
			end if;
	end for;
return test; end function;

/* Let "T=[m1,...,mr]" be the type of a generating vector of an abelian group "G"
   of length "r" at least four. Then the order of the group "GOrd" is bounded by "4(g-1)", 
   except for "T=[2,2,3,3]" (here "G" is must be cyclic of order six, 
   see Bauer, Catanese "Some new surfaces with pg=q=0" Proposition 2.1).
   The function "TestAbel" returns "true" or "false", accordingly. */ 

TestFour:=function(T,GOrd,g)
test:=true;
if T ne  [2,2,3,3] then 
	if #T ge 4 and GOrd gt 4*(g-1) then 
		test:=false;
	end if;
	else
if GOrd ne 6 then 
	test:=false;
end if;
end if;
return test;
end function;

/* The script "ListOfTypes" determines the set of types T=[m1,..., mr] such that:
   i)   the mi's are bounded by 4g + 2 (Wiman),
   ii)  the mi's are divisors of "GOrd" and "d", 
   iii) the length r of T is bounded from above by "4*(g-1)/GOrd + 4"
   iv)  Hurwitz' formula 2g-2 = GOrd*Theta(T) holds. 
   iiv) the functions "TestAbel(T)" and "TestFour(T)" return "true" */

ListOfTypes:=function(GOrd,g,d) 
Types:= {}; Orders:={k: k in Divisors(Gcd(GOrd,d)) | 2 le k and k le 4*g +2}; 
MaxRam:= Floor(4*(g-1)/GOrd + 4);        
	for r in [3..MaxRam] do M:=Multisets({x : x in Orders},r);
		for mset in M do T:=OrderedSeq(mset);
			if GOrd*Theta(T) eq (2*g-2) then  
			    if TestAbel(T) and TestFour(T,GOrd,g) then 
					Include(~Types,T);
				end if; 
			end if;
	end for; end for;
return Types; end function;      

/* The script "Threetup" determines all possible triples of the form 
   [GOrd,g1,g2] such that:   

i)   chi*GOrd >=(g1-1)*(g2-1), where chi:=1+pg
ii)  8(g1-1)(g2-1)/GOrd is an integer 
ii)  g2 <= 1+2*chi + Sqrt(4*chi^2 + 8*chi)
iii) g1 <= 4*chi*(g2+1)/(g2-1)+1
iv)  GOrd <= Min(4*(g1+1),4*(g2+1))

*/

Threetup:=function(pg)
chi:=pg+1;
Set:={};
g2max:=Floor(1+2*chi + Sqrt(4*chi^2 + 8*chi));
for g2 in [2..g2max] do
	g1max:=Floor(4*chi*(g2+1)/(g2-1)+1);
		for g1 in [g2..g1max] do
			for GOrd in [2..Minimum(4*(g1+1),4*(g2+1))] do 
			if 8*(g1-1)*(g2-1)/GOrd in IntegerRing() and 
				GOrd*chi ge (g1-1)*(g2-1) then 
					Include(~Set,[GOrd,g1,g2]);
			end if; 
			end for;
		end for;
end for;
return OrderedSeq(Set); end function;

/* For each triple "[GOrd,g1,g2]" in the List "Threetup(pg)" we determine 
the possible types "T1=(m1,...,mr)" and "T2=(n1,...,ns)" such that
lcm(m1,...,mr)=lcm(n1,...,ns). This number is equal to the exponent of the group. 
For each of these pairs, we compute all abelian groups admitting that 
exponent. The tuples [T1,T2,"Group Id"] are returned. 
Moreover, the fact that "K_X" is a line bundle implies that we just need the types where 
"mi" divides "2(g2-1)" and "nj" divides "2(g1-1)". */
 
TypesGroupCand:=function(pg)
Cand:=[];
Three:=Threetup(pg);
for S in Three do
	GOrd:=S[1]; g1:=S[2]; g2:=S[3];
	Types1:=ListOfTypes(GOrd,g1,2*(g2-1)); 
	Types2:=ListOfTypes(GOrd,g2,2*(g1-1));
		for T1 in Types1 do
			ex:=Lcm(T1);
			for T2 in Types2 do
				if ex eq Lcm(T2) then  
					for G in SmallGroups(GOrd: Warning:=false) do
						if IsAbelian(G) and Exponent(G) eq ex then
							IdG:=IdentifyGroup(G);
							Tup:=[T1,T2,[IdG[1],IdG[2]]];
							Append(~Cand,Tup);
						end if;
					end for;
				end if;
			end for;
		end for;
end for;
return Cand; end function;

/* *******************************************************************************
   ******************************* PART II ***************************************
   ******************* Generating Vectors, group actions ************************** 
   ******************** and the classification function **************************
   ******************************************************************************* */

/* Given a multiset of elements "V" in an abelian group "G", we determine their 
   product and return "true" in case it is the identity and "false" otherwise. */

ProdO:=function(V,G)
test:=false; h:=Id(G);
	for g in V do
		h:=g*h;
	end for;
if Order(h) eq 1 then 
	test:=true;
end if;
return test; end function;

/* The next function turns a type "T=[m1,...,mr]" into the multiset "{* m1,...,mr *}" */

MSet:=function(T)
MultS:={* *};
	for g in T do
		Include(~MultS,g);
	end for;
return MultS; end function;

/* 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}>;

/* The function "VectGensAbel" determines all generating vectors of type "T" for the abelian 
   group "G". We view and save them as multisets and not as vectors, since we assume that 
   the group is abelian which means that the order is not important. */
	
VectGensAbel:=function(G,T)
MsetVect:={}; n:=#G;

Mu:=MSet(T); Se:=Set(Mu); List:=[];
for i in Se do
	Append(~List,Multisets(ElsOfOrd(G,i),Multiplicity(Mu,i)));
end for;
Cart:=CartesianProduct(List);

for C in Cart do
	V:={* *};
		for c in C do 
			V:=V join c;
		end for;
			if ProdO(V,G) then 
				if #sub<G|Set(V)> eq  n  then 
					Include(~MsetVect,V);
				end if; 
			end if; 
end for; 
return MsetVect; end function;				

/* To find the components of the moduli space we have to determine the orbits of the action 
   of "Aut(G)" on the set of pairs "{*V1,V2*}" of generating "multisets" or generating "Vectors". 
   The function "ConjPair" returns the conjugate "{*phi(V1),phi(V1)*}" for a given pair "{*V1,V2*}" 
   and a given automorphism "phi". */

ConjPair:=function(phi,Pair);
NewPair:={* *};
for V in Pair do
	Vnew:={* *};
		for g in V do 
			Include(~Vnew,phi(g)); 
		end for;
	Include(~NewPair,Vnew);
end for;
return NewPair;
end function;

/* Given two sets of spherical generators, the singular points of the resulting surface are the 
   image of points in the product of curves "A x B" having nontrivial stabilizer. These correspond 
   to tuples (g_1,n_1,g_2,n_2), where
   
   i)   g_1 is a generator of the first set
   ii)  g_2 is a generator of the second set
   iii) 1<=n_1<=ord(g_1)
   iv)  1<=n_2<=ord(g_2) 
   v)   g_1^n_1=g_2^n_2

   The first function "BasketByAPairOfGens" computes the singular points coming from a fixed pair 
   (g1,g2). The second function "BasketOfSings" determines the basket of singularities of every 
   constructed surface. */

BasketByAPairOfGens:=function(group,gen1,gen2)
basket:={* *}; RC:={ }; delta:=GCD(Order(gen1),Order(gen2));
alpha1:=IntegerRing()!(Order(gen1)/delta);
alpha2:=IntegerRing()!(Order(gen2)/delta);
RC2,f2:=RightTransversal(group,sub<group | gen2 >);
for g2 in RC2 do test:=true;
for g in sub<group| gen1 > do
if f2(g2*g) in RC then test:=false; break g;
end if; end for;
if test then Include(~RC, g2);
end if; end for;
for g in RC do
for d1 in [1..delta-1] do
for d2 in [1..delta-1] do
if (gen1^(d1*alpha1)) eq (gen2^(d2*alpha2))^g then
Include(~basket,d2/delta); break d1;
end if; end for; end for; end for;
return basket;
end function;

BasketOfSings:=function(Pair,G)
basket:={* *};
V1:=Pair[1];
V2:=Pair[2];
for gen1 in V1 do
for gen2 in V2 do 
	basket:=basket join BasketByAPairOfGens(G,gen1,gen2); 
end for; end for;
return basket;
end function;

/* The function "Recomputepg" recomputes the geometric genus of a constructed surfaces, from the 
   genera of the curves and the basket of singularities. */ 

Recomputepg:=function(GOrd,g1,g2,Basket)
pg:=(g1-1)*(g2-1)/GOrd -1;
for s in Basket do
n:=Denominator(s);
pg:=pg + (n^2-1)/(12*n);
end for;
return pg;
end function;

/* Still, some of the surfaces may have non-canonical singularities. We use the function 
   "Test_Can" to verify if the surface has canonical singularities, it returns and "true" or 
   "false" accordingly. */

Test_Can:=function(Basket) 
testCan:=true; testTer:=true; 
for s in Basket do
a:=Numerator(s);
n:=Denominator(s);
	if (1+a)/n notin IntegerRing() then 
		testCan:=false;
		break s;
	end if;
end for;
return testCan;
end function;

/* ******************************* PART III **************************************
**************** ************The character sheaves *******************************
********************************************************************************** */

/* The next function is used to determine the tensor product of "Lchi" with the 
   canonical sheaf of "Y:=P^1 \times P^1" */ 


LchiTensorKY:=function(V1,V2,chi)
m:=0; n:=0; d:=Order(chi);
F:=CyclotomicField(d); z:=F.1; 
	for g1 in V1 do
			for u1 in [0..d-1] do
				if chi(g1) eq z^u1 then 
					m:=m+u1;
				end if;	
			end for;
	end for;
	
	for g2 in V2 do
			for u2 in [0..d-1] do
				if chi(g2) eq z^u2 then 
					n:=n+u2;
				end if;	
			end for;
	end for;	

	return m/d-2, n/d-2;
end function;

/* The next step is to analyse the canonical system. First, the sheaves 
   "M_chi:=L_chi \otimes K_Y" with non-zero cohomology are determined. By the general theory, 
   the canonical system is generated by the subsystems "|s_chi \cdot \psi^*H^0(M_chi)|", where 
   the support of the divisor "\Phi_chi" of the section "s_chi" is contained in the ramification 
   locus of "\psi". The function "Sections_xy" determines the "vertical" and "horizontal" 
   components of this divisor for all "chi". */ 

LchiNonZeroCoh:=function(G,V1,V2)
CT:=CharacterTable(G);
Liste:=[];
for i in [1..#CT] do
    chi:=CT[i];
        a,b:=LchiTensorKY(V1,V2,chi);
        if a ge 0 and b ge 0 then 
            Append(~Liste,[a,b,i]);
        end if;
end for;
return Liste;
end function;

PartOfBasis:=function(V,chi) 
List:=[]; d:=Order(chi); 
F:=CyclotomicField(d); z:=F.1; 
	for g in V do
		m:=Order(g);
			for u in [0..d-1] do
				if chi(g) eq z^u then
					Append(~List,m-1-u*m/d);
				end if;
			end for;
	end for;
return List; end function;

Sections_xy:=function(G,V1,V2)
S:=[]; CT:=CharacterTable(G);
List:=LchiNonZeroCoh(G,V1,V2);
for D in List do
    dime:=IntegerRing()!(D[1]+1)*(D[2]+1);
        i:=IntegerRing()!D[3];
        chi:=CT[i];
        Append(~S,[PartOfBasis(V1,chi),PartOfBasis(V2,chi)]);
end for;
return S;
end function;

/* Since the linear systems "|M_chi|" on "P^1 \times P^1" are base-point-free, the 
   canonical system of the product-quotient surface "S" has at most isolated base-points 
   if and only if the divisors  "\Phi_chi" for which "|M_chi|" is non-empty meet only in 
   finitely many points (i.e. the canonical system of the quotient model "Y" is b.p.f.)
   and non of the intersection points is contained in the singular locus of the 
   quotient model "Y" (i.e. no -2 curve is a base curve of |K_X|). 
   The function "TestBaseCurves" is used to decide if "|K_Y|" has base curves. For the 
   surfaces "Y" with no base curves, the function "IsoBasePoints" determines the base-points 
   of "|K_Y|", more precisely it produces a list of triples "[i,j,nr]", where the intersection 
   of the vertical and horizontal components "E_i" and "E_j" yield "nr" base-points. Apart 
   from this list, it returns the total number of base-points and a boolean variable, that tells 
   us if "|K_X|" has isolated base-points. */
   
TestBaseCurves:=function(G,V1,V2)
test:=false;
Bxy:=Sections_xy(G,V1,V2);

for i in [1..#V1] do
		for A in Bxy do
		cou:=0;
			if A[1][i] ge 1 then 
				cou:=cou+1;
			end if;
		if cou eq #Bxy then 
			test:=true;
			break i; 
		end if;
		end for;
end for;
if test eq false then 
	for j in [1..#V2] do
		for B in Bxy do
		cou:=0;
			if B[2][j] ge 1 then 
				cou:=cou+1;
			end if;
		if cou eq #Bxy then 
			test:=true;
			break j; 
		end if;
		end for;
	end for;
end if;
		
return test;
end function;


IsoBasePoints:=function(G,V1,V2)
iso:=true;
Bxy:=Sections_xy(G,V1,V2);
S:=[];
for i in [1..#V1] do
    for j in [1..#V2] do
        test:=true;
            for B in Bxy do
                if B[1][i] eq 0 and  B[2][j] eq  0 then 
                    test:=false;
                    break B;
                end if;
            end for;
        if test eq true then 
			H:=sub<G|V1[i]>; K:=sub<G|V2[j]>; 
			nr:=#G/(#H*#K);
			Stab:=H meet K; 
			if #Stab ne 1 then  // i.e. the intersection point is non-singular 
				iso:=false;
			end if;
			 Append(~S,[i,j,nr]);
		end if;
end for; end for;

b:=0;
for base in S do 
	b:=b+base[3];
end for;


return S, iso, b;
end function;

/* We use the function "OrbitsCanSing" search for pairs of spherical generators which
correspond to canonical product quotient surfaces with group "G" and types "T_1" and "T_2". 
We need to take care of repetitive output i.e. pairs corresponding to isomorphic surfaces. 
By the general theory, two surfaces are isomorphic if the pairs of spherical generators are 
equivalent for the equivalence relation generated by Hurwitz moves (on each set of generators 
separately) and the automorphisms of "G" (simultaneously on both). Since we work in the abelian 
category, the Hurwitz moves just act as permutations, hence regarding a generating vector as a 
set instead of an ordered tuple, we don't need to implement the Hurwitz moves and just consider 
simultaneous conjugation by the automorphisms of "G". Unfortunately, to iterate directly over 
"Aut(G)" is not possible. A solution to overcome this problem is to collect all automorphism 
as a set. However, this can take very long, as "Aut(G)" is in general very large. Our solution 
is to construct a permutation representation of "Aut(G)" over which we can iterate. In comparison 
this method is much faster then the previous. the output of the function consists of all non-
equivalent pairs for which the corresponding surface has only canonical singularities. */

OrbitsCanSing:=function(G,T1,T2)
Orb:={};
	Vec2:=VectGensAbel(G,T2);
		if not IsEmpty(Vec2) then 
			Vec1:=VectGensAbel(G,T1);
			C:=CartesianProduct([Vec1,Vec2]);
			Pairs:={{*V[1],V[2]*} : V in C};
			Aut:=AutomorphismGroup(G); BadSet:={};
			f,P:= PermutationRepresentation(Aut); // Since we cannot iterate over Aut(G) directely,
				while not IsEmpty(Pairs) do	      // we construct a permutation representation. 
					Tup:=Rep(Pairs);		      // This is much faster than the 
					Include(~Orb,Tup);		      // previous method of memorizing Aut(G) as a set!
						for p in P do
							phi:=p @@ f; 
							Exclude(~Pairs, ConjPair(phi,Tup));
								if IsEmpty(Pairs) then 
									break p;
								end if;
						end for;
				end while;
		end if;
Orb:={[V : V in Pair] : Pair in Orb};
Orb:={Pair : Pair in Orb | Test_Can(BasketOfSings(Pair,G))};
return Orb; 
end function; 							

/* ********************************* PART IV ************************************
**************************** The Classification function ************************
******************************************************************************* */

/* This is the main function, it classifies all Product Quotients with canonical 
   singularities, abelian Galois group and a fixed value of "pg" such that the 
   canonical map has at most "nrb" base points. 
*/

ListSurf:=function(pg,nrb)	
DGRS:={};
F:="CanBase_pg_" cat IntegerToString(pg) cat "_nrb_" cat IntegerToString(nrb) cat ".txt";
fprintf F, "Canonical Product Quotients with abelian Galois group and pg=%o
with canonical system having at most %o isolated base-points. 
\n\n\n\n\n", pg, nrb;
count:=1;
LNDG:=TypesGroupCand(pg);
ListNum:=[]; LOrb:=[* *];
for D in LNDG do
	T1:=D[1]; T2:=D[2]; GOrd:=D[3][1];
	G:=SmallGroup(D[3][1],D[3][2]);
	g1:=Theta(T1)*GOrd/2 +1;
	g2:=Theta(T2)*GOrd/2 +1;
	Orbits:=OrbitsCanSing(G,T1,T2);
		for Pair in Orbits do 
			Basket:=BasketOfSings(Pair,G);
			Vol:=8*(g1-1)*(g2-1)/GOrd;
				if Recomputepg(GOrd,g1,g2,Basket) eq pg then 
					V1:=[g : g in Pair[1]];
					V2:=[g^-1 : g in Pair[2]];
					NonSheaves:=LchiNonZeroCoh(G,V1,V2);
						if TestBaseCurves(G,V1,V2) eq false then 
							BPoints, iso, b:=IsoBasePoints(G,V1,V2);
								degr:=(Vol-b)/(pg-2);
									if iso and b le nrb then
										fprintf F, "%o) Group: G= %o\n",count, IdentifyGroup(G);
										fprintf F, "\n\n Abelian Invariants: %o\n", AbelianInvariants(G);
										fprintf F, "\n\n T1=%o,  T2=%o\n", T1, T2;
										fprintf F, "\n\n Generating-Vectors and the character sheaves F_chi \n";
										fprintf F, "\n V1=%o, \n V2=%o  \n", Pair[1], Pair[2];
										fprintf F, "\n The sheaves F_chi with non-zero H^0: %o   \n", NonSheaves;
										fprintf F, "\n Base-points: %o   \n", BPoints;
										fprintf F, "\n Singularities: %o   \n", Basket;
										fprintf F, "\n Upper bound for the degree of the canonical map: %o   \n", degr;
										fprintf F, "\n\n\n\n\n";
										count:=count+1;
									end if;
						end if;
				end if;
		end for; 
end for;
return "Classification finished for pg=" cat IntegerToString(pg);
end function;