/* With this code, we classify all rigid hyperelliptic fourfolds with holonomy group He(3), in particular we use it to perform the 
   computation in the proof of Proposition 6.11. 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: diffeomorphic classification */

/* PART I: generation of all free rigid He(3)-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;
t:=(1+2*ze)/3;

/* According to Remark 3.2, the analytic representation of He(3) is given by: */

g:=Matrix(F, 3, 3,[0,0,1,1,0,0,0,1,0]);
h:=DiagonalMatrix([1,ze^2,ze]);
I3:=IdentityMatrix(F, 3);
k:=g^-1*h^-1*g*h;

/* Next we define the linear and translation parts, of all possible biholomorphisms
   f(z)=Az+d. Recall that A=Diag(c,D) is block diagonal, where c is a unit in the 
   Eisenstein integers and the 3x3 matrix "D" belongs to the group <D1,D2>, see Corollary 6.8. 
   The translation parts, on the other hand, are determined according to Proposition 6.10: 
   "Listd1" contains the possible first-components "d1" of the translation vector d=(d1,d2,d3,d4). 
   These are precisely the elements in the kernel of the multiplication map 3(zeta_3-1): E --> E. 
   The last three components d':=(d2,d3,d4) of d are the elements in the kernel of 
   (zeta_3-1):C^3/Lambda' --> C^3/Lambda'. This kernel depends on the choice of the lattice Lambda_i. 
   For i=1, we denote it by "Listd234Small" and for i=2 by "Listd234Big". */

D1:=Matrix(F, 3, 3,[ze,0,0,0,ze^2,0,0,0,1]);
D2:=(ze+2)/3*Matrix(F, 3, 3,[1,ze^2,ze^2,ze^2,1,ze^2,ze^2,ze^2,1]);
S:=sub<GL(3,F)|D1,D2>;

Listd1:=[a/3 + b*t/3 : a in [0..2], b in [0..8]];
Listd234Small:=[a*Matrix(F, 3,1,[1/3,1/3,1/3])+ b*Matrix(F, 3,1,[t,0,0]) + c*Matrix(F, 3,1,[0,t,0]): a,b,c in [0..2]];
Listd234Big:=[a*Matrix(F, 3,1,[1/3,1/3,1/3])+ b*Matrix(F, 3,1,[t,0,0]) + c*Matrix(F, 3,1,[1/3,2/3,0]): a,b,c in [0..2]];
Listd234:=[Listd234Small,Listd234Big];

/* 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)^n (n=3,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 "InLattSmall" takes as input a complex vector (w1,w2,w3) and checks if (0,w1,w2,w3) belongs to the small lattice Lambda_1, the function "InLattBig" checks 
   if (0,w1,w2,w3) belongs to the big lattice Lambda_2. Recall that each "Lambda_i" is a direct sum of a lattice of rank two and a lattice of rank six. So the function 
   checks if (w1,w2,w3) belongs to the rank six summand. */ 

InLattSmall:=function(w)
ve:=Matrix(F, 3,1,[t,t,t]);
return IntegralTest(w) or IntegralTest(w-ve) or IntegralTest(w-2*ve);
end function;

InLattBig:=function(w)
ve2:=Matrix(F,3,1,[t,-t,0]);
return InLattSmall(w) or InLattSmall(w-ve2) or InLattSmall(w-2*ve2);
end function;

// The function "InLatt" is equal to "InLattSmall", if the parameter j=1 and equal to "InLattBig" if j=2. 

InLatt:=function(w,j)
return [InLattSmall(w),InLattBig(w)][j];
end function;

/* The function "RigidActions" determines all rigid and free He(3)-actions (in standard form) on C^4/Lambda_j, for each j=1,2. 
   The second output is a representative for each of the special cohomology classes in H^1(He(3),C^4/Lambda_j). */

RigidActions:=function(j)
E3:={0,1/3,ze/3,2/3,(ze+1)/3,(ze+2)/3,2*ze/3,2*(ze+1)/3,(2*ze+1)/3};
ListOfActions:={};
a1:=1/3;
b2:=1/3;
for j3 in [0..2] do
for j4 in [0..2] do 
c1:=(1-ze)*a1;
b3:=ze^j3*b2;
b4:=ze^j4*b2;
for a2 in E3 do
for a3 in E3 do
for a4 in E3 do 
v1:=Matrix(F, 3,1,[(1-ze)*b2,(1-ze)*b3,(1-ze)*b4]);
v4:=Matrix(F, 3,1,[b4-ze*b2+(1-ze)*a2,b2-ze*b3,b3-ze*b4+(1-ze^2)*a4]);
v5:=Matrix(F, 3,1,[a2+a3+a4,a2+a3+a4,a2+a3+a4]);
v6:=Matrix(F, 3,1,[(1-ze)*a2,(1-ze)*a3,(1-ze)*a4]);
if InLatt(v1,j) and InLatt(v4,j) and InLatt(v5,j) and InLatt(v6,j) and 
IsIntegral((1-ze)*a1-c1) and IsIntegral((ze-1)*c1) 
and 
	not IsIntegral(ze^2*(a2+a3)+a4+ze^2*(b2+b4)+b3) and 
	not IsIntegral(-ze*(a2+a3)-a4+b2+b3+ze^2*b4) and 
	not IsIntegral(ze*(a2+a3)+a4-ze*(b2+b3)-b4) and 
	not IsIntegral(ze^2*(a2+a3)+a4 +b2+ze*b3+b4) and 
	not IsIntegral((1-ze)*b2) and 
	not IsIntegral((1-ze)*b3) and 
	not IsIntegral((1-ze)*b4)
then 

G:=Matrix(F, 5, 5,
[1,0,0,0,a1,
0,0,0,1,a2,
0,1,0,0,a3,
0,0,1,0, a4,
0,0,0,0,1]);

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]);

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

RefListeEx:=ListOfActions;
SpecialClasses:={};
while not IsEmpty(ListOfActions) do 
	RefListeEx:=ListOfActions;
	W1:=Rep(ListOfActions); Include(~SpecialClasses,W1); 
	d1G:=Submatrix(W1[1],2,5,3,1);
	d1H:=Submatrix(W1[2],2,5,3,1);
	for W2 in RefListeEx do 
		d2G:=Submatrix(W2[1],2,5,3,1);
		d2H:=Submatrix(W2[2],2,5,3,1);
			for d234 in Listd234[j] do 
				if InLatt((g-I3)*d234-d1G+d2G,j) and InLatt((h-I3)*d234-d1H+d2H,j) 
				then 
				for d1 in Listd1 do
					if IsIntegral(W1[1][1][5]-W2[1][1][5]) and IsIntegral((ze-1)*d1)  
					then Exclude(~ListOfActions, W2); break d234 ;
					end if;
				end for;
				end if;
			end for;
	end for; 
end while; 

return LA, SpecialClasses;
end function;

/* PART II: The functions to check the cocycle-condition (2) in Remark 4.6, for the last three coordinates.   
   
   Given an affinity f(x)=Dx+d as a 4x4 matrix "B", the function "PartsAff" returns 
   the translation vector "d" and the 3x3 matrix "D". */

PartsAff:=function(B);
D:=Submatrix(B,1,1,3,3);
d:=Submatrix(B,1,4,3,1);
return D,d;
end function;

/* Given a matrix "D=D_phi" in the group S=<D1,D2>  and "the last three coordinates" of two rigid actions [G1,H1] and [phiG2,phiH2], where the    
   second one is conjugated by the automorphism "phi" corresponding to "D", the function "Existsd234" is used to check the existence of a vector
   (d2,d3,d4) in "Listd234" which fulfills the cocycle-condition (2) of Remark 4.6 for the last three coordinates. If j=1, then we perform the
   computation with the small lattice, for j=2 we use the big lattice instead. Here we use again the fact that each "Lambda_i" is a direct sum 
   of a lattice of rank two and a lattice of rank six. So the condition (2) in Remark 4.6 can be checked separately for the first coordinate and 
   for the last three coordinates. This is the main difference between the code for He(3) and Z_3^2. */ 

Existsd234:=function(G1,H1,phiG2,phiH2,D,j)
test:=false;
I3:=Matrix(F, 3, 3,[1,0,0,0,1,0,0,0,1]);
MPphiG2, TPphiG2:=PartsAff(phiG2); 
MPphiH2, TPphiH2:=PartsAff(phiH2);
MPG1, TPG1:=PartsAff(G1); 
MPH1, TPH1:=PartsAff(H1);
for d in Listd234[j] do 
	if InLatt((MPphiG2-I3)*d-D*TPG1+TPphiG2,j) and InLatt((MPphiH2-I3)*d-D*TPH1+TPphiH2,j) then 
		test:=true; break d;
	end if;
end for; 	
return test;
end function;

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

   Given an affinity f(z)=diag(c,D)z+(d_1,d') as a 5x5 matrix, the function "DBlock" returns the affinity z--> Dz+d' 
   as a 4x4 matrix. */

DBlock:=function(A)
return Submatrix(A,2,2,4,4);
end function;

/* Given a matrix "D=D_phi" in the group S=<D1,D2> the function "IdAuto" returns the corresponding 
   automorphism "phi". More precisely, it returns the exponents [n1,n2,n3] and [m1,m2,m3] such that 
   phi(g)=g^n1*h^n2*k^n3 and phi(h)=g^m1*h^m2*k^m3. */

IdAuto:=function(D,g,h,k)
phig_ex:=[]; phih_ex:=[];
for n1 in [0..2] do
for n2 in [0..2] do
for n3 in [0..2] do
	if D*g*D^(-1) eq g^(n1)*h^(n2)*k^(n3) then
		phig_ex:=[n1,n2,n3];
	end if;
	if D*h*D^(-1) eq g^(n1)*h^(n2)*k^(n3) then 
		phih_ex:=[n1,n2,n3];
	end if;
end for; end for; end for;
return phig_ex, phih_ex;
end function;

/* The function "TestIso" is used to check if two given actions "W1" and "W2" define biholomorphic quotients. 
   As an input we need the group "S=<D1,D2>" and the parameter "j". If "j=1", then we work with the small 
   lattice, if "j=2", then we use the big lattice. */

TestIso:=function(W1,W2,S,j)
G1:=W1[1]; H1:=W1[2];
G2:=W2[1]; H2:=W2[2];
K2:=G2^-1*H2^-1*G2*H2;
test:=false;
for D in S do
	phig_ex, phih_ex:= IdAuto(D,g,h,k);
	phiG2:=G2^(phig_ex[1])*H2^(phig_ex[2])*K2^(phig_ex[3]);
	phiH2:=G2^(phih_ex[1])*H2^(phih_ex[2])*K2^(phih_ex[3]);
	if Existsd234(DBlock(G1),DBlock(H1),DBlock(phiG2),DBlock(phiH2),D,j) then 
	for c in [1,ze,ze^2,-1,-ze,-ze^2] do 
		for d1 in Listd1 do
			if IsIntegral(c/3-phiG2[1][5]) and IsIntegral((phiH2[1][1]-1)*d1+phiH2[1][5]) then 
				test:=true;
				break D;
			end if;
		end for;
	end for;
	end if;
end for;

return test;
end function;

/* Here the magic happens: we compare all actions in both cases "j=1" and "j=2" (small and big lattice) to get a complete list of
   biholomorphism classes. For each lattice, we obtain two classes. 
   On our machine "Intel(R) Core(TM)2 Duo CPU 3.00GHz" the complete classification for both lattices takes "8.820" seconds. 
   The output is directly on the terminal, in contrast to the Z_3^2 case, no file is created. */

ClassHeis:=function(j)

LA, SpecialClasses:=RigidActions(j);
RefListeEx:=SpecialClasses;

IsoClasses:=[];
while not IsEmpty(SpecialClasses) do 
	RefListeEx:=SpecialClasses;
	W1:=Rep(SpecialClasses); Append(~IsoClasses,W1); 
	for W2 in RefListeEx do 
		if TestIso(W1,W2,S,j) then 
		Exclude(~SpecialClasses, W2);
		end if;
	end for; 
end while; 
return IsoClasses;
end function;

// Here we run the classification function for "j=1,2"

for j in [1,2] do
ClassHeis(j);
end for;


/*********************** PART IV: 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;

// The function "RealMat" turns a complex 3x3 matrix "D" into a real 6x6 matrix under the canonical embedding  

RealMat:=function(D)
return Matrix(F, 6, 6,
[RI(D[1][1])[1],-RI(D[1][1])[2],RI(D[1][2])[1],-RI(D[1][2])[2],RI(D[1][3])[1],-RI(D[1][3])[2],   
RI(D[1][1])[2],RI(D[1][1])[1],RI(D[1][2])[2],RI(D[1][2])[1],RI(D[1][3])[2],RI(D[1][3])[1],
RI(D[2][1])[1],-RI(D[2][1])[2],RI(D[2][2])[1],-RI(D[2][2])[2],RI(D[2][3])[1],-RI(D[2][3])[2],   
RI(D[2][1])[2],RI(D[2][1])[1],RI(D[2][2])[2],RI(D[2][2])[1],RI(D[2][3])[2],RI(D[2][3])[1],
RI(D[3][1])[1],-RI(D[3][1])[2],RI(D[3][2])[1],-RI(D[3][2])[2],RI(D[3][3])[1],-RI(D[3][3])[2],   
RI(D[3][1])[2],RI(D[3][1])[1],RI(D[3][2])[2],RI(D[3][2])[1],RI(D[3][3])[2],RI(D[3][3])[1]]);
end function;

// These are the matrices corresponding to the involutions phi_3 and phi_4 in the table of Lemma 6.5.

D3:=Matrix(F, 6, 6,[1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,1,0,0]);
D4:=Matrix(F, 6, 6,[1,0,0,0,0,0,0,-1,0,0,0,0,0,0,1,0,0,0,0,0,0,-1,0,0,0,0,0,0,1,0,0,0,0,0,0,-1]);

// The group of semilinearities SR=<D1,...,D4> as a subgroup of GL(6,F). 

SR:=sub<GL(6,F)|RealMat(D1),RealMat(D2),D3,D4>;

/* Given a C-linear or C-antilinear map "D" as a 6x6 matrix and a complex vector "v" with three components, the function 
   "RealAction" determines the product "w=D*v"and returns it as a complex vector with three components */ 

RealAction:=function(D,v)
vRe:=Matrix(F, 6, 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]]);
wRe:=D*vRe;
w:=Matrix(F, 3, 1,[wRe[1][1]+i*wRe[2][1],wRe[3][1]+i*wRe[4][1],wRe[5][1]+i*wRe[6][1]]);
return w;
end function;

/* The function "Existsd234Real" performs the same computations as "Existsd234", except that "D" belongs to the larger group 
   SR=<D1,...,D4>. */

Existsd234Real:=function(G1,H1,phiG2,phiH2,D,j)
test:=false;
I3:=Matrix(F, 3, 3,[1,0,0,0,1,0,0,0,1]);
MPphiG2, TPphiG2:=PartsAff(phiG2); 
MPphiH2, TPphiH2:=PartsAff(phiH2);
MPG1, TPG1:=PartsAff(G1); 
MPH1, TPH1:=PartsAff(H1);
for d in Listd234[j] do 
	if InLatt((MPphiG2-I3)*d-RealAction(D,TPG1)+TPphiG2,j) and InLatt((MPphiH2-I3)*d-RealAction(D,TPH1)+TPphiH2,j) then 
		test:=true; break d;
	end if;
end for; 	
return test;
end function;

/* The function "TestDiffeo" checks in analogy to "TestIso" if two given actions "W1" and "W2" define diffeomorphic
   quotients. As an input we need the group SR=<D1,...,D4> and the parameter "j". If "j=1", then we use the small lattice, if "j=2" 
   then we use the big lattice. */

TestDiffeo:=function(W1,W2,SR,j)
G1:=W1[1]; H1:=W1[2];
G2:=W2[1]; H2:=W2[2];
K2:=G2^-1*H2^-1*G2*H2;
test:=false;
for D in SR do
	phig, phih:= IdAuto(D,RealMat(g),RealMat(h),RealMat(k));
	phiG2:=G2^(phig[1])*H2^(phig[2])*K2^(phig[3]);
	phiH2:=G2^(phih[1])*H2^(phih[2])*K2^(phih[3]);
	if Existsd234Real(DBlock(G1),DBlock(H1),DBlock(phiG2),DBlock(phiH2),D,j) then 
	for c in [1,ze,ze^2,-1,-ze,-ze^2] do 
		for d1 in Listd1 do
			if IsIntegral(c/3-phiG2[1][5]) and IsIntegral((phiH2[1][1]-1)*d1+phiH2[1][5]) then 
				test:=true;
				break D;
			end if;
		end for;
	end for;
	end if;
end for;

return test;
end function;

// To start the classification, please run the following code. Replace "j=1" by "j=2" for the classification in the case of Lambda_2.

j:=1;
W:=ClassHeis(j);
W;
W1:=W[1]; W2:=W[2];
TestDiffeo(W1,W2,SR,j);


Cputime();