/* ====================================================================================================
In this file we provide a function to determine the basket of singularities of a product quotient 
threefold X from an algebraic datum of X. We also provide a function to decide whether a cyclic 
quotient singularity is canonical or terminal and a function to determine the basket of a crepant 
terminalisation Xhat of a product quotient threefold X with canonical singularities (cf. Remark 4.3.2).
===================================================================================================== */

/* The script "Sings_gen_of_stab" is used to determine the singularities of X which are contained in 
the intersection of the fibres of X --> Ci/G over the branch points qi of the covers Ci --> Ci/G (see 
Proposition 4.4.3). The input of the function consists of the group G and the hi's in the generating 
vectors Vi which correspond to the branch points qi. The singularities are returned as a "multiset", where 
an element [m,a,b] denotes a singularity of type 1/m(1,a,b). The script "Basket_of_Sing" allows us to 
determine the basket of singularities of X: we simply run over all h1, h2 and h3 in (the ramification parts) 
of the generating vectors V1, V2 and V3, apply the previous function and collect all the singularities. */
 
Sings_gen_of_stab:=function(G,h1,h2,h3)
basket:={* *}; RC:={ }; 
RC2,f2:=RightTransversal(G,sub<G|h2 >); RC3,f3:=RightTransversal(G,sub<G|h3 >);
Cart:=CartesianProduct(RC2,RC3);
for ca in Cart do test:=true; 
	for g in sub<G | h1 > do
		if <f2(ca[1]*g),f3(ca[2]*g)> in RC then 
			test:=false; 
			break g;
		end if; 
	end for;
	if test then 
		Include(~RC, ca);
    end if; 
end for;
	for b in RC do
		for n1 in [1..Order(h1)-1] do
		for n2 in [1..Order(h2)-1] do
		for n3 in [1..Order(h3)-1] do
			x:=(h1^n1);
				if x eq (h2^n2)^b[1] and x eq (h3^n3)^b[2] then 
					m:=Order(x);
						Include(~basket,[m, IntegerRing()!(m*n2/Order(h2)), IntegerRing()!(m*n3/Order(h3)) ]); 
						break n1;
				end if; 
		end for; end for; end for;end for;
return basket; end function;

Basket_of_Sing:=function(G,V1,V2,V3,T1,T2,T3)
Basket:={* *};
	for n1 in [1..(#T1-1)] do
	for n2 in [1..(#T2-1)] do
	for n3 in [1..(#T3-1)] do
		Basket:=Basket join Sings_gen_of_stab(G,V1[n1],V2[n2],V3[n3]);
	end for; end for; end for;
return Basket; end function;

/* The function "Test_Can_Ter" allows us to decide if a singularity "S" is canonical or terminal(see Proposition 
4.2.11). The output consists of two "booleans" i.e. true or false. The first value tells us whether the singularity 
is canonical, the other value whether the singularity is terminal. Recall that a terminal singularity is also 
canonical. Therefore, the output "false, true" is not possible. */

Test_Can_Ter:=function(S) 
m:=IntegerRing()!S[1]; a:=IntegerRing()!S[2]; b:=IntegerRing()!S[3]; 
testCan:=true; testTer:=true; 
	for k in [1..m-1] do
		x:=k mod m; y:=k*a mod m; z:=k*b mod m;
			if x + y + z lt m then 
				testCan:=false; testTer:=false; break k;
			elif x + y + z eq m then 
				testTer:=false; 
			end if;
	end for;
return testCan, testTer; end function;

/* The script "Test_Basket" verifies if the basket of singularities of a product quotient threefold contains 
only canonical singularities. Here, we simply run over the basket and apply the function "Test_Can_Ter" to each 
singularity. */

Test_Basket:=function(Basket)
test:=true;
	for S in Basket do
		testCan, testTer:=Test_Can_Ter(S);
			if testCan then 
				test:= true;
			else 
				test:=false; break S;
			end if;
	end for;
return test; end function;

/* Let a and m be coprime elements. The function "InvMod" determines the inverse of a modulo m. */ 

InvMod:=function(a,m)
	c, b:=XGCD(a,m); b:=b mod m;  
return b; end function; 

/* Given a terminal cyclic quotient singularity 1/m(1,x,y), the function "ter" returns, among all isomorphic 
singularities of the form 1/m(1,a,m-a) (see Example 4.2.2) the one, where a is minimal. The output is given 
as a fraction a/m, which is possible since a and m are coprime. 
The function "gor" does the same in the case of a Gorenstein singularity. Recall that Gorenstein cyclic quotient
singularities are of the form 1/m(1,a,m-a-1). */

ter:=function(S)
m:=S[1]; x:=S[2]; y:=S[3];
Set:={[x,y],[InvMod(x,m) mod m, y*InvMod(x,m) mod m], [InvMod(y,m) mod m, x*InvMod(y,m) mod m]};
a:=m+1;
	for Z in Set do
		x:=Z[1]; y:=Z[2];
			if y eq m-x then          // for a terminal singularity 1/m(1,x,y) it does not automatically hold that y=m-x
				a:=Minimum({a,x,y});
			end if;
	end for;
return a/m; end function; 

gor:=function(S)
m:=S[1]; x:=S[2]; y:=S[3];
Set:={[x,y],[InvMod(x,m) mod m, y*InvMod(x,m) mod m ], [InvMod(y,m) mod m, x*InvMod(y,m) mod m]};
a:=m+1;
	for Z in Set do
		x:=Z[1]; y:=Z[2]; a:=Minimum({a,x,y});
	end for;
return a/m; end function;

/* According to Theorem 4.2.13 an isolated canonical cyclic quotient singularity is either terminal, Gorenstein 
or it is of type III. The script "Split_up_Basket" divides the Basket of singularities of our product quotient 
threefold in three disjoint subsets according to these possibilities. */

Split_up_Basket:=function(Basket)
BaskTer:={* *};  BaskGor:={* *}; BaskTypeIII:={* *}; 
	for S in Basket do	
		testCan, testTer:=Test_Can_Ter(S);
			if testTer then 
				Include(~BaskTer,ter(S));
			elif (1+S[2]+S[3])/S[1] in IntegerRing() then      
				Include(~BaskGor,gor(S));
			else 
				Include(~BaskTypeIII,S);
			end if;
	end for;
return BaskTer, BaskGor, BaskTypeIII; end function;

/* The script "Xhat_Basket" determines the basket of a crepant terminalisation of a product quotient threefold X 
with canonical singularities. According to Remark 4.3.2 it consists of the terminal singularities of X together
with the singularities that we obtain when we partially resolve the singularities of type III: recall that  

* 1/9(1,4,7) contributes with three terminal singularities of type 1/3(1,1,2) and 

* 1/14(1,9,11) with seven nodes 

(cf. Lemma 4.2.16 and Lemma 4.2.17). */

Xhat_Basket:=function(BaskTer,BaskTypeIII)
BaskHat:=BaskTer;
for S in BaskTypeIII do 			     				
	if S[1] eq 9 then                    				
		BaskHat:= BaskHat join {*(1/3)^^3*};			
	else 
		BaskHat:= BaskHat join {*(1/2)^^7*};
	end if;
end for;
return BaskHat; end function;