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

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

Kernels:=function(dim)
List:=[];
Orb:=OrbitsOfSpaces(Nor,dim);
V:=VectorSpace(F3,4);
UnitVec:=[V![1,0,0,0],V![0,1,0,0],V![0,0,1,0], V![0,0,0,1]];
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 eq true then 
List:=Append(List,W[2]);
end if;
end for;
return List;
end function;

// There are precisely 13 orbits: 1 of dim 0, 5 of dim 1, 6 of dim 2 and 1 of dim 3

#Kernels(0);
#Kernels(1);
#Kernels(2);
#Kernels(3);

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

   With this code, we classify all rigid hyperelliptic fourfolds with holonomy group Z_3^2, in particular we use it to perform the 
   computations in the proofs of Proposition 6.14 and 6.18. The code consists of four parts:

   PART I: generation of all free rigid He(3)-actions and the initialisation of all potential biholomorphisms
   PART II: the functions to check the cocycle-condition (2) in Remark 4.6  for the last three coordinates 
   PART III: the classification functions for the biholomorphic classification
   PART IV: verification of Proposition 6.18 i.e. the diffeomorphic classification */



/* PART I: generation of all free and rigid Z_3^2-actions and the initialisation of all potential biholomorphisms

   We work with the 12-th cyclotomic field. In this case ze:=F.1^4 is the first primitive third root of unity and i=F.1^3 is the 
   imaginary unit.  */

F:=CyclotomicField(12);
ze:=F.1^4;
i:=F.1^3;
E3:=[0,(2*ze+1)/3,(ze+2)/3,1/3,ze/3,2/3,(ze+1)/3,2*ze/3,2*(ze+1)/3];
t:=(1+2*ze)/3;

// These are some "nice" representatives for the 13 orbits i.e. the kernels K_i of the addition map that we will use

K1:={Matrix(F,4,1,[0,0,0,0])};
K2:={a*Matrix(F,4,1,[0,0,t,t]): a in {0,1,2}};
K3:={a*Matrix(F,4,1,[0,t,t,t]): a in {0,1,2}};
K4:={a*Matrix(F,4,1,[t,0,0,t]): a in {0,1,2}};
K5:={a*Matrix(F,4,1,[t,0,t,t]): a in {0,1,2}};
K6:={a*Matrix(F,4,1,[t,t,t,t]): a in {0,1,2}};
K7:={a*Matrix(F,4,1,[0,t,t,t])+ b*Matrix(F,4,1,[0,t,-t,0]): a,b in {0,1,2}};
K8:={a*Matrix(F,4,1,[0,0,t,t])+ b*Matrix(F,4,1,[t,0,-t,0]): a,b in {0,1,2}};
K9:={a*Matrix(F,4,1,[0,0,t,t])+ b*Matrix(F,4,1,[t,t,0,0]): a,b in {0,1,2}};
K10:={a*Matrix(F,4,1,[0,0,t,t])+ b*Matrix(F,4,1,[t,t,0,t]): a,b in {0,1,2}};
K11:={a*Matrix(F,4,1,[t,0,0,t])+ b*Matrix(F,4,1,[t,t,t,-t]): a,b in {0,1,2}};
Kexc:={a*Matrix(F,4,1,[t,-t,0,t])+ b*Matrix(F,4,1,[t,t,-t,0]): a,b in {0,1,2}};
K12:={a*Matrix(F,4,1,[-t,t,0,0])+ b*Matrix(F,4,1,[t,0,t,t])+ c*Matrix(F,4,1,[t,t,t,0]): a,b,c in {0,1,2}};

/* The function "IntegralTest", checks if all the entries of a vector "w" are cyclotomic 
   integers i.e. integral over Z. 
   Applied to elements in Q(ze)^4 we use it to check if the entries are Eisenstein. */ 

IntegralTest:=function(w)
test:=true;
for l in [1..NumberOfRows(w)] do
if not IsIntegral(w[l][1]) then 
test:=false; break l;
end if;
end for;
return test;
end function;

/* The function "InLatt" takes as input a vector "w" and a parameter j=1,...,13 and checks if "w" belongs to the lattice 
   "Lambda_j=Z[ze]^4+K_j" (here K_13=K_exc). */

InLatt:=function(v,j)
Kernels:=[K1,K2,K3,K4,K5,K6,K7,K8,K9,K10,K11,K12,Kexc];
K:=Kernels[j];
test:=false;
for l in K do
	if IntegralTest(v-l) then
			test:=true;
			break l;
	end if;
end for;
	return test; 
end function;

/* The function "FreenessCheck" is used to verify the freeness of a rigid action of Z_3^2 in normal form, as
   explained in Lemma 5.19. */

FreenessCheck:=function(c1,b2,b3,b4,j)
Kernels:=[K1,K2,K3,K4,K5,K6,K7,K8,K9,K10,K11,K12,Kexc];
K:=Kernels[j];
test:=true;
for k in K do
	if IsIntegral(c1-k[1][1]) or IsIntegral(b2-k[2][1]) or IsIntegral(b3-k[3][1]) or IsIntegral(b4-k[4][1]) then 
		test:=false; break k;
	end if;
end for;
return test;
end function;

/* Here we define the possible translation parts of potential biholomorphisms between two rigid hyperelliptic fourfolds, as explained in 
   Proposition 6.13. We also use them to define the 1-coboundaries, which we need to compute a set of representatives for the special cohomology classes. 
   The input is the parameter "j", which refers to the lattice, resp. the kernel, that is used in the computation. */
   
Listd:=function(j)
Ld234:={};Repres:={};
Ld:={};
Kernels:=[K1,K2,K3,K4,K5,K6,K7,K8,K9,K10,K11,K12,Kexc];
K:=Kernels[j];

for d2 in E3 do
	for d3 in E3 do
		for d4 in E3 do
			d:=Matrix(F, 4, 1,[0,d2,d3,d4]);
				v:=(ze-1)*d;
				for k in K do 
				w:=Matrix(F, 3, 1,[k[2][1]-v[2][1],k[3][1]-v[3][1],k[4][1]-v[4][1]]);								
				if IntegralTest(w) then 
					Ld234:=Include(Ld234,d);
					break k;
				end if;
				end for;

		end for;
	end for;
end for;

Listd1:=[a/3 + b*t/3 : a in [0..2], b in [0..8]];

for d1 in Listd1 do
	for d234 in Ld234 do 
		d:=Matrix(F, 4, 1,[d1,d234[2][1],d234[3][1],d234[4][1]]);
		Ld:=Include(Ld,d);
	end for;
end for;

refer:=Ld;
while not IsEmpty(Ld) do 
	refer:=Ld;
	d:=Rep(Ld); Include(~Repres,d); 
	for v in refer do 
		if InLatt(d-v,j) then 
		Exclude(~Ld, v);
		end if;
	end for; 
end while; 
return Repres;
end function;

/* The function "RigidActions" determines all rigid and free actions (in standard form) on E^4/K_j, for each j=1,..,13. 
   The second output is a representative for each of the special cohomology classes in H^1(Z_3^2,E^4/K_j). 
   As an input it takes "j" and the list of translation vectors Ld computed with "Listd". */

RigidActions:=function(j,Ld)
ListOfActions:={};

for b2 in E3 do
for b3 in E3 do
for b4 in E3 do 
for c1 in E3 do

v:=Matrix(F,4,1,[(1-ze)*c1, (ze-1)*b2, (ze-1)*b3, (ze-1)*b4]);

if InLatt(v,j) and FreenessCheck(c1,b2,b3,b4,j) then 

H:=Matrix(F, 5, 5,
[ze,0,0,0,0,
0,1,0,0,b2,
0,0,ze^2,0,b3,
0,0,0,ze, b4,
0,0,0,0,1]);

K:=Matrix(F, 5, 5,
[1,0,0,0,c1,
0,ze,0,0,0,
0,0,ze,0,0,
0,0,0,ze,0,
0,0,0,0,1]);

ListOfActions:=Include(ListOfActions,[H,K]);
end if;
end for; end for; end for; end for;
LA:=ListOfActions;

h:=DiagonalMatrix([ze,1,ze^2,ze]);
k:=DiagonalMatrix([1,ze,ze,ze]);
I4:=IdentityMatrix(F, 4);

RefListeEx:=ListOfActions;
SpecialClasses:={};
while not IsEmpty(ListOfActions) do 
	RefListeEx:=ListOfActions;
	W1:=Rep(ListOfActions); Include(~SpecialClasses,W1); 
	d1H:=Submatrix(W1[1],1,5,4,1);
	d1K:=Submatrix(W1[2],1,5,4,1);
	for W2 in RefListeEx do 
		d2H:=Submatrix(W2[1],1,5,4,1);
		d2K:=Submatrix(W2[2],1,5,4,1);
			for d in Ld do 
				if InLatt((h-I4)*d-d1H+d2H,j) and InLatt((k-I4)*d-d1K+d2K,j) then 
					Exclude(~ListOfActions, W2); break d;
				end if;
			end for;
	end for; 
end while; 
return LA, SpecialClasses;
end function;

// With the function "LattIso" we can check if a 4x4-matrix "A" defines an automorphism of the lattice Lambda_j

LattIso:=function(j,A)
GeneratorsKernels:=
[[Matrix(F,4,1,[0,0,0,0])], 
[Matrix(F,4,1,[0,0,t,t])],
[Matrix(F,4,1,[0,t,t,t])],
[Matrix(F,4,1,[t,0,0,t])],
[Matrix(F,4,1,[t,0,t,t])],
[Matrix(F,4,1,[t,t,t,t])],
[Matrix(F,4,1,[0,t,t,t]), Matrix(F,4,1,[0,t,-t,0])],
[Matrix(F,4,1,[0,0,t,t]), Matrix(F,4,1,[t,0,-t,0])],
[Matrix(F,4,1,[t,t,0,0]), Matrix(F,4,1,[0,0,t,t])],
[Matrix(F,4,1,[0,0,t,t]), Matrix(F,4,1,[t,t,0,t])],
[Matrix(F,4,1,[t,0,0,t]), Matrix(F,4,1,[t,t,t,-t])],
[Matrix(F,4,1,[-t,t,0,0]), Matrix(F,4,1,[t,0,t,t]),Matrix(F,4,1,[t,t,t,0])],
[Matrix(F,4,1,[t,-t,0,t]), Matrix(F,4,1,[t,t,-t,0])]];
test:=true;
for ve in GeneratorsKernels[j] do
if  InLatt(A*ve,j) and InLatt(A^-1*ve,j) 
then  test:=true;
else test:=false; break ve;
end if;
end for;
return test;
end function;

// The function "Normal" determines the normalizers N_C(Lambda_j) from Proposition 6.12.

Normal:=function(j)
A1:=DiagonalMatrix([-ze,1,1,1]);
A2:=DiagonalMatrix([1,-ze,1,1]);
A3:=Matrix(F, 4, 4,[1,0,0,0,0,0,1,0,0,0,0,1,0,1,0,0]);
Nor:=sub<GL(4,F)|A1,A2,A3>;                             // this is the normalizer N_{Aut(E^4)}(Z_3^2) (Cor. 5.18)
List:=[];
	for A in Nor do
		if  LattIso(j,A) then 
			List:=Append(List,A);
		end if; 
	end for;
return List;
end function;

/* PART II: The functions to check the cocycle-condition (2) in Remark 4.6  
   
   Given an affinity "alpha(x)=Ax+d" as a 5x5 matrix "B", the function "PartsAff" returns the the 4x4 matrix "A" and the 
   translation vector "d". */

PartsAff:=function(B);
A:=Submatrix(B,1,1,4,4);
d:=Submatrix(B,1,5,4,1);
return A,d;
end function;

/* Given a matrix "A=A_{phi}" in the normalizer N_C(Lambda_j) and two rigid actions [G1,H1] and [phiG2,phiH2], where the second one is 
   conjugated by the automorphism "phi", the function "Existsd" is used to check the existence of a vector
   d=(d1,d2,d3,d4) in "Ld=Listd(j)" which fulfills the cocycle-condition (2) of Remark 4.6. */ 

Existsd:=function(A,H1,K1,phiH2,phiK2,j,Ld)
test:=false;
I4:=IdentityMatrix(F,4);
MPH1,TPH1:=PartsAff(H1);
MPK1,TPK1:=PartsAff(K1);
MPphiH2,TPphiH2:=PartsAff(phiH2);
MPphiK2,TPphiK2:=PartsAff(phiK2);
for d in Ld do
if InLatt((MPphiH2-I4)*d-A*TPH1+TPphiH2,j) and 
	InLatt((MPphiK2-I4)*d-A*TPK1+TPphiK2,j) then 
	test:=true; break d;
end if;
end for;
return test;
end function;


/* PART III: the classification functions for the biholomorphic classification. 
   
   Each matrix "A" in the normalizer N_{Aut(E^4)}(Z_3^2) determines an automorphism phi of Z_3^2 such that phi(k)=k 
   and phi(h)=hk^m. The function "Power" returns the integer "m". */

Power:=function(A)
m:=0;
H:=DiagonalMatrix([ze,1,ze^2,ze]);
K:=DiagonalMatrix([1,ze,ze,ze]);
Hnew:=A*H*A^-1;
for l in [0..2] do
	if Hnew eq H*K^l 
		then m:=l;
	end if;
end for; 
return m;
end function;

/* The function "TestIso" is used to check if two given actions "W1" and "W2" yield biholomorphic quotients. As an input we need 
   the normalizer N=N_C(Lambda_j) and the list "Ld" of translation vectors computed with the function "Listd" */ 

TestIso:=function(W1,W2,j,N,Ld)
test:=false;
H1:=W1[1]; K1:=W1[2];
H2:=W2[1]; K2:=W2[2];
for A in N do
m:=Power(A);
	phiH2:=H2*K2^m;
	phiK2:=K2;
		if Existsd(A,H1,K1,phiH2,phiK2,j,Ld) then 
				test:=true;
				break A;
		end if;
end for;
return test;
end function;

/* The function "ClassZ3xZ3" is the main classification function. It determines for each lattice all biholomorphism classes of rigid 
   hyperelliptic fourfolds with holonomy Z_3^2. On our machine "Intel(R) Core(TM)2 Duo CPU 3.00GHz" the complete classification for all 13 kernels 
   takes about "54" minutes. 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 rigid and free actions, 
   2) the number of special 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^4/Z_3^2. */

ClassZ3xZ3:=function(j);
F:="Z3xZ3.txt";
fprintf F, "Kernel %o)\n \n", j;
IsoClasses:=[];
Ld:=Listd(j);
LA, SpecialClasses:=RigidActions(j,Ld);
fprintf F, "Number of rigid and free actions:  %o \n \n", #LA;
fprintf F, "Number of special cohomology classes:  %o \n \n", #SpecialClasses;
RefListeEx:=SpecialClasses;
	N:=Normal(j);
	fprintf F, "Size of the normalizer:  %o \n \n", #N;
		while not IsEmpty(SpecialClasses) do 
			RefListeEx:=SpecialClasses;
			W1:=Rep(SpecialClasses); Append(~IsoClasses,W1); 
			for W2 in RefListeEx do 
				if TestIso(W1,W2,j,N,Ld) then 
					Exclude(~SpecialClasses, W2);
				end if;
			end for; 
		end while;
fprintf F, "Number of biholomorphism classes:  %o \n \n", #IsoClasses;
fprintf F, "Actions [H,K]: \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,...,13". 

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



/* PART IV: verification of Proposition 6.18 i.e. the diffeomorphic classification 

   The function RI returns real and imaginary part of a complex number "c". */

RI:=function(c)
re:=(c+ComplexConjugate(c))/2;
im:=-i*(c-re);
return [re, im];
end function;

// Next we define the matrices M1, M2 and M3 that generate the normalizer "N=<M1,M2,M3>" in Remark 6.16. 

M1:=Matrix(F, 8, 8,
[1,0,0,0,0,0,0,0,
0,-1,0,0,0,0,0,0,
0,0,1,0,0,0,0,0,
0,0,0,1,0,0,0,0,
0,0,0,0,0,0,1,0,
0,0,0,0,0,0,0,1,
0,0,0,0,1,0,0,0,
0,0,0,0,0,1,0,0]);

M2:=Matrix(F, 8, 8,
[0,0,1,0,0,0,0,0,
0,0,0,1,0,0,0,0,
0,0,0,0,1,0,0,0,
0,0,0,0,0,-1,0,0,
0,0,0,0,0,0,1,0,
0,0,0,0,0,0,0,1,
1,0,0,0,0,0,0,0,
0,1,0,0,0,0,0,0]);

Rze:=RI(ze)[1];
Ize:=RI(ze)[2];

M3:=Matrix(F, 8, 8,
[-Rze,Ize,0,0,0,0,0,0,
-Ize,-Rze,0,0,0,0,0,0,
0,0,1,0,0,0,0,0,
0,0,0,1,0,0,0,0,
0,0,0,0,1,0,0,0,
0,0,0,0,0,1,0,0,
0,0,0,0,0,0,1,0,
0,0,0,0,0,0,0,1]);

N:=sub<GL(8,F)|M1, M2, M3>;

/* Given a semilinear map "A" as a 8x8 matrix and a complex vector "v" with four components, the function 
   "RealAction" determines the product "w=A*v" and returns it as a complex vector */ 

RealAction:=function(A,v)
vRe:=Matrix(F, 8, 1,[RI(v[1][1])[1],RI(v[1][1])[2],RI(v[2][1])[1],RI(v[2][1])[2],RI(v[3][1])[1],RI(v[3][1])[2],
RI(v[4][1])[1],RI(v[4][1])[2]]);
wRe:=A*vRe;
w:=Matrix(F, 4, 1,[wRe[1][1]+i*wRe[2][1],wRe[3][1]+i*wRe[4][1],wRe[5][1]+i*wRe[6][1],wRe[7][1]+i*wRe[8][1]]);
return w;
end function;

// The function "LattIsojk" is used to check if a matrix "A" in the normalizer "N=<M1,M2,M3>" defines an isomorphism between the lattices Lambda_j and Lambda_k
// for j,k=1..12

LattIsojk:=function(A,j,k)
GeneratorsKernels:=
[[Matrix(F,4,1,[0,0,0,0])], 
[Matrix(F,4,1,[0,0,t,t])],
[Matrix(F,4,1,[0,t,t,t])],
[Matrix(F,4,1,[t,0,0,t])],
[Matrix(F,4,1,[t,0,t,t])],
[Matrix(F,4,1,[t,t,t,t])],
[Matrix(F,4,1,[0,t,t,t]), Matrix(F,4,1,[0,t,-t,0])],
[Matrix(F,4,1,[0,0,t,t]), Matrix(F,4,1,[t,0,-t,0])],
[Matrix(F,4,1,[t,t,0,0]), Matrix(F,4,1,[0,0,t,t])],
[Matrix(F,4,1,[0,0,t,t]), Matrix(F,4,1,[t,t,0,t])],
[Matrix(F,4,1,[t,0,0,t]), Matrix(F,4,1,[t,t,t,-t])],
[Matrix(F,4,1,[-t,t,0,0]), Matrix(F,4,1,[t,0,t,t]),Matrix(F,4,1,[t,t,t,0])]];
test:=true;
for vj in GeneratorsKernels[j] do
	for vk in GeneratorsKernels[k] do
		if InLatt(RealAction(A,vj),k) and	
			InLatt(RealAction(A^-1,vk),j) 
 				then  test:=true;
				else test:=false; break vj;
		end if;
	end for;
end for;
return test;
end function;

/* We run over all lattices Lambda_j and Lambda_k and check if there are isomorphisms. The indices j and k as well as the number of 
   orientation preserving and the number   of orientation reversing linear diffeomorphisms are returned. This computation produces 
   the table in Proposition 6.18. */

for j in [1..12] do
for k in [j+1..12] do
	SetAOP:={};SetAOR:={};
		for A in N do
			if LattIsojk(A,j,k) then 
				det:=Determinant(A);
				if det eq  1 then 
				SetAOP:=Include(SetAOP,A);
					else 
				SetAOR:=Include(SetAOR,A);
				end if; 
			end if;
		end for;
	[j,k,#SetAOP,#SetAOR];
end for;
end for;