//  This is the code accompanying the paper 
// 'Bhabha scattering and a special family of K3 surfaces' 
// by Dino Festi and Duco van Straten. 
// It follows the order of the sections in the paper

////////////////////
//   Section 1    //
////////////////////


////////////////////
//   Section 2    //
////////////////////


K<s>:=FunctionField(Rationals());
R<u,v,t>:=FieldOfFractions(PolynomialRing(K,3));
X:=1+s+1/(u^2-1)+1/(v^2-1)+1/(t^2-1);
assert R!Numerator(X) eq u^2*v^2*t^2-u^2-v^2-t^2+2+s*(u^2-1)*(v^2-1)*(t^2-1);
X:=R!Numerator(X);
assert X eq t^2*(u^2*v^2-1+s*(u^2*v^2-u^2-v^2+1))-(u^2+v^2-2+s*(u^2*v^2-u^2-v^2+1));

P2<x,y,z>:=ProjectivePlane(K);
T:=z^4-x^2*y^2; S:=(x^2-z^2)*(y^2-z^2);
f1h:=T+(s-1)*S;
f0h:=T+s*S;
Q1:=Curve(P2,f1h);
Q0:=Curve(P2,f0h);
SingularPoints(Q1);
SingularPoints(Q0);
RationalPoints(Q1 meet Q0);

// Here we perform the Cremona transformation
Ls:=[(z-y)*(z-x), (z-y)*z, z*(z-x)];
Ct:=map<P2 -> P2 | Ls>;

RationalPoints(BaseScheme(Ct));

g0:=Evaluate(f0h,[x+z, y+z, z]);
g1:=Evaluate(f1h,[x+z, y+z, z]);

B0:=Scheme(P2,g0);
gg0:=DefiningEquation(QuadraticTransformation(B0));
h0:=Evaluate(gg0,[x-z,y-z,z]);

B1:=Scheme(P2,g1);
gg1:=DefiningEquation(QuadraticTransformation(B1));
h1:=Evaluate(gg1,[x-z,y-z,z]);


T1:=(x^2+y^2)*z-2*x*y*(x+y);
S1:=(2*x-z)*(2*y-z)*z;
g1:=T1+s*S1;
g0:=T1+(s+1)*S1;
g0 eq h0;
g1 eq h1;




////////////////////
//   Section 3    //
////////////////////


// Proof of Lemma 3.1

K<s>:=FunctionField(Rationals());
L<a>:=ext<K|Polynomial([-s^2+s,0,1])>;
PP<x,y,z,w>:=ProjectiveSpace(L,[1,1,1,3]);
T:=(x^2+y^2)*z-2*x*y*(x+y);
S:=(2*x-z)*(2*y-z)*z;
f1:=T+s*S;
f2:=T+(s+1)*S;
f:=f1*f2;

ll:=[0,2*x, 2*y, x+y, x/(s+a), x/(s-a), y/(s+a), y/(s-a)];
/*
for l in ll do
  l;
  Evaluate(f1*f2,[x,y,l,0]);
  Factorisation($1);
end for;
*/
Ls:=[
[z, w-2*x*y*(x+y)],
[z-2*x, w-2*x^2*(x-y)],
[z-2*y, w-2*y^2*(y-x)],
[z-x-y, w-a*(x-y)^2*(x+y)],

[z-(1-1/s*a)*x, w-x*y*(x-(a+s)*y)/s],
[z-(1+1/s*a)*x, w-x*y*(x+(a-s)*y)/s],
[z-(1-1/s*a)*y, w-x*y*(y-(a+s)*x)/s],
[z-(1+1/s*a)*y, w-x*y*(y+(a-s)*x)/s]

];

for l in Ls do
  l;
  assert Evaluate(f,[x,y,-l[1]+z,w]) eq (l[2]-w)^2;
end for;





////////////////////
//   Section 4    //
////////////////////



// Proof of Proposition 4.3

ListDivs:=[Scheme(PP,eqn) : eqn in Ls];
SP:=[[1,0,0,0],[0,1,0,0],[1,1,2,0],[1,-1,0,0]];
SP:=[PP!P : P in SP];

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;
	else
	    M[i,j]:=Degree(Int);
	    if SP[1] in Int or SP[2] in Int or SP[3] in Int or SP[4] in Int then
	        M[i,j]:=0;
	    end if;

	end if;
	M[j,i]:=M[i,j];
    end for;
end for;

return M;

end function;

GramMatrixD(ListDivs);


////////////////////
//   Section 5    //
////////////////////

// Proof of proposition 5.1


for t in CartesianPower([0,1],7) do
L3E31:=t[1];
L3E3m1:=1-L3E31;
L5E22:=t[2];
L5E2m2:=1-L5E22;
L6E22:=t[3];
L6E2m2:=1-L6E22;
L7E12:=t[4];
L7E1m2:=1-L7E12;
L8E12:=t[5];
L8E1m2:=1-L8E12;
L2E21:=t[6];
L2E2m1:=1-L2E21;
L3E11:=t[7];
L3E1m1:=1-L3E11;


A:=Matrix([
[-2, 1, 0, 0, 0,  0, 0, 0, 0, 0,  0, 0, 0,  0,  0,  1, 0, 0, 0,  0, 0, L7E12, L8E12],
[1, -2, 1, 0, 0,  0, 0, 0, 0, 0,  0, 0, 0,  0,  0,  0, 0, L3E11, 0,  0, 0, 0, 0],
[0, 1, -2, 1, 0,  0, 0, 0, 0, 0,  0, 0, 0,  0,  0,  0, 0, 0, 0,  0, 0, 0, 0],
[0, 0, 1, -2, 1,  0, 0, 0, 0, 0,  0, 0, 0,  0,  0,  0, 0, L3E1m1, 0,  0, 0, 0, 0],
[0, 0, 0, 1, -2,  0, 0, 0, 0, 0,  0, 0, 0,  0,  0,  0, 0, 0, 0,  0, 0, L7E1m2, L8E1m2],

[0, 0, 0, 0, 0,  -2, 1, 0, 0, 0,  0, 0, 0,  0,  0,  1, 0, 0, 0,  L5E22, L6E22, 0, 0],
[0, 0, 0, 0, 0,  1, -2, 1, 0, 0,  0, 0, 0,  0,  0,  0, L2E21, 0, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 1, -2, 1, 0,  0, 0, 0,  0,  0,  0, 0, 0, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 0, 1, -2, 1,  0, 0, 0,  0,  0,  0, L2E2m1, 0, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 0, 0, 1, -2,  0, 0, 0,  0,  0,  0, 0, 0, 0,  L5E2m2, L6E2m2, 0, 0],

[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  -2, 1, 0,  0,  0,  0, 1, L3E31, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  1, -2, 1,  0,  0,  0, 0, 0, 1,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  0, 1, -2,  0,  0,  0, 0, L3E3m1, 0,  0, 0, 0, 0],

[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  0, 0, 0,  -2,  0,  1, 0, 0, 1,  0, 0, 0, 0],

[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  0, 0, 0,  0,   2,  1, 1, 1, 1,  1, 1, 1, 1],

[1, 0, 0, 0, 0,  1, 0, 0, 0, 0,  0, 0, 0,  1,  1,  -2, 0, 0, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, L2E21, 0, L2E2m1, 0,  1, 0, 0,  0,  1,  0, -2, 0, 0,  0, 0, 0, 0],
[0, L3E11, 0, L3E1m1, 0,  0, 0, 0, 0, 0,  L3E31, 0, L3E3m1,  0,  1,  0, 0, -2, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  0, 1, 0,  1,  1,  0, 0, 0, -2,  0, 1, 0, 1],

[0, 0, 0, 0, 0,  L5E22, 0, 0, 0, L5E2m2,  0, 0, 0,  0,  1,  0, 0, 0, 0,  -2, 0, 1, 0],
[0, 0, 0, 0, 0,  L6E22, 0, 0, 0, L6E2m2,  0, 0, 0,  0,  1,  0, 0, 0, 1,  0, -2, 0, 1],
[L7E12, 0, 0, 0, L7E1m2,  0, 0, 0, 0, 0,  0, 0, 0,  0,  1,  0, 0, 0, 0,  1, 0, -2, 0],
[L8E12, 0, 0, 0, L8E1m2,  0, 0, 0, 0, 0,  0, 0, 0,  0,  1,  0, 0, 0, 1,  0, 1, 0, -2]


]);

assert IsSymmetric(A);
//Parent(A);
rk:=Rank(A);
if rk le 20 then
t;
L0:=Lattice(A);
B0:=Basis(L0);
Sol0:=Solution(A,[Vector(b) : b in B0]);
M0:=Matrix(Sol0);
Gram0:=M0*A*Transpose(M0);
Rank(Gram0);
Factorisation(Determinant(Gram0));
L:=LatticeWithGram(Gram0: CheckPositive:= false);
AL,_,_:=DualQuotient(L);
AL;
end if;
end for;


\\ Proof of Corollary 5.2

L3E31:=1;
L3E3m1:=1-L3E31;
L5E22:=1;
L5E2m2:=1-L5E22;
L6E22:=1;
L6E2m2:=1-L6E22;
L7E12:=1;
L7E1m2:=1-L7E12;
L8E12:=1;
L8E1m2:=1-L8E12;
L2E21:=1;
L2E2m1:=1-L2E21;
L3E11:=1;
L3E1m1:=1-L3E11;


A:=Matrix([
[-2, 1, 0, 0, 0,  0, 0, 0, 0, 0,  0, 0, 0,  0,  0,  1, 0, 0, 0,  0, 0, L7E12, L8E12],
[1, -2, 1, 0, 0,  0, 0, 0, 0, 0,  0, 0, 0,  0,  0,  0, 0, L3E11, 0,  0, 0, 0, 0],
[0, 1, -2, 1, 0,  0, 0, 0, 0, 0,  0, 0, 0,  0,  0,  0, 0, 0, 0,  0, 0, 0, 0],
[0, 0, 1, -2, 1,  0, 0, 0, 0, 0,  0, 0, 0,  0,  0,  0, 0, L3E1m1, 0,  0, 0, 0, 0],
[0, 0, 0, 1, -2,  0, 0, 0, 0, 0,  0, 0, 0,  0,  0,  0, 0, 0, 0,  0, 0, L7E1m2, L8E1m2],

[0, 0, 0, 0, 0,  -2, 1, 0, 0, 0,  0, 0, 0,  0,  0,  1, 0, 0, 0,  L5E22, L6E22, 0, 0],
[0, 0, 0, 0, 0,  1, -2, 1, 0, 0,  0, 0, 0,  0,  0,  0, L2E21, 0, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 1, -2, 1, 0,  0, 0, 0,  0,  0,  0, 0, 0, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 0, 1, -2, 1,  0, 0, 0,  0,  0,  0, L2E2m1, 0, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 0, 0, 1, -2,  0, 0, 0,  0,  0,  0, 0, 0, 0,  L5E2m2, L6E2m2, 0, 0],

[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  -2, 1, 0,  0,  0,  0, 1, L3E31, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  1, -2, 1,  0,  0,  0, 0, 0, 1,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  0, 1, -2,  0,  0,  0, 0, L3E3m1, 0,  0, 0, 0, 0],

[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  0, 0, 0,  -2,  0,  1, 0, 0, 1,  0, 0, 0, 0],

[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  0, 0, 0,  0,   2,  1, 1, 1, 1,  1, 1, 1, 1],

[1, 0, 0, 0, 0,  1, 0, 0, 0, 0,  0, 0, 0,  1,  1,  -2, 0, 0, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, L2E21, 0, L2E2m1, 0,  1, 0, 0,  0,  1,  0, -2, 0, 0,  0, 0, 0, 0],
[0, L3E11, 0, L3E1m1, 0,  0, 0, 0, 0, 0,  L3E31, 0, L3E3m1,  0,  1,  0, 0, -2, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  0, 1, 0,  1,  1,  0, 0, 0, -2,  0, 1, 0, 1],

[0, 0, 0, 0, 0,  L5E22, 0, 0, 0, L5E2m2,  0, 0, 0,  0,  1,  0, 0, 0, 0,  -2, 0, 1, 0],
[0, 0, 0, 0, 0,  L6E22, 0, 0, 0, L6E2m2,  0, 0, 0,  0,  1,  0, 0, 0, 1,  0, -2, 0, 1],
[L7E12, 0, 0, 0, L7E1m2,  0, 0, 0, 0, 0,  0, 0, 0,  0,  1,  0, 0, 0, 0,  1, 0, -2, 0],
[L8E12, 0, 0, 0, L8E1m2,  0, 0, 0, 0, 0,  0, 0, 0,  0,  1,  0, 0, 0, 1,  0, 1, 0, -2]


]);

assert IsSymmetric(A);
rk:=Rank(A);
L0:=Lattice(A);
B0:=Basis(L0);
Sol0:=Solution(A,[Vector(b) : b in B0]);
M0:=Matrix(Sol0);
Gram0:=M0*A*Transpose(M0);
Rank(Gram0);
Factorisation(Determinant(Gram0));
Lambda:=LatticeWithGram(Gram0: CheckPositive:= false);
ALambda,DLambda,phi:=DualQuotient(Lambda);
pSignature(Lambda,-1) eq -18;
ALambda;





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

mE8:=-Matrix(GramMatrix(Lattice("E",8)));
mE8:=ChangeRing(mE8,Integers());
U:=Matrix([[0,1],[1,0]]);
m12:=Matrix([[-12]]);


Dec1:=<  <mE8,2>, <U,1>, <m12,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);






// Proof of Theorem 5.6


function QuadInt(v,w,M)

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

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

end function;


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


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

M2;



////////////////////
//   Section 6    //
////////////////////


// Subsection 6.1: s=1

s:=1;
PP2<x,y,z>:=ProjectivePlane(Rationals());
T:=(x^2+y^2)*z-2*x*y*(x+y);
S:=(2*x-z)*(2*y-z)*z;
g0:=T+s*S;
g1:=T+(s+1)*S;
g:=g0*g1;
B0:=Curve(PP2,g0);
B1:=Curve(PP2,g1);
Int01:=B0 meet B1;
for P in Points(Int01) do
  P,Multiplicity(P);
end for;

fB0:=Factorisation(g0);
g0p:=fB0[1,1];
g0pp:=fB0[2,1];

B0p:=Curve(PP2,g0p);
B0pp:=Curve(PP2,g0pp);
Int0:=B0p meet B0pp;
for P in Points(Int0) do
  P,Multiplicity(P);
end for;



B:=Curve(PP2,g);
BSPts:=SingularPoints(B);
for P in BSPts do
  Multiplicity(B,P);
  IsOrdinarySingularity(B,P);
end for;



g0px:=Derivative(g0p,x);
g0py:=Derivative(g0p,y);
g0pz:=Derivative(g0p,z);
dg0p:=[g0p,g0px,g0py,g0pz];

g0ppx:=Derivative(g0pp,x);
g0ppy:=Derivative(g0pp,y);
g0ppz:=Derivative(g0pp,z);
dg0pp:=[g0pp,g0ppx,g0ppy,g0ppz];

g1x:=Derivative(g1,x);
g1y:=Derivative(g1,y);
g1z:=Derivative(g1,z);
dg1:=[g1,g1x,g1y,g1z];

for P in BSPts do
  P;
  e:=Eltseq(P);
  Evaluate(dg0p,e);
  Evaluate(dg0pp,e);
  Evaluate(dg1,e);
end for;





for a in CartesianPower([0,1],3) do
 
  L3E31:=a[1];
  L3E3m1:=1-a[1];
  L2L5:=a[2];
  L3L4:=a[3];


A:=Matrix([
[-2, 1, 0, 0, 0,      0, 0, 0, 0, 0,      0, 0, 0,               0,  0,  0,  0,  0,  1, 0, 0, 0, 1],
[1, -2, 1, 0, 0,      0, 0, 0, 0, 0,      0, 0, 0,               0,  0,  0,  0,  0,  0, 0, 1, 0, 0],
[0, 1, -2, 1, 0,      0, 0, 0, 0, 0,      0, 0, 0,               0,  0,  0,  0,  0,  0, 0, 0, 0, 0],
[0, 0, 1, -2, 1,      0, 0, 0, 0, 0,      0, 0, 0,               0,  0,  0,  0,  0,  0, 0, 0, 0, 0],
[0, 0, 0, 1, -2,      0, 0, 0, 0, 0,      0, 0, 0,               0,  0,  0,  0,  0,  0, 0, 0, 0, 0],

[0, 0, 0, 0, 0,      -2, 1, 0, 0, 0,      0, 0, 0,               0,  0,  0,  0,  0,  1, 0, 0, 1, 0],
[0, 0, 0, 0, 0,      1, -2, 1, 0, 0,      0, 0, 0,               0,  0,  0,  0,  0,  0, 1, 0, 0, 0],
[0, 0, 0, 0, 0,      0, 1, -2, 1, 0,      0, 0, 0,               0,  0,  0,  0,  0,  0, 0, 0, 0, 0],
[0, 0, 0, 0, 0,      0, 0, 1, -2, 1,      0, 0, 0,               0,  0,  0,  0,  0,  0, 0, 0, 0, 0],
[0, 0, 0, 0, 0,      0, 0, 0, 1, -2,      0, 0, 0,               0,  0,  0,  0,  0,  0, 0, 0, 0, 0],

[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       -2, 1, 0,              0,  0,  0,  0,  0,  0, 1, L3E31, 0, 0],
[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       1, -2, 1,              0,  0,  0,  0,  0,  0, 0, 0, 0, 0],
[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       0, 1, -2,              0,  0,  0,  0,  0,  0, 0, L3E3m1, 0, 0],


[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       0, 0, 0,               -2,  0,  0,  0,  0,  1, 0, 0, 0, 0],

[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       0, 0, 0,               0,  -2,  0,  0,  0,  0, 0, 0, 1, 0],

[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       0, 0, 0,               0,  0,  -2,  0,  0,  0, 0, 0, 0, 1],

[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       0, 0, 0,               0,  0,  0,  -2,  0,  0, 0, 0, 1, 1],

[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       0, 0, 0,               0,  0,  0,  0,   2,  1, 1, 1, 1, 1],

[1, 0, 0, 0, 0,      1, 0, 0, 0, 0,       0, 0, 0,               1,  0,  0,  0,  1,  -2, 0, 0, 0,  0],
[0, 0, 0, 0, 0,      0, 1, 0, 0, 0,       1, 0, 0,               0,  0,  0,  0,  1,  0, -2, 0, 0,  L2L5],
[0, 1, 0, 0, 0,      0, 0, 0, 0, 0,       L3E31, 0, L3E3m1,      0,  0,  0,  0,  1,  0, 0, -2, L3L4,  0],
[0, 0, 0, 0, 0,      1, 0, 0, 0, 0,       0, 0, 0,               0,  1,  0,  1,  1,  0, 0, L3L4,  -2, 0],
[1, 0, 0, 0, 0,      0, 0, 0, 0, 0,       0, 0, 0,               0,  0,  1,  1,  1,  0,  L2L5, 0, 0,  -2]

]);


assert IsSymmetric(A);
rk:=Rank(A);
if rk eq 20 then
a;
rk;
Determinant(A);
end if;

end for;





L3E31:=1;
L3E3m1:=0;
L3L4:=0;
L2L5:=0;

A:=Matrix([
[-2, 1, 0, 0, 0,      0, 0, 0, 0, 0,      0, 0, 0,               0,  0,  0,  0,  0,  1, 0, 0, 0, 1],
[1, -2, 1, 0, 0,      0, 0, 0, 0, 0,      0, 0, 0,               0,  0,  0,  0,  0,  0, 0, 1, 0, 0],
[0, 1, -2, 1, 0,      0, 0, 0, 0, 0,      0, 0, 0,               0,  0,  0,  0,  0,  0, 0, 0, 0, 0],
[0, 0, 1, -2, 1,      0, 0, 0, 0, 0,      0, 0, 0,               0,  0,  0,  0,  0,  0, 0, 0, 0, 0],
[0, 0, 0, 1, -2,      0, 0, 0, 0, 0,      0, 0, 0,               0,  0,  0,  0,  0,  0, 0, 0, 0, 0],

[0, 0, 0, 0, 0,      -2, 1, 0, 0, 0,      0, 0, 0,               0,  0,  0,  0,  0,  1, 0, 0, 1, 0],
[0, 0, 0, 0, 0,      1, -2, 1, 0, 0,      0, 0, 0,               0,  0,  0,  0,  0,  0, 1, 0, 0, 0],
[0, 0, 0, 0, 0,      0, 1, -2, 1, 0,      0, 0, 0,               0,  0,  0,  0,  0,  0, 0, 0, 0, 0],
[0, 0, 0, 0, 0,      0, 0, 1, -2, 1,      0, 0, 0,               0,  0,  0,  0,  0,  0, 0, 0, 0, 0],
[0, 0, 0, 0, 0,      0, 0, 0, 1, -2,      0, 0, 0,               0,  0,  0,  0,  0,  0, 0, 0, 0, 0],

[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       -2, 1, 0,              0,  0,  0,  0,  0,  0, 1, L3E31, 0, 0],
[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       1, -2, 1,              0,  0,  0,  0,  0,  0, 0, 0, 0, 0],
[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       0, 1, -2,              0,  0,  0,  0,  0,  0, 0, L3E3m1, 0, 0],


[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       0, 0, 0,               -2,  0,  0,  0,  0,  1, 0, 0, 0, 0],

[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       0, 0, 0,               0,  -2,  0,  0,  0,  0, 0, 0, 1, 0],

[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       0, 0, 0,               0,  0,  -2,  0,  0,  0, 0, 0, 0, 1],

[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       0, 0, 0,               0,  0,  0,  -2,  0,  0, 0, 0, 1, 1],

[0, 0, 0, 0, 0,      0, 0, 0, 0, 0,       0, 0, 0,               0,  0,  0,  0,   2,  1, 1, 1, 1, 1],

[1, 0, 0, 0, 0,      1, 0, 0, 0, 0,       0, 0, 0,               1,  0,  0,  0,  1,  -2, 0, 0, 0,  0],
[0, 0, 0, 0, 0,      0, 1, 0, 0, 0,       1, 0, 0,               0,  0,  0,  0,  1,  0, -2, 0, 0,  L2L5],
[0, 1, 0, 0, 0,      0, 0, 0, 0, 0,       L3E31, 0, L3E3m1,      0,  0,  0,  0,  1,  0, 0, -2, L3L4,  0],
[0, 0, 0, 0, 0,      1, 0, 0, 0, 0,       0, 0, 0,               0,  1,  0,  1,  1,  0, 0, L3L4,  -2, 0],
[1, 0, 0, 0, 0,      0, 0, 0, 0, 0,       0, 0, 0,               0,  0,  1,  1,  1,  0,  L2L5, 0, 0,  -2]

]);



assert IsSymmetric(A);
Parent(A);
Rank(A);
rk:=Rank(A);

L0:=Lattice(A);
B0:=Basis(L0);
Sol0:=Solution(A,[Vector(b) : b in B0]);
M0:=Matrix(Sol0);
Gram0:=M0*A*Transpose(M0);
Rank(Gram0);
Factorisation(Determinant(Gram0));
Lambda:=LatticeWithGram(Gram0: CheckPositive:= false);
ALambda,DLambda,phi:=DualQuotient(Lambda);
pSignature(Lambda,-1) eq -18;
ALambda;



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

mE8:=-Matrix(GramMatrix(Lattice("E",8)));
mE8:=ChangeRing(mE8,Integers());
U:=Matrix([[0,1],[1,0]]);
m2:=Matrix([[-2]]);
m4:=Matrix([[-4]]);


Dec1:=<  <mE8,2>, <U,1>, <m2,1>, <m4,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);


function QuadInt(v,w,M)

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

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

end function;


F2:=GF(2);
V2:=VectorSpace(F2,20);
A2:=ChangeRing(Gram0,F2);
k2:=sub<V2 | Kernel(A2)>;
M2:=[v : v in k2 | QuadInt(v,v,Gram0) mod 8 eq 0];
M2;




// Subsection 6.3: s=-1



s:=-1;
K:=Rationals();
L<a>:=ext<K|Polynomial([-s^2+s,0,1])>;
PP<x,y,z,w>:=ProjectiveSpace(L,[1,1,1,3]);
T:=(x^2+y^2)*z-2*x*y*(x+y);
S:=(2*x-z)*(2*y-z)*z;
f1:=T+s*S;
f2:=T+(s+1)*S;
f:=f1*f2;

ll:=[0,2*x, 2*y, x+y, x/(s+a), x/(s-a), y/(s+a), y/(s-a)];
/*
for l in ll do
  l;
  Evaluate(f1*f2,[x,y,l,0]);
  Factorisation($1);
end for;
*/
Ls:=[
[z, w-2*x*y*(x+y)],
[z-2*x, w-2*x^2*(x-y)],
[z-2*y, w-2*y^2*(y-x)],
[z-x-y, w-a*(x-y)^2*(x+y)],

[z-(1-1/s*a)*x, w-x*y*(x-(a+s)*y)/s],
[z-(1+1/s*a)*x, w-x*y*(x+(a-s)*y)/s],
[z-(1-1/s*a)*y, w-x*y*(y-(a+s)*x)/s],
[z-(1+1/s*a)*y, w-x*y*(y+(a-s)*x)/s]

];

for l in Ls do
  l;
  assert Evaluate(f,[x,y,-l[1]+z,w]) eq (l[2]-w)^2;
end for;


ListDivs:=[Scheme(PP,eqn) : eqn in Ls];
SP:=[[1,0,0,0],[0,1,0,0],[1,1,2,0],[1,-1,0,0]];
SP:=[PP!P : P in SP];

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;
	else
	    M[i,j]:=Degree(Int);
	    if SP[1] in Int or SP[2] in Int or SP[3] in Int or SP[4] in Int then
	        M[i,j]:=0;
	    end if;

	end if;
	M[j,i]:=M[i,j];
    end for;
end for;

return M;

end function;

GramMatrixD(ListDivs);




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




mE8:=-Matrix(GramMatrix(Lattice("E",8)));
mE8:=ChangeRing(mE8,Integers());
U:=Matrix([[0,1],[1,0]]);
m2:=Matrix([[-2]]);
m12:=Matrix([[-12]]);


Dec1:=<  <mE8,2>, <U,1>, <m2,1>, <m12,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);

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




for t in CartesianPower([0,1],7) do
L3E31:=t[1];
L3E3m1:=1-L3E31;
L5E22:=t[2];
L5E2m2:=1-L5E22;
L6E22:=t[3];
L6E2m2:=1-L6E22;
L7E12:=t[4];
L7E1m2:=1-L7E12;
L8E12:=t[5];
L8E1m2:=1-L8E12;
L2E21:=t[6];
L2E2m1:=1-L2E21;
L3E11:=t[7];
L3E1m1:=1-L3E11;


A:=Matrix([
[-2, 1, 0, 0, 0,  0, 0, 0, 0, 0,  0, 0, 0,  0,  0,  1, 0, 0, 0,  0, 0, L7E12, L8E12],
[1, -2, 1, 0, 0,  0, 0, 0, 0, 0,  0, 0, 0,  0,  0,  0, 0, L3E11, 0,  0, 0, 0, 0],
[0, 1, -2, 1, 0,  0, 0, 0, 0, 0,  0, 0, 0,  0,  0,  0, 0, 0, 0,  0, 0, 0, 0],
[0, 0, 1, -2, 1,  0, 0, 0, 0, 0,  0, 0, 0,  0,  0,  0, 0, L3E1m1, 0,  0, 0, 0, 0],
[0, 0, 0, 1, -2,  0, 0, 0, 0, 0,  0, 0, 0,  0,  0,  0, 0, 0, 0,  0, 0, L7E1m2, L8E1m2],

[0, 0, 0, 0, 0,  -2, 1, 0, 0, 0,  0, 0, 0,  0,  0,  1, 0, 0, 0,  L5E22, L6E22, 0, 0],
[0, 0, 0, 0, 0,  1, -2, 1, 0, 0,  0, 0, 0,  0,  0,  0, L2E21, 0, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 1, -2, 1, 0,  0, 0, 0,  0,  0,  0, 0, 0, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 0, 1, -2, 1,  0, 0, 0,  0,  0,  0, L2E2m1, 0, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 0, 0, 1, -2,  0, 0, 0,  0,  0,  0, 0, 0, 0,  L5E2m2, L6E2m2, 0, 0],

[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  -2, 1, 0,  0,  0,  0, 1, L3E31, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  1, -2, 1,  0,  0,  0, 0, 0, 1,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  0, 1, -2,  0,  0,  0, 0, L3E3m1, 0,  0, 0, 0, 0],

[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  0, 0, 0,  -2,  0,  1, 0, 0, 1,  0, 0, 0, 0],

[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  0, 0, 0,  0,   2,  1, 1, 1, 1,  1, 1, 1, 1],

[1, 0, 0, 0, 0,  1, 0, 0, 0, 0,  0, 0, 0,  1,  1,  -2, 0, 0, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, L2E21, 0, L2E2m1, 0,  1, 0, 0,  0,  1,  0, -2, 0, 0,  0, 0, 0, 0],
[0, L3E11, 0, L3E1m1, 0,  0, 0, 0, 0, 0,  L3E31, 0, L3E3m1,  0,  1,  0, 0, -2, 0,  0, 0, 0, 0],
[0, 0, 0, 0, 0,  0, 0, 0, 0, 0,  0, 1, 0,  1,  1,  0, 0, 0, -2,  0, 1, 0, 1],

[0, 0, 0, 0, 0,  L5E22, 0, 0, 0, L5E2m2,  0, 0, 0,  0,  1,  0, 0, 0, 0,  -2, 0, 1, 0],
[0, 0, 0, 0, 0,  L6E22, 0, 0, 0, L6E2m2,  0, 0, 0,  0,  1,  0, 0, 0, 1,  0, -2, 0, 1],
[L7E12, 0, 0, 0, L7E1m2,  0, 0, 0, 0, 0,  0, 0, 0,  0,  1,  0, 0, 0, 0,  1, 0, -2, 0],
[L8E12, 0, 0, 0, L8E1m2,  0, 0, 0, 0, 0,  0, 0, 0,  0,  1,  0, 0, 0, 1,  0, 1, 0, -2]


]);

A:=DiagonalJoin(A,Matrix([[-2]]));

assert IsSymmetric(A);
//Parent(A);
rk:=Rank(A);
if rk le 20 then
t;
L0:=Lattice(A);
B0:=Basis(L0);
Sol0:=Solution(A,[Vector(b) : b in B0]);
M0:=Matrix(Sol0);
Gram0:=M0*A*Transpose(M0);
Rank(Gram0);
Factorisation(Determinant(Gram0));
L:=LatticeWithGram(Gram0: CheckPositive:= false);
AL,DL,phi:=DualQuotient(L);
AL;

q:={* Rep2(Norm((phi^-1)(g))) : g in AL *};
assert pSignature(L1,-1) eq pSignature(L,-1);
assert Rank(L1) eq Rank(L);
assert Factorisation(Determinant(L1)) eq Factorisation(Determinant(L));
assert #Generators(AL1) eq #Generators(AL);
assert q eq q1;
assert IsIsomorphic(AL,AL1);

end if;
end for;