// The Case G=C3^2

// The function "Kernels" is used to determine the orbits of the action of the normalizer N_{Aut(E^3)}(Z_3^2) on the 
// set $\mathcal K$ of possible kernels "K < Z_3^3 = Fix(ze)^3=<t>^3" of the addition map. As an argument it takes the dimension "dim" 
// of the kernels which is between 0 and 2. In the computation we identify t with 1 and work with the finite field F3

F3:=FiniteField(3);
Nor:=sub<GL(3,F3)|DiagonalMatrix([-1,1,1]), Matrix(F3,3,3, [[0,1,0],[1,0,0],[0,0,1]]), Matrix(F3,3,3, [[0,1,0],[0,0,1],[1,0,0]])>;

V:=VectorSpace(F3,3);
UnitVec:=[V![1,0,0],V![0,1,0],V![0,0,1]];

function Kernels(dim)
	List:=[];
	Orb:=OrbitsOfSpaces(Nor,dim);
	for W in Orb do
		test:=true;
		for e in UnitVec do
			if e in W[2] then
				test:=false;
				break e;
			end if;
		end for;
		if test then
			List:=Append(List,W[2]);
		end if;
	end for;
	return List;
end function;


// There are precisely 4 orbist: 1 of dim 0, 2 of dim 1, 1 of dim 2

/*
#Kernels(0);
#Kernels(1);
#Kernels(2);
*/


/* ********* Below is the main part of our code ****************************

With this code, we classify all 3-dimensional quotients of abelian varieties by an action of Z_3^2 that has a non-empty and finite fixed locus and preserves the volume form.
In particular, we use it to perform the computations in the proof of Proposition 4.26.
The code consits of three parts:
 	PART I: generation of all Z_3^2-actions with isolated fixed points and the initialization of all potential biholomorphisms
	PART II: the functions to check the cocycle-condition in Remark 3.7(b)
	PART III: the classification functions for the biholomorphic classification */

/***************************************************************************/

/* PART I: generation of all Z_3^2-actions with isolated fixed points and the initialization of all potential biholomorphisms. */

//We work with the 3rd cyclotomic field, ze is a 3rd primitive root of unity.

F<ze>:=CyclotomicField(3);

/*  We choose generators h and k of Z_3^2. According to Theorem 2.3, the analytic representation is given by:*/

h:=DiagonalMatrix([1,ze^2,ze]);
k:=DiagonalMatrix([ze,ze,ze]);


//E3 denotes the group of 3-torsion points of the Fermat elliptic curve E.
E3:={1/3*(a+b*ze): a, b in {0,1,2}};

// Generator of the fixed locus of ze in E:
t:=1/3+2/3*ze;

I3:=DiagonalMatrix([1,1,1]);

// These are representatives for the 4 orbits, i.e., the kernels K_i of the addition map, that we will use.

K1:={Matrix(F,3,1,[0,0,0])};
K2:={a*Matrix(F,3,1,[t,t,0]): a in {0,1,-1}};
K3:={a*Matrix(F,3,1,[t,t,t]): a in {0,1,-1}};
K4:={a*Matrix(F,3,1,[t,t,t])+b*Matrix(F,3,1,[t,-t,0]): a,b in {0,1,-1}};

Kernels:=[K1,K2,K3,K4];

/* The function "IntegralTest" checks if the entries of a 3-vector v (given as 3x1-matrix) are Eisenstein integers, i.e. integral over Z. */

function IntegralTest(v)
	return IsIntegral(v[1][1]) and IsIntegral(v[2][1]) and IsIntegral(v[3][1]);
end function;

/* The function "InLatt" takes as input a vector "v" and a parameter "j=1,...,4" and decides whether "v" belongs to the lattice "Lambda_j=Z[ze]^3+K_j". */

function InLatt(v,j)
	K:=Kernels[j];
	for l in K do
		if IntegralTest(v-l) then
			return true;
		end if;
	end for;
	return false; 
end function;

/* The function "GoodCond" checks whether a cocycle in standard form is good (cf. Lemma 4.21). */

function GoodCond(a1,a2,a3,j)
	for t1 in [0,t,-t] do
		for t2 in [0,t,-t] do 
			w1:=Matrix(F,3,1,[a1,t1,t2]);
			w2:=Matrix(F,3,1,[t1,a2,t2]);
			w3:=Matrix(F,3,1,[t1,t2,a3]);
			if InLatt(w1,j) or InLatt(w2,j) or InLatt(w3,j) then 
				return false;				
			end if; 
		end for; 
	end for;
	return true;
end function;


/* Here, we define the possible translation parts of potential biholomorphisms between two quotients, as explained in Proposition 4.25. Note that all candidates are contained in the set E[9]^3. We check with the code which of these elements are in fact fiexed points of multiplication with \ze_3 and then, we choose one representative of each class of the elements in C^3 mod the lattice.
   We also use them to define the 1-coboundaries which we need to compute a set of represenatatives for the good cohomology classes.
   The input is the parameter "j" that refers to the lattice, resp. the kernel, that is used in the computation. */

function Listd(j)
	ListdMax:={a/9 + b/9*ze: a in [0..8], b in [0..8]};
	listd:=[];
	for d1 in ListdMax do
	for d2 in ListdMax do
	for d3 in ListdMax do
		d:=Matrix(F,3,1,[d1,d2,d3]);
		if  InLatt((ze-1)*I3*d,j) then
			Append(~listd,d);
		end if;
	end for; end for; end for;
	Reflistd:=listd;
	shortlistd:=[];
	while not IsEmpty(listd) do
		Reflistd:=listd;
		d:=Rep(listd); Append(~shortlistd,d);
		for e in Reflistd do
			if InLatt(d-e,j) then
				Exclude(~listd,e);
			end if;
		end for;
	end while;
	return shortlistd;
end function;

/* With the function "TestCohom" we test if two given cocycles give the same cohomology class in H^1(Z_3^2,E^3/K_j).
   Input: "v1", "v2": two translation vectors; "j": parameter for the kernel; "Ld": list of translation parts of the biholomorphisms computed with "Listd(j)" */

function TestCohom(v1,v2,j,Ld)
	for d in Ld do
		if InLatt(v1-v2-(h*d-d),j) then
			return true;
		end if;
	end for;
	return false;
end function;

/* The function "Actions" determines all actions with only isolated fixed points with translation part in standard form on E^3/K_j, for each j=1,...,4 (first output).
   The second output is a list consisting of one representative for each good comohomology class in H^1(Z_3^2,E^3/K_j).
   As an input, it takes the parameter "j" corresponding to the kernel K_j and the list of translation vectors "Ld=Listd(j)".*/

function Actions(j,Ld)
	ListTransVec:={};
	ListOfActions:={};
	for a1 in E3  do 
	for a2 in E3  do
	for a3 in E3  do
		v:=Matrix(F,3,1,[(ze-1)*a1,(ze-1)*a2,(ze-1)*a3]);
		if InLatt(v,j) and GoodCond(a1,a2,a3,j) then
			ListTransVec:=Include(ListTransVec, Matrix(F,3,1,[a1,a2,a3]));
			ListOfActions:=Include(ListOfActions, Matrix(F,4,4, [[1,0,0,a1],[0,ze^2,0,a2],[0,0,ze,a3],[0,0,0,1]]));
		end if;
	end for; end for; end for;
	GoodClasses:={};
	RefListTransVec:=ListTransVec;
	while not IsEmpty(ListTransVec) do
		RefListTransVec:=ListTransVec;
		v1:=Rep(ListTransVec);
		GoodClasses:=Include(GoodClasses,Matrix(F,4,4, [[1,0,0,v1[1][1]],[0,ze^2,0,v1[2][1]],[0,0,ze,v1[3][1]],[0,0,0,1]]));
		for v2 in RefListTransVec do
			if TestCohom(v1,v2,j,Ld) then
				Exclude(~ListTransVec,v2);
			end if;
		end for;
	end while;
	return ListOfActions, GoodClasses;
end function;


/* The function "Normal" determines the normalizers N_C(Lambda_j) from Proposition 4.24 */

/* First, we define the normalizer N_{Aut(E^3)}(rho(Z_3^2)):*/
Nor:=MatrixGroup<3,F| DiagonalMatrix([-ze,1,1]), Matrix([[0,1,0],[1,0,0],[0,0,1]]), Matrix([[0,1,0],[0,0,1],[1,0,0]])>; 

function Normal(j)
	List:={};
	K:=Kernels[j];
	for C in Nor do
		test:= true;
		for l in K do
			if not InLatt(C*l,j) or not InLatt(C^-1*l,j) then
				test:=false; break l;
			end if;
		end for;
		if test then
			List:=Include(List,C);
		end if;
	end for;
	return List;
end function;

/*******************************************************************************/

/* PART II: The functions to check the cocycle-condition in Remark 3.7(b)*/

/* Given an affinity "f(z)=Cz + d" as a 4x4 matrix "B", the function "PartsAff" returns the 3x3 matrix "C" and the translation vector "d". */

function PartsAff(B)
	C:=Submatrix(B,[1,2,3],[1,2,3]);
	d:=Submatrix(B,[1,2,3],[4]);
	return C,d;
end function;


/* Given a matrix "C=C_{phi}" in the normalizer N_C(Lambda_j) and two actions "H1" and "phiH2", where the second one is composed with the automorphism "phi", the function "Existsd" checks the existence of a vector d in "Ld=Listd(j)" fulfillling the cocylcle-condition.*/

function Existsd(C,H1,phiH2,j,Ld)
	MPH1, TPH1:=PartsAff(H1);
	MPphiH2,TPphiH2:=PartsAff(phiH2);
	for d in Ld do
		if InLatt((MPphiH2-I3)*d-C*TPH1+TPphiH2,j) then
			return true;
		end if;
	end for;
	return false;
end function;

/******************************************************************************/

/* PART III: the classification functions for the biholomorphic classification. */

/* Given a matrix C in the normalizer N_{Aut(E^3)}(rho(Z_3^2)), the function "IdAuto" returns the corresponding automorphism. More precisley, it returns the image of h as a tuple [a,b] s.t. phi(h)=h^a*k^b. (Remark that k is fixed by phi).
 The possible images of h in terms of [a,b] s.t. phi(h)=h^a*k^b are listed at the beginning of the function. */

function IdAuto(C)
	ListImages:=[[1,0],[2,0],[1,1],[2,2],[1,2],[2,1]];
	for n in ListImages do
		if C*h*C^-1 eq h^n[1]*k^n[2] then
			return n;
		end if;
	end for;
end function;

/* The function "TestIso" checks if two given actions "W1" and "W2" lead to biholomorphic quotients.
   As an input it gets the actions "W1" and "W2", the paramter "j", the normalizer "N=N_C(Lambda_j)", and the list of translation vectors "Ld=Listd(j)". */

function TestIso(W1,W2,j,N,Ld)
	H1:=W1; 
	H2:=W2;
	for C in N do
		phi:=IdAuto(C);
		phiH2:=H2^phi[1]*Matrix(F,4,4,[[ze,0,0,0],[0,ze,0,0],[0,0,ze,0],[0,0,0,1]])^phi[2];
		if Existsd(C,H1,phiH2,j,Ld) then
			return true;
		end if;
	end for;
	return false;
end function;

/* The function "ClassZ3xZ3" is the main classification function. It determines for each lattice all biholomorphism classes of 
   three-dimensional quotients of abelian varities by an action of Z_3^2 that has finite and non-empty fixed locus and preserves the volume form. Running this function finishes the proof of Proposition 4.26.
   Because of the size of the output and the length of the computation a file "Z3xZ3.txt" is created. It contains for each "j" 
   1) the number of actions, 
   2) the number of  good cohomology classes, 
   3) the size of the normalizer N_C(Lambda_j), 
   4) the number of biholomorphism classes and 
   5) for each biholomorphism class a corresponding action on E^3/K_j. */

function ClassZ3xZ3(j)
	F:="Z3xZ3.txt";
	fprintf F, "Kernel %o)\n \n", j;
	IsoClasses:=[];
	Ld:=Listd(j);
	LA,GoodClasses:=Actions(j,Ld);
	fprintf F, "Number of actions with isolated fixed points: %o \n \n", #LA;
	fprintf F, "Number of good cohomology classes: %o\n\n", #GoodClasses;
	RefListAct:=GoodClasses;
	N:=Normal(j);
	fprintf F, "Size of the normalizer: %o \n \n", #N;
	while not IsEmpty(GoodClasses) do
		RefListAct:=GoodClasses;
		W1:=Rep(GoodClasses);
		Append(~IsoClasses,W1);
		for W2 in RefListAct do
			if TestIso(W1,W2,j,N,Ld) then
				Exclude(~GoodClasses,W2);
			end if;
		end for;
	end while;
	fprintf F, "Number of biholomorphism classes: %o \n \n", #IsoClasses;
	fprintf F, "Actions [H]: \n %o \n \n \n \n", IsoClasses;
	return "Classification for kernel", j, "is completed!";
end function;

/* Here, we run the classification function for each "j=1,...,4". */

for j in [1..4] do
	ClassZ3xZ3(j);
end for;