// This is the MAGMA file accompanying the paper 'On the arithmetic of a family of degree two K3 surfaces'.
// In here one can find the codes that are complementary to the proofs 
// in the thesis.
// In the code, the surface X_eta will be here denoted with S

// In order to perform any part of the code, 
// is often require to run also some previous parts.
// To avoid errors we therefore encourage to run the code from the very beginning to 
// the desired point.

/***************************************************************/
/*                                                             */
/*      Proof of Remark 2.15                                   */
/*                                                             */
/***************************************************************/



QQ:=Rationals();

Pol<b,c,a0,a1,a2,a3,D>:=PolynomialRing(QQ,7);

Polt<t>:=PolynomialRing(Pol);

P2<x,y,z>:=PolynomialRing(Polt,3);
f:=x^6+y^6+z^6+D*x^2*y^2*z^2;

fg:=Evaluate(f,[-b*t-c,t,1]);
sq:=(a0*t^3+a1*t^2+a2*t+a3)^2;
df:=fg-sq;
Coeff:=Coefficients(df);
I:=Ideal([Pol!g : g in Coeff]);
1 in I;

GB:=GroebnerBasis(I);

Polw<w>:=PolynomialRing(QQ);

n:=#GB;

h:=GB[n];

h:=Polw!Evaluate(h,[0,0,0,0,0,0,w]);
h;
Factorisation(h);

K1<z3>:=CyclotomicField(3);
Roots(h,K1);

/***************************************************************/
/*                                                             */
/*      End of proof of Remark 2.15                            */
/*                                                             */
/***************************************************************/


// We define the field of definition of the geometric Picard lattice of X_eta, denoted by L,
// together with the weighted projective space PP_L and
// the surface itself

preK<zeta12>:=CyclotomicField(12);
zeta3:=zeta12^4;
zeta4:=zeta12^3;
zeta6:=zeta12^2;
K<t>:=FunctionField(preK);
R<bb0,bb1,bb2>:=PolynomialRing(K,3);
K1<b0,b1,b2>:=quo<R|Ideal([bb0^2 - t-3, bb1^2 - t - 3*zeta3, bb2^2-(t+3*zeta3^2)])>;
S<cc0>:=PolynomialRing(K1);
L<c0>:=quo<S|cc0^3+t*cc0^2+4>;
delta:=4*zeta4*b0*b1*b2;
eps:=delta/(c0*(3*c0+2*t));
c1:=(-t-c0+eps)/2;
c2:=-t-c0-c1;
assert c2 eq (-t-c0-eps)/2;

PT<x,y,z,w>:=ProjectiveSpace(L,[1,1,1,3]);
T:=CoordinateRing(PT);

fS:=x^6+y^6+z^6+t*x^2*y^2*z^2-w^2;
S:=Scheme(PT,fS);

// In order to compute the Gram matrix, we will need 
// to perform some computations over a finite field FF_p^n.
// Let L_t0 be the number field given by specialising  L to t=t0=7.
// One can see that p=79 is a prime of good reduction for X_7.
// Let pp be a prime of L,t0 lying above p=79,
// then the residue field of pp is FF_{79^2}.


p:=79;
F<alpha>:=GF(p^2);
tp:=F!7;

Pol<v>:=PolynomialRing(F);

zeta12p:=Roots(v^12-1)[12][1];
zeta6p:=zeta12p^2;
zeta4p:=zeta12p^3;
zeta3p:=zeta12p^4;
h:=v^3+tp*v^2+4;
cp:=Roots(h);
c0p:=cp[1][1];
c1p:=cp[2][1];
c2p:=cp[3][1];
_,b0p:=IsSquare(3+tp);
_,b1p:=IsSquare(3*zeta3p+tp);
_,b2p:=IsSquare(3*zeta3p^2+tp);

PP<X,Y,Z,W>:=ProjectiveSpace(F,[1,1,1,3]);
PPW:=CoordinateRing(PP);

// We define the reduction maps mod p
// namely the reduction map of fields and polynomial rings.
// and also the reduction of S mod p.
// Using these maps we define the reduction mod p=79
// of X_7.

prePsi:=hom<preK->F | zeta12p>;
Psi:=hom<K->F| prePsi, tp>;
Psi1:=hom<K1->F | Psi, [b0p,b1p,b2p]>;
Psi2:=hom<L->F | Psi1, c0p>;
PsiP:=hom<T->PPW | Psi2, [X,Y,Z,W]>;


assert Psi2(c1) eq c1p;
assert Psi2(c2) eq c2p;

fSp:=PsiP(fS);
Sp:=Scheme(PP,fSp);


// We define the subgroup H of the automorphism of X_eta

Var:={x,y,z};
Perm:=Permutations(Var);


idT:=hom<T->T | x,y,z,w>;
pxy:=hom<T->T | y,x,z,w>;
pxyz:=hom<T->T | y,z,x,w>;
p105:=hom<T->T | zeta6*x,y,zeta6^(-1)*z,w >;
p015:=hom<T->T | x,zeta6*y,zeta6^(-1)*z,w >;
p003:=hom<T->T | x,y,-z,w >;
genH:=[pxy,pxyz,p105,p015,p003];

S3:=[hom<T -> T | perm[1],perm[2],perm[3],w> : perm in Perm];
C2a:=[idT,p003];
C6a:=[hom<T->T | zeta6^i*x, y, zeta6^(-i)*z,w >: i in [0..5]];
C6b:=[hom<T->T | x, zeta6^i*y, zeta6^(-i)*z,w >: i in [0..5]];
C2b:=[p015,idT];
HH:=car<S3,C6a,C2b,C2a>;
H:=[s[1]*s[2]*s[3]*s[4] : s in HH];
assert #H eq 36*4;

/*
// we check that H keeps S fixed

for s in H do
    assert s(fS) eq fS;
end for;
*/

// Now we define the divisors B_1,...,B_5 inside Omega
// and the set H*Omega

a4:=(9*c0+6*t)*delta/4/(t^3+27);
b4:=-c0^2-t*c0;
c4:=(18-3*t^2*c0-3*t*c0^2)/8/(t^3+27)*delta;

a5:=zeta12*(-zeta6+2)*(b0*b1+b0*b2+b1*b2+t)/9;
c5:=zeta12*(zeta6-2)/3;
r5:=zeta12*(zeta6-2)*(2*b0*b1*b2+(2*t-3)*b0+(2*t-3*zeta3)*b1+(2*t+3*zeta6)*b2)/9;
v5:=-b0-b1-b2;

Omega:=
[

[x^2 + y^2 + zeta3*z^2, w-b1*x*y*z],

[x^2 + zeta3*y^2 + zeta3^2*z^2, w-b0*x*y*z],

[2*x*y-c1*z^2, x^3-y^3-w],

[c0*x^2 - 2*(9 * c0^2 + 3 * t * c0 - 2 * t^2) / delta *x*y + 2*y^2 - z^2,
(x^3+a4*x^2*y+ b4*x*y^2 + c4*y^3)*(-c0^2*c1/2+1)-w],

[a5*x^2 +c5*(y^2+z^2)+y*z,
r5*x^3+ v5*x*y*z-w] 

];

HOmega:={};

for eqn in Omega do
    for t in H do
	eqn11:=t(eqn[1]);
	eqn12:=t(eqn[2]);
	if Psi2(LeadingCoefficient(eqn11)) ne F!0 and Psi2(LeadingCoefficient(eqn12)) ne F!0 then
	    Include(~HOmega, [eqn11/LeadingCoefficient(eqn11), eqn12/LeadingCoefficient(eqn12)]);
	else
	    Include(~HOmega, [eqn11, eqn12]);
	end if;
    end for;
end for;

#HOmega;






/***************************************************************/
/*                                                             */
/*      Computations 3.3                                       */
/*                                                             */
/***************************************************************/


function GramMatrixD(List)

m:=#List;
M:=ScalarMatrix(Integers(),m,-2);

for i in [1..m] do
    for j in [1..i] do
	Int:=List[i] meet List[j];
        if Dimension(Int) eq 1 then
	    M[i,j]:=2*ArithmeticGenus(Int)-2;
	else
	    M[i,j]:=Degree(Int);
	end if;
	M[j,i]:=M[i,j];
    end for;
end for;

return M;

end function;




DivsH:=HOmega;
DivsHd:={[eqn[1]/LeadingCoefficient(eqn[1]),eqn[2]/LeadingCoefficient(eqn[2])] : eqn in DivsH};
assert #DivsH eq #DivsHd;


DivsHp:=[[PsiP(eqn[1]),PsiP(eqn[2])]:eqn in HOmega];
DivsHpp:={[eqn[1]/LeadingCoefficient(eqn[1]),eqn[2]/LeadingCoefficient(eqn[2])] : eqn in DivsHp};
assert #DivsHp eq #DivsHpp;

ListDivsHp:=[Scheme(PP,eqn) : eqn in DivsHp];


AG:=GramMatrixD(ListDivsHp);
LG:=Lattice(AG);
BG:=Basis(LG);
rk:=#BG;
SolG:=Solution(AG,[Vector(b) : b in BG]);
MG:=Matrix(SolG);
GramG:=MG*AG*Transpose(MG);

assert Rank(GramG) eq 19;
assert Determinant(GramG) eq 2^5*3^3;

Lambda:=LatticeWithGram(GramG: CheckPositive:= false);
ALambda,DLambda,phi:=DualQuotient(Lambda);
pSignature(Lambda,-1) eq -17;
ALambda;

/*   auxiliary function   */

function Rep2(x)
    return x-2*Round(x/2);
end function;

/*   end auxiliary function    */



mE8:=-Matrix(GramMatrix(Lattice("E",8)));
mE8:=ChangeRing(mE8,Integers());
U:=Matrix([[0,1],[1,0]]);
mA5:=-Matrix(GramMatrix(Lattice("A",5)));
mA2:=-Matrix([[2,-1],[-1,2]]); 
m4A2:=-4*Matrix([[2,-1],[-1,2]]);


Dec1:=<  <mE8,1>, <U,1>, <mA2,1>, <mA5,1> , <m4A2,1>>;
L1:=Dec1[1,1];
for i in [1..#Dec1] do
    if i eq 1 and Dec1[i,2] ne 1 then
	for j in [2..Dec1[1,2]] do
	    L1:=DiagonalJoin(L1,Dec1[i,1]);
        end for;
    end if;
    if i gt 1 then
	for j in [1..Dec1[i,2]] do
	    L1:=DiagonalJoin(L1,Dec1[i,1]);
	end for;
    end if;
end for;

L1:=LatticeWithGram(L1: CheckPositive:= false);
AL1,DL1,phi1:=DualQuotient(L1);

q:={* Rep2(Norm((phi^-1)(g))) : g in ALambda *};
q1:={* Rep2(Norm((phi1^-1)(g1))) : g1 in AL1 *};

assert pSignature(L1,-1) eq pSignature(Lambda,-1);
assert Rank(L1) eq Rank(Lambda);
assert Factorisation(Determinant(L1)) eq Factorisation(Determinant(Lambda));
assert #Generators(AL1) eq #Generators(ALambda);
assert q eq q1;
assert IsIsomorphic(ALambda,AL1);



/***************************************************************/
/*                                                             */
/*      End of Computations 3.3                                */
/*                                                             */
/***************************************************************/

// We define the Galois group Gal(L/K)

Gal:=[];
EPhi12:=[1,5,7,11];
Ipm:=[1,-1];
GalK1K:=CartesianPower(Ipm,3);
cc:=[L!c0,L!c1,L!c2];
GalLE:=CartesianPower(Ipm,2);


for i in EPhi12 do
    presigma := hom<preK->L| zeta12^i>;
    sigma := hom<K->L | presigma, t>;
    for s in GalK1K do
        if sigma(zeta3) eq L!zeta3 then
	    sigma1:=hom<K1->L | sigma, [s[1]*b0,s[2]*b1,s[3]*b2]>;
	else
	    sigma1:=hom<K1->L | sigma, [s[1]*b0,s[2]*b2,s[3]*b1]>;
	end if;
	for j in [1..#cc] do
	    c:=cc[j];
	    sigma2:=hom<L->L | sigma1, c>;
	    Append(~Gal,sigma2);
	end for;
    end for;
end for;

// we name a set of generators genGal:=[tau1,...,tau5] of Gal(L/K) 


gen:=[zeta12,b0,b1,b2,c0];
ImTau:=[
[zeta12^7, b0, b1, b2, c0],
[zeta12, b0, b1, b2, c1],
[zeta12^7, -b0, b1, b2, c0],
[zeta12^11, b0, -b2, b1, c0],
[zeta12^7, b0, b1, -b2, c0]
];

genGal:=[];

for i in [1..#ImTau] do
    for s in Gal do
	if [s(gg) : gg in gen] eq ImTau[i] then
	    Append(~genGal,s);
	end if;
    end for;
end for;

GalT:=[hom<T->T | s,[x,y,z,w]> : s in Gal];
genGalT:=[hom<T->T | s,[x,y,z,w]> : s in genGal]; 

// We let GalT act on HOmega
// getting a list of divisor that is invariant 
// under the action of H and Gal.

HGalOmega:={};
for eqn in HOmega do
    for s in GalT do
	eqn11:=s(eqn[1]);
	eqn12:=s(eqn[2]);
	if Psi2(LeadingCoefficient(eqn11)) ne F!0 and Psi2(LeadingCoefficient(eqn12)) ne F!0 then
	    Include(~HGalOmega, [eqn11/LeadingCoefficient(eqn11), eqn12/LeadingCoefficient(eqn12)]);
	else
	    Include(~HGalOmega, [eqn11, eqn12]);
	end if;    
    end for;
end for;


GOmega:=[eqn : eqn in HGalOmega];

#GOmega;

// We skim the list GOmega deleting the repetitions,
// getting two different lists of equations of divisors:
// SetDivs5 and SetDivsp.
// The elements of the list at the same position are the same up to a scalar
// SetDivsp will be used to reduce the equation mod p
// SetDivsG will be used to compute the matrices 
// corresponding to the action of G=<H,Gal> on Pic(S)

SetDivsNP:=[eqn : eqn in GOmega | {LeadingCoefficient(eqn[1]),LeadingCoefficient(eqn[2])} eq {1}];

SetDivsP:=[eqn : eqn in GOmega | {LeadingCoefficient(eqn[1]),LeadingCoefficient(eqn[2])} ne {1}];

#SetDivsP;
#SetDivsNP;
assert #SetDivsP+#SetDivsNP eq #GOmega;

SetDivsP1:=[[eqn[1]/LeadingCoefficient(eqn[1]),eqn[2]/LeadingCoefficient(eqn[2])]: eqn in SetDivsP];
SetDivsP2:=Set(SetDivsP1);
IndProb:=[Index(SetDivsP1,eqn) : eqn in SetDivsP2 ];
SetDivsP3:=[SetDivsP[i] : i in IndProb];
SetDivsP4:=[[eqn[1]/LeadingCoefficient(eqn[1]),eqn[2]/LeadingCoefficient(eqn[2])]: eqn in SetDivsP3];



SetDivs5:= SetDivsNP cat SetDivsP4;
SetDivsG:=SetDivs5;
SetDivsp:=SetDivsNP cat SetDivsP3;
assert #SetDivsp eq #SetDivs5;
assert #SetDivsp eq #GOmega;

Divsp:=[[PsiP(eqn[1]),PsiP(eqn[2])]:eqn in SetDivsp];
ListDivsp:=[Scheme(PP,eqn) : eqn in Divsp];



/***************************************************************/
/*                                                             */
/*      Proof of Lemma 3.7                                     */
/*                                                             */
/***************************************************************/

// In order to prove that G preserves Lambda
// We shot that the lattice generated by the elements of GOmega
// is exactly Lambda.


A0:=GramMatrixD(ListDivsp);
L0:=Lattice(A0);
B0:=Basis(L0);
rk:=#B0;
Sol0:=Solution(A0,[Vector(b) : b in B0]);
M0:=Matrix(Sol0);
Gram0:=M0*A0*Transpose(M0);

assert Rank(GramG) eq Rank(Gram0);
assert Determinant(GramG) eq Determinant(Gram0);


// By construction of the list GOmega,
// the group G acts as permutation on it.
// We want to compute the 414x414 permutation matrices that correspond to the 
// elements of G.
// Notice that they will all be 414x414 permutation matrices.
// Since these divisors generate a rank 19 lattice
// we can write these matrices as 19x19 matrices


nD:=#SetDivsG;
GalGraphs:=[];

for s in genGalT do
    Graph:=[];
    for i in [1..nD] do
	eqn:=SetDivsG[i];
	eqn1:=[s(eqn[1]),s(eqn[2])];
	j:=Index(SetDivsG,[eqn1[1]/LeadingCoefficient(eqn1[1]), 
eqn1[2]/LeadingCoefficient(eqn1[2])]);
	Append(~Graph,j);
    end for;
    Append(~GalGraphs,Graph);
end for;
 
MgenGal:=[];         // 414x414 matrices corresponding to the generators 
                     // of the Galois group

for Gr in GalGraphs do
    Ms:=[[0: i in [1..nD]] : j in [1..nD]];
    for i in [1..nD] do
	Ms[i][Gr[i]]:=1;
    end for;
    Append(~MgenGal,Matrix(Ms));
end for;


HGraphs:=[];

for t in genH do
    Graph:=[];
    for i in [1..nD] do
	eqn:=SetDivsG[i];
	eqn1:=[t(eqn[1]),t(eqn[2])];
	j:=Index(SetDivsG,[eqn1[1]/LeadingCoefficient(eqn1[1]), 
eqn1[2]/LeadingCoefficient(eqn1[2])]);
	Append(~Graph,j);
    end for;
    Append(~HGraphs,Graph);
end for;
 

MgenH:=[];                   // 414x414 matrices corresponding to the generators
			      // of H

for Gr in HGraphs do
    Ms:=[[0: i in [1..nD]] : j in [1..nD]];
    for i in [1..nD] do
	Ms[i][Gr[i]]:=1;
    end for;
    Append(~MgenH,Matrix(Ms));
end for;

MgenG:=MgenH cat MgenGal;   // 414x414 matrices corresponding to the generators
			        // of G


// The matrix M0 defines a map from ZZ^19 to ZZ^414 
// We need a section of this map, represented by the matrix N0
// With this matrix we compute the 19x19 matrices 
// representing the generators of H and Gal

BB0:=Matrix(B0);
Id19:=ScalarMatrix(19,1);
N0:=Matrix(Solution(BB0,[A0[i] : i in [1..#ListDivsp]]));

assert M0*N0 eq Id19;
assert N0*Gram0*Transpose(N0) eq A0;

M19genG:=[M0*M*N0 : M in MgenG];       // 19x19 matrices corresponding
           				       // to the generators of
					       // G, that is,
					       // the union of the generators
					       // of H and Gal
 
M19genH:=[M0*M*N0 : M in MgenH];             // 19x19 matrices corresponding
           				       // to the generators of
					       // H. 

M19genGal:=[M0*M*N0 : M in MgenGal];         // 19x19 matrices corresponding
           				       // to the generators of
					       // Gal.

G19:=MatrixGroup<19,Integers()|M19genG>;   // The group of 19x19 matrices
					       // corresponding to <H,Gal>.

H19:=MatrixGroup<19,Integers()|M19genH>;     // The group of 19x19 matrices
					       // corresponding to H.

Gal19:=MatrixGroup<19,Integers()|M19genGal>; // The group of 19x19 matrices
					       // corresponding to Gal.

// In order to prove that there is no element of G acting as the identity on Lambda
// it is enough to check that the cardinality of G19 equals the product of the cardinalities
// of H and Gal, or equivalently, that 
// the cardinality of Gal19 equals the cardinality of Gal 

assert #G19 eq #H*#Gal;
assert #Gal19 eq #Gal;


/***************************************************************/
/*                                                             */
/*      End of Proof of Lemma 3.7                              */
/*                                                             */
/***************************************************************/

// We give names to the matrices corresponding to the generators of H and Gal


tM1:=M19genGal[1];
tM2:=M19genGal[2];
tM3:=M19genGal[3];
tM4:=M19genGal[4];
tM5:=M19genGal[5];

pxyM:=M19genH[1];
pxyzM:=M19genH[2];
p105M:=M19genH[3];
p015M:=M19genH[4];
p003M:=M19genH[5];



/***************************************************************/
/*                                                             */
/*          Proof of Lemma 3.9                                 */
/*                                                             */
/***************************************************************/


/* auxuliary function  */
// This function compute the self intersection of 
// the element v of a lattice L having M as intersection form  

function QuadInt(v,w,M)

v1:=ChangeRing(Matrix(v),Integers());
w1:=ChangeRing(Matrix(w),Integers());

return (v1*M*Transpose(w1))[1,1];

end function;

/* end auxiliary function */




// First we prove the injectivity of iota_p for p=2.
// We compute M_p, denoted by 

F2:=GF(2);
V2:=VectorSpace(F2,19);
Gram2:=ChangeRing(Gram0,F2);
k2:=sub<V2 | Kernel(Gram2)>;
k2;

GenH2:=[ChangeRing(M,F2) : M in M19genH];
H2:=MatrixGroup<19,F2|GenH2>;
H2t:=[M : M in Set(H2)];

M2:=[v : v in k2 | QuadInt(v,v,Gram0) mod 8 eq 0];

VM2:=[v : v in M2 | Dimension(sub<V2 | Orbit(H2,v)>) le 2];
W:=sub<V2 | VM2>;
assert Set(W) eq Set(VM2);

// we compute a basis of W
v1:=Basis(W)[1];
v2:=Basis(W)[2];

// We see which representative of v1 has been used by MAGMA

elv1:=Eltseq(v1);
u1:=&+[elv1[i]*B0[i] : i in [1..#B0]];                // We write w1 as combination
							 // of the elements of the basis
							 // of L0.

uu,N:=Solution(A0,u1);                                // we write it as combination
							 // of -2 curves.
							 // It does not need to be minimal

// We pick a particular representative w1 of v1, 
// such that w1^2 = -8
// and we check that ut is indeed a representative of v1

ww:=[ 0  : i in [1..414]];
ww[37]:=1;
ww[70]:=-1;

ww:=Parent(uu)!(ww);
w1:=ww*N0;

assert QuadInt(w1,w1,Gram0) eq -8;
assert V2!w1 eq v1;



eqn:=Omega[4];
B4:=[eqn[1]/LeadingCoefficient(eqn[1]),eqn[2]/LeadingCoefficient(eqn[2])];



p030:=hom<T->T | x,-y,z,w >;
eqn:=[p030(B4[1]),p030(B4[2])];
p030B4:=[eqn[1]/LeadingCoefficient(eqn[1]),eqn[2]/LeadingCoefficient(eqn[2])];


// We set E0:=\tau_2^2\sigma_{(xy)}p_{0,3,0}(B4)
// and E2:= \tau_2^2\sigma_{(xy)}(B4)

tau2:=genGalT[2];

eqn:=[tau2(tau2(pxy(p030B4[1]))),tau2(tau2(pxy(p030B4[2])))];
E0:=[eqn[1]/LeadingCoefficient(eqn[1]),eqn[2]/LeadingCoefficient(eqn[2])];
e0:=Index(SetDivsG,E0);

eqn:=[tau2(tau2(pxy(B4[1]))),tau2(tau2(pxy(B4[2])))];
E2:=[eqn[1]/LeadingCoefficient(eqn[1]),eqn[2]/LeadingCoefficient(eqn[2])];
e2:=Index(SetDivsG,E2);



ww:=[ 0  : i in [1..414]];
ww[e0]:=1;
ww[e2]:=-1;

ww:=Parent(uu)!(ww);
w2:=ww*N0;

assert QuadInt(w2,w2,Gram0) eq -8;
assert V2!w2 eq v2;


// Then we prove the injectivity for p=3

F3:=GF(3);
V3:=VectorSpace(F3,19);
Gram3:=ChangeRing(Gram0,F3);
k3:=sub<V3 | Kernel(Gram3)>;
k3;

GenG3:=[ChangeRing(M,F3) : M in M19genG];
G3:=MatrixGroup<19,F3|GenG3>;

M3:=[v : v in k3 | QuadInt(v,v,Gram0) mod 18 eq 0];

VM3:=[v : v in M3 | Dimension(sub<V3 | Orbit(G3,v)>) le 1];

assert VM3 eq [V3![0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]];


/***************************************************************/
/*                                                             */
/*          End of proof of Lemma 3.9                          */
/*                                                             */
/***************************************************************/

/***************************************************************/
/*                                                             */
/*      Proof of Computation 4.1                               */
/*                                                             */
/***************************************************************/



// Using the matrix representation of the action of G on L0
// We can compute the H^0(Gal,Pic), 

gmpic:=GModule(Gal19);
cmpic:=CohomologyModule(Gal19,gmpic);
cgpic:=CohomologyGroup(cmpic,0);
cgpic;


// Using the matrix representation of the action of GG on L0
// We can compute the H^1(Gal,Pic), 
// and H^1(H,Pic) for any H<Gal subgroup, not necessarily normal, subgroup of G

gmpic:=GModule(Gal19);
cmpic:=CohomologyModule(Gal19,gmpic);
cgpic:=CohomologyGroup(cmpic,1);
cgpic;

subs:=[H`subgroup : H in Subgroups(Gal19) | #H`subgroup ne 1];
nzsubs:=[];
zsubs:=[];
indxs:={};

for H in subs do
  gmpicH:=GModule(H);
  cmpicH:=CohomologyModule(H,gmpicH);
  cgpicH:=CohomologyGroup(cmpicH,1);
  cgpicH;
  Order(H);
  Degree(cgpicH);
  if Degree(cgpicH) ne 0 then
    Append(~nzsubs,H);
    Include(~indxs,Degree(cgpicH));
  else
    Append(~zsubs,H);
  end if;
  "----";
end for;

// We check how many normal subgroups of Gal
// have non-zero or zero H^1.

zNsubs:=[H : H in zsubs | IsNormal(Gal19,H)];
nzNsubs:=[H : H in nzsubs | IsNormal(Gal19,H)];


// Using the matrix representation of the action of GG on L0
// We can compute the H^0(Gal,Pic), 

gmpic:=GModule(Gal19);
cmpic:=CohomologyModule(Gal19,gmpic);
cgpic:=CohomologyGroup(cmpic,0);
cgpic;


/***************************************************************/
/*                                                             */
/*      End of proof of Computation 4.1                        */
/*                                                             */
/***************************************************************/