/* ****************** Cohomology ********************* */

Z := Integers();
M := GL(9,Z);
A := RSpace(Z,9);
G1 := W!SigmaAction;
G2 := W!TauAction;
G3 := W!Phi1Action;
G4 := W!Phi2Action;

/* The parameter WhichField below is used to specify if the ground field
over which we are working contains a cube root of 2 (which we'll call s) 
and/or a primitive sixth root of unity (which we'll call zeta):

Value of WhichField     Property
         1              Contains neither s nor zeta 
         2              Contains s but not zeta
         3              Contains zeta but not s
         4              Contains s and zeta
         
The Generic group of a diagonal dP1 changes with the value of WhichField
as indicated below. The user is invited to change the value as he/she pleases. */

WhichField := 3;

if WhichField eq 1 then
	GenericGroup := sub<W | G1,G2,G3,G4>;
elif WhichField eq 2 then
	GenericGroup := sub<W | G2,G3,G4>;
elif WhichField eq 3 then
	GenericGroup := sub<W | G1,G3,G4>;
else
	GenericGroup := sub<W | G3,G4>;
end if;
	
Q, quotmap := quo<GenericGroup | G3, G4>;

/* Let G be the GenericGroup, and let G0 be a subgroup of G. The composite map
G_0 --->  G  --->  G/<G3,G4>
must be surjective in order for G0 to correspond to the Galois group of a diagnoal
dP1 (see the paper). The variable QuotientSize records the cardinality of G/<G3,G4>
as the parameter WhichField changes. */

if WhichField eq 1 then
	QuotientSize := 6;
elif WhichField eq 2 then
	QuotientSize := 2;
elif WhichField eq 3 then
	QuotientSize := 3;
else
	QuotientSize := 1;
end if;

GM := GModule(GenericGroup);
CM := CohomologyModule(GenericGroup,GM);

SubGpLat := Subgroups(GenericGroup);
print "The generic group has size", #GenericGroup;
print "It has", #SubGpLat, "conjugacy classes of subgroups."; 
// print GenericGroup;

/* We enumerate the elements of GenericGroup and store them in an array. */

ListOfGpElements := [];
for i in [0..Order(G1)-1] do
        for j in [0..Order(G2)-1] do
                for k in [0..Order(G3)-1] do
                        for l in [0..Order(G4)-1] do
                                elt := G1^i*G2^j*G3^k*G4^l;
                                ListOfGpElements := Append(ListOfGpElements,elt);
                        end for;
                end for;
        end for;
end for;

/* The following function takes an integer n as input and outputs the vector 
[i,j,k,l] such that G1^i*G2^j*G3^k*G4^l is the n-th element of ListOfGpElements */

ConvertToFactorization := function(n)
	i := Truncate((n-1)/72);
	Garbage, calculatej := Quotrem(n-1,72);
	j := Truncate(calculatej/36);
	Garbage, calculatek := Quotrem(calculatej,36);
	k := Truncate(calculatek/6);
	Garbage, calculatel := Quotrem(calculatek,6);
	Garbage, l := Quotrem(n-1,6);
	return [i,j,k,l];
end function;

/* The following function takes for input that i-th subgroup of a lattices and 
re-writes its generators using the previous two functions. */

ComputeGpGens := function(i, Lattice)
	Hi := Lattice[i]`subgroup;
	gensbad := SetToSequence(Generators(Hi));
	gensgood := [];
	for i in [1..#gensbad] do
		for j in [1..#ListOfGpElements] do
			if gensbad[i] eq ListOfGpElements[j] then
				exponent := ConvertToFactorization(j);
				// StringGen := "G1^",exponent[1],"*G2^",exponent[2],"*G3^",
				exponent[3],"*G4^",exponent[4];
				Append(~gensgood,exponent);
			end if;
		end for;
	end for;
	return gensgood;
end function;

writeGroup := function(CohoGpi)
	if CohoGpi eq [ 2 ] then
		return "Z/2Z";
	elif CohoGpi eq [ 3 ] then
		return "Z/3Z";
	elif CohoGpi eq [ 2, 2 ] then
		return "Z/2Z + Z/2Z";
        elif CohoGpi eq [ 2, 6 ] then
                return "Z/2Z + Z/6Z";
	elif CohoGpi eq [ 3, 3 ] then
		return "Z/3Z + Z/3Z";
        elif CohoGpi eq [ 3, 3, 3 ] then
                return "Z/3Z + Z/3Z + Z/3Z";
        elif CohoGpi eq [ 3, 3, 3, 3 ] then
                return "Z/3Z + Z/3Z + Z/3Z + Z/3Z";
	elif CohoGpi eq [ 2, 2, 2 ] then
		return "Z/2Z + Z/2Z + Z/2Z";
	elif CohoGpi eq [ 2, 2, 2, 2 ] then
                return "Z/2Z + Z/2Z + Z/2Z + Z/2Z";
        elif CohoGpi eq [ 2, 2, 2, 2, 2, 2 ] then
                return "Z/2Z + Z/2Z + Z/2Z + Z/2Z + Z/2Z + Z/2Z";
        elif CohoGpi eq [ 2, 2, 2, 2, 2, 2, 2, 2 ] then
                return "Z/2Z + Z/2Z + Z/2Z + Z/2Z + Z/2Z + Z/2Z+ Z/2Z + Z/2Z";
        elif CohoGpi eq [ 6, 6 ] then
                return "Z/6Z + Z/6Z";
	elif CohoGpi eq [ ] then
		return "{1}";
	else
		return "Z/6Z";
	end if;
end function;

LineToNumber := function(givenline)
	for i in [1..240] do 
		if givenline eq Line[i] then
			return i;
		end if;
	end for;
	return "ummm";
end function;

/* The following function takes the i-th subgroup Hi of a SubGpLat, and computes 
the orbit of the line j under Hi. If the orbit consists of skew-lines then the Galois 
action of Hi corresponds to that of a non-minimal surfaces. */

GaloisStable := function(Hi,i,j)
	jOrbit := { LineToNumber(Line[j]*(MatrixRing(Integers(),9)!MATR)) : MATR in Hi};
	LocalIntMatrix := MatrixRing(Integers(),#jOrbit)![ IntersectionMatrix[i,j] : 
		i in jOrbit, j in jOrbit ];
	if LocalIntMatrix eq -IdentityMatrix(Integers(),#jOrbit) then
	//NumberOfNonZeroEntries(LocalIntMatrix) ne 1 then
		/*
		print "We are in trouble!";
		print "Subgroup # ", i, "has a Galois Stable set of lines!";
		print "The orbit of line ", j, "is bad: ", jOrbit;
		genes := ComputeGpGens(i, SubGpLat);
		print "Order of group:", #Hi, "Generators:", genes;
		*/
	return true, jOrbit;
	end if;
	return false, jOrbit;
end function;

// initialize some variables and files

counter := 0;
counterForGoodNormalSubgroups := 0;
galstablegroups := 0;
notenough := 0;
PossibleH1s := [];
FL := Open("DataFile", "w");
FL2 := Open("CyclicSearchSpace", "w");

for i in [2..#SubGpLat] do
	Hi := SubGpLat[i]`subgroup;
	if #quotmap(Hi) eq QuotientSize then;
                M2i := GModule(Hi);
                CohoGpi := Moduli(CohomologyGroup(CohomologyModule(Hi,M2i),1));
		PossibleH1s := Append(PossibleH1s,CohoGpi);
			DoWeProceed := true;
			for j in [1..240] do
				trouble, jOrbit := GaloisStable(Hi,i,j);
				if trouble eq true then
					galstablegroups := galstablegroups + 1;
					DoWeProceed := false;
					break;
				end if;
			end for;
			if DoWeProceed then		
				TheGroup := writeGroup(CohoGpi);
				Put(FL, "Subgroup #" cat Sprint(i));
				Put(FL, " Order = " cat Sprint(#Hi));
				Put(FL, " Br X/Br Q = " cat TheGroup);
				Puts(FL, " # Conj. Classes = " cat Sprint(SubGpLat[i]`length));
				genes := ComputeGpGens(i, SubGpLat);
				Puts(FL, "Generators: " cat Sprint(genes));
				counter := counter + 1;
			end if;	

	/* Let M = Pic Xbar, and G = Galois Group of the splitting 
	field of the exceptional divisors. In this case
	G = Hi. The procedure below looks for normal H subgroups
	of G with nontrivial cyclic quotient. Say this quotient is
	generated by rho and has order n. We compute H^1(G/H,M^H) 
	via Tate cohomology as
		kernel(1 + rho + ... + rho^{n-1})/im(1 - rho).
	The result is stored in the file CyclicSearchSpace.
	The parameter CohoModuli below allows us to restrict to
	the case where H^1(G, Pic Xbar) = CohoModuli. */
		
	// CohoModuli := [ 2, 2 ];

		if DoWeProceed and CohoGpi ne [] then
			counterForGoodNormalSubgroups := counterForGoodNormalSubgroups + 1;
			print "Examine Group ", i;
			NewSubGpLat := Subgroups(Hi);
			for j in [1..#NewSubGpLat] do
				Hij := NewSubGpLat[j]`subgroup;
				if IsNormal(Hi,Hij) and IsCyclic(Hi/Hij) and Order(Hi/Hij) ne 1 then
					IntersectionSoFar := RSpace(Integers(),9);
					gens := Setseq(Generators(Hij));
					for k in [1..#Generators(Hij)] do
						Am := Eigenspace(M!gens[k],1);
						IntersectionSoFar := (IntersectionSoFar meet Am);
					end for;
					// Note that (Pic Xbar)^Hij = IntersectionSoFar;
					/* since the canonical class is always in (Pic Xbar)^Hij
					we check to make sure there is something else in this
					group. Otherwise H^1(Hi/Hij,(Pic Xbar)^Hij) = 0 and this
					case is of no interest. */
					if Dimension(IntersectionSoFar) gt 1 then
						Puts(FL2, "Surface X with Br X/Br k = " cat Sprint(writeGroup(CohoGpi)));
						Puts(FL2, "Order of G = " cat Sprint(Order(Hi)));
                                        	Puts(FL2, "This group has generators");
                                        	Puts(FL2, Sprint(ComputeGpGens(i, SubGpLat)));
                                        	Puts(FL2, "Order of H = " cat Sprint(Order(Hij)));
                                        	Puts(FL2, "This subgroup has generators");
                                        	Puts(FL2, Sprint(ComputeGpGens(j, NewSubGpLat)));
						endo := EndomorphismAlgebra(IntersectionSoFar);
						I := IdentityMatrix(Integers(),9);
						Id := endo!I;
						Qu, quomap := quo<Hi | Generators(Hij)>;
						for h in Hi do
							if Order(h) eq Order(Hi/Hij) then
								rho := h;
								break;
							end if;
						end for;
						rho := endo!(rho);
						ker := Kernel(&+[rho^i : i in [1..Order(rho)]]) meet IntersectionSoFar;
						NewGen := Setseq(Generators(IntersectionSoFar));
						im := sub< IntersectionSoFar | [ i*(Id - rho) : i in NewGen ]>;
						Put(FL2, "(Pic Xbar)^H =");
						Puts(FL2, Sprint(IntersectionSoFar));
						Puts(FL2, "Kernel(Norm of generator) = " cat Sprint(ker));
						Puts(FL2, "Image(Id - generator)" cat Sprint(im));
						Puts(FL2, "----------------------------------------------------------------------");
					end if;
				end if;
			end for;
		end if;
	else
		notenough := notenough + 1;
        end if;
end for;
print "The possible groups Br X/ Br k are:", Seqset(PossibleH1s);
delete FL;
delete FL2;
print "Number of conjugacy classes of subgroups with nontrivial H^1 and no Galois Stable orbit of points:", counter;
print "Number of potential groups for counterexample to Weak Approximation:", counterForGoodNormalSubgroups;
print "# of groups with galois stable sets of skew lines:", galstablegroups;
print "# of subgroups G0 for which G0 --> G/<G3,G4> is not surjective:", notenough;