//Section 3

//Proposition 3.1
print "Proposition 3.1 verifications.";
F0<A,B>:=FunctionField(Rationals(),2);
F<t>:=FunctionField(F0);
EAB:=EllipticCurve([0,A*t^6+B]);
assert Discriminant(EAB) eq -432*(A*t^6+B)^2;
print "Checked: Discriminant of y^2=x^3+A*t^6+B is -432*(A*t^6+B)^2";

Loc:=LocalInformation(EAB);
print "Places of bad reduction <place,type>: ", [<x[1],x[5]>: x in Loc];

print "Proposition 3.2 verifications.";
//Proposition 3.2
CoordR<a1,a2,a3,a4,a,b,c>:=PolynomialRing(F0,7);

eq1:=[-a4^2 + B + c^3, -2*a3*a4 + 3*b*c^2, -a3^2 - 2*a2*a4 + 3*b^2*c + 3*a*c^2, 
 -2*a2*a3 - 2*a1*a4 + b^3 + 6*a*b*c, -a2^2 - 2*a1*a3 + 3*a*b^2 + 3*a^2*c, -2*a1*a2 + 3*a^2*b, 
 a^3 + A - a1^2];
I:=ideal<CoordR|eq1>;
print "Computing the elimination ideal J...";
J:=EliminationIdeal(I,{a,b,c});
print "Computing the primary decomposition of J...";
pr:=PrimaryDecomposition(J);

Ili:=[GroebnerBasis(Basis(x) cat eq1): x in pr];
Aff:=AffineSpace(CoordR);
Sli:=[Scheme(Aff,g): g in Ili];
 
//orbits sorted by size (like in the paper)
Orbs:=[Sli[i]: i in [8,3,7,1,2,5,6,4]];

print "Orbits with equations and degrees.";

for i in [1..#Orbs] do
	print "Degree of the orbit O_",i," = ",Degree(Orbs[i]);
	print "Equations of the orbit:";
	print DefiningPolynomials(Orbs[i]);
end for;


//verification of the calculations for orbits in Section 3

Q3<zeta3>:=CyclotomicField(3);
Cub2<cub2>:=NumberField(Polynomial([Q3!(-2),0,0,1]));
F2<A6,B6>:=FunctionField(Cub2,2);
//F2<B6>:=FunctionField(F1);
//A:=A6^6;
//B:=B6^6;
FF<t>:=FunctionField(F2);
E:=EllipticCurve([0,A6^6*t^6+B6^6]);
P:=E![-B6^2,A6^3*t^3]; 
Ps:=E![-zeta3*B6^2,A6^3*t^3];
Ps2:=E![-zeta3^2*B6^2,A6^3*t^3]; 
Q:=E![-A6^2*t^2,B6^3]; 
Qs:=E![-zeta3*A6^2*t^2,B6^3];
Qs2:=E![-zeta3^2*A6^2*t^2,B6^3];
s:=cub2^2*A6^2*B6^2;
R:=E![2*A6^3*B6^3/s*t,A6^3*t^3+B6^3];
Rs:=E![2*A6^3*B6^3/(zeta3*s)*t,A6^3*t^3+B6^3]; 
Rs2:=E![2*A6^3*B6^3/(zeta3^2*s)*t,A6^3*t^3+B6^3]; 
S:=E![-2*A6^3*B6^3/s*t,A6^3*t^3-B6^3]; 
Ss:=E![-2*A6^3*B6^3/(zeta3*s)*t,A6^3*t^3-B6^3];
Ss2:=E![-2*A6^3*B6^3/(zeta3^2*s)*t,A6^3*t^3-B6^3];  
bas:=[P,Ps,Q,Qs,R,Rs,S,Ss];

//F0toF2:=hom<F0->F2|[(F2!A6)^6,B6^6]>;
F0toF2:=hom<F0->F2|[A6^6,B6^6]>;

CoordRF2<a1F2,a2F2,a3F2,a4F2,aF2,bF2,cF2>:=PolynomialRing(F2,7);
RtoRF2:=hom<CoordR->CoordRF2|[a1F2,a2F2,a3F2,a4F2,aF2,bF2,cF2]>;

//create a tuple of coordinates in the ring CoordR/CoordRF2 for a given point PP
function ConvertHeight2PointToOrbitCoords(PP)
	x,y,_:=Explode(Eltseq(PP));
	xnum:=Numerator(x);
	ynum:=Numerator(y);
	
	return Reverse([Coefficient(ynum,i): i in [0..3]]) cat Reverse([Coefficient(xnum,i): i in [0..2]]);
end function;

//produce orbits equations in the coordinate ring CoordRF2
OrbsInRF2:=[];
AffRF2:=AffineSpace(CoordRF2);
for orb in Orbs do
	defpols:=DefiningPolynomials(orb);
	tupPol:=[];
	for Pol in defpols do
		coeffs,monoms:=CoefficientsAndMonomials(Pol);
		coeffsInRF2:=[F0toF2(x):x in coeffs];
		monomsInRF2:=[RtoRF2(x):x in monoms];
		Append(~tupPol,&+[coeffsInRF2[i]*monomsInRF2[i]: i in [1..#monomsInRF2]]);
	end for;
	Append(~OrbsInRF2,Scheme(AffRF2,tupPol));
end for;

//orbit 1 computations

print "Verifying Proposition 3.5...";
//points P,Ps, Ps2 and their opposites belonge to the orbit no. 1
orb1ps:=[P,Ps,Ps2,-P,-Ps,-Ps2];
for x in orb1ps do
	assert ConvertHeight2PointToOrbitCoords(x) in OrbsInRF2[1];
end for;

assert IsZero(P+Ps+Ps2);
assert Height(P) eq 2;
assert HeightPairingMatrix([P,Ps]) eq Matrix(Integers(),2,[2,-1,-1,2]);

//orbit 2 computations
print "Verifying Proposition 3.6...";
//points Q,Qs, Qs2 and their opposites belonge to the orbit no. 2
orb2ps:=[Q,Qs,Qs2,-Q,-Qs,-Qs2];
for x in orb2ps do
	assert ConvertHeight2PointToOrbitCoords(x) in OrbsInRF2[2];
end for;

assert IsZero(Q+Qs+Qs2);
assert Height(Q) eq 2;
assert HeightPairingMatrix([Q,Qs]) eq Matrix(Integers(),2,[2,-1,-1,2]);

print "Verifying Proposition 3.8...";
orb3ps:=[R,Rs,Rs2,-R,-Rs,-Rs2,S,Ss,Ss2,-S,-Ss,-Ss2];
for x in orb3ps do
	assert ConvertHeight2PointToOrbitCoords(x) in OrbsInRF2[3];
end for;

assert IsZero(R+Rs+Rs2);
assert Height(R) eq 2;
assert IsZero(S+Ss+Ss2);
assert Height(S) eq 2;
assert HeightPairingMatrix([R,Rs]) eq Matrix(Integers(),2,[2,-1,-1,2]);
assert HeightPairingMatrix([S,Ss]) eq Matrix(Integers(),2,[2,-1,-1,2]);

print "Verifying Proposition 3.9...";
bas:=[P,Ps,Q,Qs,R,Rs,S,Ss];
Mmat:=Matrix(Integers(),2,[2,-1,-1,2]);
assert HeightPairingMatrix(bas) eq KroneckerProduct(IdentityMatrix(Integers(),4),Mmat);

//orbit 4 computations

Uptf:=func<j,k,m|E![-(zeta3^j*(A6^2*cub2^2*t^2 + B6^2*zeta3^k)),(A6*t*(A6^2*cub2^2*t^2 + 2*B6^2*zeta3^k)*(1 + 2*zeta3^m))/cub2^2]>;
//j in {0,1,2}, k in {0,1,2}, m in {1,2}

print "Verifying Proposition 3.10...";
orb4ps:=[Uptf(j,k,m): j,k in [0..2], m in [1..2]];
for x in orb4ps do
	assert ConvertHeight2PointToOrbitCoords(x) in OrbsInRF2[4];
end for;

assert &and[IsZero(Uptf(i,j,1)+Uptf(i,j,2)): i,j in [0,1,2]];
assert &and[IsZero(Uptf(j,0,1)+Uptf(j,1,1)+Uptf(j,2,1)): j in [0,1,2]];
assert &and[IsZero(Uptf(0,k,1)+Uptf(1,k,1)+Uptf(2,k,1)): k in [0,1,2]];

bas4:=[Uptf(0,0,1),Uptf(0,1,1),Uptf(1,0,1),Uptf(1,1,1)];
Mbas4:=HeightPairingMatrix(bas4);
assert IsIsomorphic(LatticeWithGram(Mbas4),Lattice("D",4));
bas4alt:=[-bas4[3],bas4[2],-bas4[1]-bas4[2]-bas4[3],bas4[1]+bas4[2]+bas4[3]+bas4[4]];
Mbas4alt:=HeightPairingMatrix(bas4alt);
IsIsom,iso:=IsIsomorphic(LatticeWithGram(Mbas4alt),Lattice("D",4)); 
assert IsIsom;
assert iso eq IdentityMatrix(Rationals(),4);
//hence orbit no. 4 of length 18 carries 4 linearly independent points which form a D4-type lattice


//orbit 5 computations

Vptf:=func<j,k,m|E![-(zeta3^j*(B6^2*cub2^2 + A6^2*t^2*zeta3^k)),B6*(B6^2 + A6^2*cub2*t^2*zeta3^k)*(1 + 2*zeta3^m)]>;
//j in {0,1,2}, k in {0,1,2}, m in {1,2}

print "Verifying Proposition 3.11...";
orb5ps:=[Vptf(j,k,m): j,k in [0..2], m in [1..2]];
for x in orb5ps do
	assert ConvertHeight2PointToOrbitCoords(x) in OrbsInRF2[5];
end for;

assert &and[IsZero(Vptf(j,k,1)+Vptf(j,k,2)): j,k in [0,1,2]];
assert &and[IsZero(Vptf(j,0,1)+Vptf(j,1,1)+Vptf(j,2,1)): j in [0,1,2]];
assert &and[IsZero(Vptf(0,k,1)+Vptf(1,k,1)+Vptf(2,k,1)): k in [0,1,2]];

bas5:=[Vptf(0,0,1),Vptf(0,1,1),Vptf(1,0,1),Vptf(1,1,1)];
Mbas5:=(HeightPairingMatrix(bas5)); 
assert IsIsomorphic(LatticeWithGram(Mbas5),Lattice("D",4));
bas5alt:=[bas5[1]+bas5[2]+bas5[3]+bas5[4],bas5[2],bas5[4],bas5[3]];
Mbas5alt:=HeightPairingMatrix(bas5alt);
IsIsom,iso:=IsIsomorphic(LatticeWithGram(Mbas5alt),Lattice("D",4)); 
assert IsIsom;
assert iso eq IdentityMatrix(Rationals(),4);

//orbit 6 computations

Wptf:=func<j,k,m,n|E![-(B6^2*cub2^2*zeta3^j) + 2*(-1)^m*A6*B6*cub2*t*zeta3^k + (-1)^m*A6*B6*cub2*t*zeta3^(1 + k) - 2*A6^2*t^2*zeta3^(2*j + 2*k) - 2*A6^2*t^2*zeta3^(1 + 2*j + 2*k),(-1)^n*(B6^3 + 3*(-1)^m*A6^3*t^3 + 2*B6^3*zeta3 - 3*(-1)^m*A6*B6^2*cub2^2*t*zeta3^(1 + 2*j + k) - 
2*A6^2*B6*cub2*t^2*zeta3^(j + 2*k) + 2*A6^2*B6*cub2*t^2*zeta3^(1 + j + 2*k))]>;
//j,k in {0,1,2}; m,n in {0,1}
print "Verifying Proposition 3.12...";
orb6ps:=[Wptf(j,k,m,n): j,k in [0..2], m,n in [0..1]];
for x in orb6ps do
	assert ConvertHeight2PointToOrbitCoords(x) in OrbsInRF2[6];
end for;

assert &and[IsZero(Wptf(j,k,m,0)+Wptf(j,k,m,1)): j,k in [0..2], m in [0,1]];
assert &and[IsZero(&+[Wptf(j,j+s,m,n):j in [0..2]]):s in [0..2],m,n in [0,1]];
assert &and[IsZero(&+[Wptf(j,k,m,n):k in [0..2]]):j in [0..2],m,n in [0,1]];

assert Wptf(1,1,1,0) eq Wptf(0,1,1,0)+Wptf(1,1,0,0)-Wptf(0,1,0,0);
assert Wptf(1,0,1,0) eq Wptf(0,0,1,0)+Wptf(1,0,0,0)-Wptf(0,0,0,0);

bas6:=[Wptf(0,0,0,0),Wptf(0,1,0,0),Wptf(0,0,1,0),Wptf(0,1,1,0),Wptf(1,0,0,0),Wptf(1,1,0,0)];
Mbas6:=HeightPairingMatrix(bas6);
assert IsIsomorphic(LatticeWithGram(Mbas6),Lattice("E",6));
_,iso6:=IsIsomorphic(LatticeWithGram(Mbas6),Lattice("E",6));
bas6alt:=[&+[x[i]*bas6[i]: i in [1..#bas6]]: x in Rows(iso6)];

//orbit 7 computations

Xptf:=func<k,l,m,n|E![zeta3^l*(2*(-1)^n*A6*B6*cub2*t + (-1)^n*A6*B6*cub2*t*zeta3 + 2*B6^2*zeta3^k + 
  A6^2*cub2^2*t^2*zeta3^(2*k) + A6^2*cub2^2*t^2*zeta3^(1 + 2*k)),(-1)^m*(3*B6^3 + (-1)^n*A6^3*t^3 + 2*(-1)^n*A6^3*t^3*zeta3 + 
  3*A6^2*B6*cub2^2*t^2*zeta3^k + 4*(-1)^n*A6*B6^2*cub2*t*zeta3^(2*k) + 3*A6^2*B6*cub2^2*t^2*zeta3^(1 + k) 
  + 2*(-1)^n*A6*B6^2*cub2*t*zeta3^(1 + 2*k))]>;
//j,k in {0,1,2}; m,n in {0,1}
print "Verifying Proposition 3.13...";

orb7ps:=[Xptf(j,k,m,n): j,k in [0..2], m,n in [0..1]];
for x in orb7ps do
	assert ConvertHeight2PointToOrbitCoords(x) in OrbsInRF2[7];
end for;

assert &and[IsZero(Xptf(k,l,0,n)+Xptf(k,l,1,n)): k,l in [0,1,2], n in [0..1]];
assert &and[IsZero(&+[Xptf(k,l,m,0): l in [0..2], m in [1..2]]): k in [0..2]];
assert &and[IsZero(&+[Xptf(k,l,m,0): k in [0..2], m in [1..2]]): l in [0..2]];
assert &and[IsZero(&+[Xptf(j,j+s,m,n): j in [0..2]]): m,n in [0..1],s in [0..2]]; 
assert &and[IsZero(&+[Xptf(j,k,m,n): k in [0..2]]): j in [0..2], m, n in [0..1]]; 

bas7:=[Xptf(0,0,0,0),Xptf(0,1,0,0),Xptf(0,0,0,1),Xptf(0,1,0,1),Xptf(1,0,0,0),Xptf(1,1,0,0)];
Mbas7:=HeightPairingMatrix(bas7);
assert IsIsomorphic(LatticeWithGram(Mbas7),Lattice("E",6));
_,iso7:=IsIsomorphic(LatticeWithGram(Mbas7),Lattice("E",6));
bas7alt:=[&+[x[i]*bas7[i]: i in [1..#bas7]]: x in Rows(iso7)];

//checking the non-degeneracy condition for the bases bas, bas4,bas5,bas6,bas7
print "Checking the non-degenerac condition for the bases generated from the orbits...";
//produces a set of conditions under which the height of the point P might degenerate (differ from the generic value over the function field)
function DegenOfPoint(P)
//at finite places
x:=P[1];
xn:=Numerator(x);
xd:=Denominator(x);
coo:=Coefficients(xn);
cond1:=Monomials(Numerator(coo[#coo]));
cond2:=Monomials(Denominator(coo[#coo]));
r:=Resultant(xn,xd);
cond3:=Monomials(Numerator(r));
cond4:=Monomials(Denominator(r));
condfin:=[cond1,cond2,cond3,cond4];
//at infinity
x:=P[1];
R:=Parent(x);
t:=R.1;
xinf:=Evaluate(x,1/t)*t^2; //x-coordinate at the infinty chart in the integral model

x:=xinf;
xn:=Numerator(x);
xd:=Denominator(x);
coo:=Coefficients(xn);
cond1:=Monomials(Numerator(coo[#coo]));
cond2:=Monomials(Denominator(coo[#coo]));
r:=Resultant(xn,xd);
cond3:=Monomials(Numerator(r));
cond4:=Monomials(Denominator(r));
condinf:=[cond1,cond2,cond3,cond4];
return {x: x in &cat (condfin cat condinf)};
end function;

//function which produces a union of conditions under which the given basis might degenerate
function BasisDegen(ba)
return (&join{DegenOfPoint(ba[i]+ba[j]): i,j in [1..#ba]|i lt j}) join (&join{DegenOfPoint(ba[i]): i in [1..#ba]});
end function;

RR:=RingOfIntegers(F2);
AffRR:=AffineSpace(RR);

//verification of non-degeneracy of the given bases
for ba in [bas,bas4,bas5,bas6,bas7] do
	assert ReducedSubscheme(Scheme(AffRR,[Numerator(&*BasisDegen(ba))])) eq Scheme(AffRR,[Numerator(A6*B6)]); 
end for;

//orbit 8 computations
//initialise the field extensions
Cyc3<zeta3>:=CyclotomicField(3);
K6rel<cub2>:=ext<Cyc3|Polynomial([-2,0,0,1])>;
K6abs:=AbsoluteField(K6rel);
pol:=Polynomial([ 1, 3, 0, -5, 0, 3, 1 ]);
K6:=NumberField(pol);
_,iso:=IsIsomorphic(K6abs,K6);
zeta3:=iso(Cyc3.1);
cub2:=iso(K6rel.1);
Rpol<c1>:=PolynomialRing(K6);
pol:=(4 + 6*c1^2 + c1^3)*(16 - 24*c1^2 + 8*c1^3 + 36*c1^4 - 6*c1^5 + c1^6);
roo:=Roots(pol);
roo:=[x[1]: x in roo];

bpol:=4/3*(7 + 5*c1^3 + 25*c1^6);
bchoice:=[];
for el in roo do
 _,val:=IsPower(Evaluate(bpol,el),3);
 bchoice:=bchoice cat [val];
end for;

function bpow(c1)
_,ret:=IsPower(4/3*(7 + 5*c1^3 + 25*c1^6),6);
return ret;
end function;

a4pol:=1 + c1^3;
a4choice:=[];
for el in roo do
 _,val:=IsPower(Evaluate(a4pol,el),2);
 a4choice:=a4choice cat [val];
end for;

function a4pow(c1)
_,ret:=IsPower(1 + c1^3,2);
return ret;
end function;

//****
F8<A6,B6>:=FunctionField(K6,2);
FF8<t>:=FunctionField(F8);
E:=EllipticCurve([0,A6^6*t^6+B6^6]);
Eptf0:=func<a4,b,c1|E![-(b^2*c1^2*(-11504 + 1544*c1^3 + 7*c1^6))/(5184*B6^2)*t^2+b*t+B6^2*c1,
((a4*b^3*(2728 + 4340*c1^3 + 19*c1^6))/(1296*B6^6))*t^3+((a4*b^2*c1*(4160 + 2524*c1^3 + 11*c1^6))/(2592*B6^4))*t^2+(-(a4*b*c1^2*(-179 + 227*c1^3 + c1^6))/(162*B6^2))*t+a4]>;

function Yptf(o,j,m,n)
roo:=[
    -cub2^2 - 2*cub2 - 2,
    -cub2^2 - 2*zeta3*cub2 + 2*zeta3 + 2,
    -cub2^2 + (2*zeta3 + 2)*cub2 - 2*zeta3,
    -zeta3*cub2^2 - 2*cub2 + 2*zeta3 + 2,
    -zeta3*cub2^2 - 2*zeta3*cub2 - 2*zeta3,
    -zeta3*cub2^2 + (2*zeta3 + 2)*cub2 - 2,
    (zeta3 + 1)*cub2^2 - 2*cub2 - 2*zeta3,
    (zeta3 + 1)*cub2^2 - 2*zeta3*cub2 - 2,
    (zeta3 + 1)*cub2^2 + (2*zeta3 + 2)*cub2 + 2*zeta3 + 2
];

//roogen:=[-cub2^2*(zeta3^(Floor((i-1)/3))+zeta3^(i-1)*cub2^2+zeta3^(2*((i-1)+Floor((i-1)/3)))*cub2): i in [1..9]];
//assert roo eq roogen;

c1:=roo[o];
a4:=B6^3*a4pow(c1)*(-1)^n;

b:=A6*B6*bpow(c1)*zeta3^j*(-1)^m;
return Eptf0(a4,b,c1);
end function;
//o in {1..9}, j in {0,1,2}, m,n in {0,1}

F0toF8:=hom<F0->F8|[(F8!A6)^6,B6^6]>;

CoordRF8<a1F8,a2F8,a3F8,a4F8,aF8,bF8,cF8>:=PolynomialRing(F8,7);
RtoRF8:=hom<CoordR->CoordRF8|[a1F8,a2F8,a3F8,a4F8,aF8,bF8,cF8]>;


//produce orbit 8 equations in the coordinate ring CoordRF8
AffRF8:=AffineSpace(CoordRF8);
orb8:=Orbs[8];
defpols8:=DefiningPolynomials(orb8);
tupPol:=[];
for Pol in defpols8 do
	coeffs,monoms:=CoefficientsAndMonomials(Pol);
	coeffsInRF8:=[F0toF8(x):x in coeffs];
	monomsInRF8:=[RtoRF8(x):x in monoms];
	Append(~tupPol,&+[coeffsInRF8[i]*monomsInRF8[i]: i in [1..#monomsInRF8]]);
end for;
Orb8InRF8:=Scheme(AffRF8,tupPol);

print "Verifying Proposition 3.14...";

orb8ps:=[Yptf(o,k,m,n): o in [1..9],k in [0..2], m,n in [0..1]];
for x in orb8ps do
	assert ConvertHeight2PointToOrbitCoords(x) in Orb8InRF8;
end for;

//relations
assert &and[IsZero(Yptf(o,j,m,0)+Yptf(o,j,m,1)): o in [1..9], j in [1..3], m in [0..1]]; 
assert &and[IsZero(Yptf(1,j,m,0)+Yptf(3,j,m,0)+Yptf(2,j,m+1,0)): j in [0..2],m in [0..1]];
assert &and[IsZero(Yptf(4,j,m,0)+Yptf(5,j,m,0)+Yptf(6,j,m+1,0)): j in [0..2],m in [0..1]];
assert &and[IsZero(Yptf(8,j,m,0)+Yptf(9,j,m,0)+Yptf(7,j,m+1,0)): j in [0..2],m in [0..1]];
assert &and[IsZero(&+[Yptf(o,j,m,0): j in [0..2], m in [0..1]]): o in [1,3,4,5,8,9]];

admtups:={<i,j,m,n>: i in [1..9], j in [0..2], m,n in [0..1]};
admtups:=admtups diff {<i,j,m,1>: i in [1..9], j in [0..2], m in [0..1]};
admtups:=admtups diff {<2,j,m,0>: j in [0..2], m in [0..1]};
admtups:=admtups diff {<6,j,m,0>: j in [0..2], m in [0..1]};
admtups:=admtups diff {<7,j,m,0>: j in [0..2], m in [0..1]};
admtups:=admtups diff {<i,2,1,0>: i in {1,3,4,5,8,9}};
adm1:=Sort([x: x in admtups]);
assert #adm1 eq 30; //relations above reduce the set of 108 points to the set of 30 points

//looking for relations among points corresponding to tuples in the list adm1, which satisfy a relation A+B=C+D, where A,B,C,D in adm1
Eptftup:=func<tup|Yptf(tup[1],tup[2],tup[3],tup[4])>;
pairsall:=[<Eptftup(adm1[i])+Eptftup(adm1[j]),<i,j>>: i,j in [1..30]|i lt j];
relsother:=&cat[[Sort([pairsall[j][2],pairsall[i][2]]): i in [1..#pairsall]|(pairsall[j][1] eq pairsall[i][1]) and not (pairsall[j][2] eq pairsall[i][2])]: j in [1..#pairsall]];
relsother:=Sort([y:y in {x: x in relsother}]);

soo:=Sort([Reverse([x[1][1],x[1][2],x[2][1],x[2][2]]): x in relsother]); //pick element from adm1 (maximal in the ascending order of the tuple which will be eliminated

remain:={i: i in [1..30]} diff {x[1]: x in soo};//remaining indices from the list adm1


ListOfNinePointsInd:=Sort([adm1[j]: j in remain]);
ListOfNinePoints:=[Eptftup(x): x in ListOfNinePointsInd];
MHeight:=HeightPairingMatrix(ListOfNinePoints);
K:=Kernel(MHeight);
SubM:=Submatrix(MHeight,1,1,8,8);
L:=LatticeWithGram(SubM);
assert IsIsomorphic(L,Lattice("E",8));
//so we can eliminate the last element from the list ListOfNinePoints, because it depends linearly on the previous 8, we obtain thus a lattice of rank 8, which we checked is isomorphic to E8 lattice
Rank_8_basis_list_indices:=[ <1, 0, 0, 0>, <1, 0, 1, 0>, <1, 1, 0, 0>, <1, 1, 1, 0>, <1, 2, 0, 0>, <3, 0, 0, 0>, <3, 0, 1, 0>, <4, 0, 0, 0>];
bas8:=[Eptftup(x): x in Rank_8_basis_list_indices];
assert ListOfNinePointsInd[1..8] eq Rank_8_basis_list_indices;

//checking the non-degeneracy condition for the bas8
RR8:=RingOfIntegers(F8);
AffRR8:=AffineSpace(RR8);
assert ReducedSubscheme(Scheme(AffRR,[Numerator(&*BasisDegen(bas8))])) eq Scheme(AffRR,[Numerator(A6*B6)]); 

//Computation of the roots in the spans of orbits
print "Verifying Remark 3.15...";

function ExtractCoefficients(P)
 assert Height(P) eq 2;
 x:=P[1];
 y:=P[2];
 assert Denominator(x) eq 1;
 assert Denominator(y) eq 1;
 assert Degree(Numerator(x)) le 2;
 assert Degree(Numerator(y)) le 3;
 xnum:=Numerator(x);
 ynum:=Numerator(y);
 a:=Coefficient(xnum,2);
 b:=Coefficient(xnum,1);
 c:=Coefficient(xnum,0);
 a1:=Coefficient(ynum,3);
 a2:=Coefficient(ynum,2);
 a3:=Coefficient(ynum,1);
 a4:=Coefficient(ynum,0);
 return [a,b,c,a1,a2,a3,a4];
end function;

//function ExtractOrbitNumber produces the number of the orbit which corresponds to the height 2 point P given on the input
function ExtractOrbitNumber(P)
 ainv:=aInvariants(Scheme(Parent(P)));
 assert &and[ainv[i] eq 0: i in [1..4]];
 a6:=ainv[5];
 assert Denominator(a6) eq 1;
 assert Degree(Numerator(a6)) eq 6;
 coo:=Coefficients(Numerator(a6));
 assert #coo eq 7;
 assert &and[coo[i] eq 0: i in [2,3,4,5,6]];
 A:=coo[7];
 B:=coo[1];
 assert not (A eq 0);
 assert not (B eq 0);
 tup:=ExtractCoefficients(P);
 //tup:=[a,b,c,a1,a2,a3,a4];
 Isa1zero:=func<t|t[4] eq 0>;
 Isa2zero:=func<t|t[5] eq 0>;
 Isa3zero:=func<t|t[6] eq 0>;
 Isa4zero:=func<t|t[7] eq 0>;
 alpha:=func<t|t[3]^3-8*B eq 0>;
 beta:=func<t|t[3]^3+4*B eq 0>;
 
 ans:={1,2,3,4,5,6,7,8};
 if Isa1zero(tup) then
	 if Isa2zero(tup) then
		 return 2;
	 else
		 return 5;
	 end if;
 else
	 if Isa3zero(tup) then
		 if Isa4zero(tup) then
			 return 1;
		 else
			 return 3;
		 end if;
	 else
		 if Isa4zero(tup) then
			 return 4;
		 else
			 if alpha(tup) then
				 return 7;
			 else
				 if beta(tup) then
					 return 6;
				 else
					 return 8;
				 end if;
			 end if;
		 end if;
	 end if;
 end if;
 
 print "error";
 
 assert false;
end function;
		  
 
//roots in the span of orbit: 4
D4:=Lattice("D",4);
svD4:=ShortestVectors(D4);
VD4,mVD4:=CoordinateSpace(D4);
liO4:=[mVD4(svD4[i]): i in [1..#svD4]] cat [mVD4(-svD4[i]): i in [1..#svD4]];
rooO4:=[ExtractOrbitNumber(&+[Integers()!x[i]*bas4alt[i]: i in [1..4]]): x in liO4];
assert #rooO4 eq 24;
assert &+[1:x in rooO4| x eq 4] eq 18;
assert &+[1:x in rooO4| x eq 1] eq 6;

 //roots in the span of orbit: 5
 D4:=Lattice("D",4);
 svD4:=ShortestVectors(D4);
 VD4,mVD4:=CoordinateSpace(D4);
 liO5:=[mVD4(svD4[i]): i in [1..#svD4]] cat [mVD4(-svD4[i]): i in [1..#svD4]];
 rooO5:=[ExtractOrbitNumber(&+[Integers()!x[i]*bas5alt[i]: i in [1..4]]): x in liO5];
 assert #rooO5 eq 24;
 assert &+[1:x in rooO5| x eq 5] eq 18;
 assert &+[1:x in rooO5| x eq 2] eq 6;
 
 //roots in the span of orbit: 6
 E6:=Lattice("E",6);
 svE6:=ShortestVectors(E6);
 VE6,mVE6:=CoordinateSpace(E6);
 liO6:=[mVE6(svE6[i]): i in [1..#svE6]] cat [mVE6(-svE6[i]): i in [1..#svE6]];
 rooO6:=[ExtractOrbitNumber(&+[Integers()!x[i]*bas6alt[i]: i in [1..6]]): x in liO6];
 assert #rooO6 eq 72;
 assert &+[1:x in rooO6| x eq 6] eq 36;
 assert &+[1:x in rooO6| x eq 5] eq 18;
 assert &+[1:x in rooO6| x eq 3] eq 12;
 assert &+[1:x in rooO6| x eq 2] eq 6; 
 
 //roots in the span of orbit: 7
 E6:=Lattice("E",6);
 svE6:=ShortestVectors(E6);
 VE6,mVE6:=CoordinateSpace(E6);
 liO7:=[mVE6(svE6[i]): i in [1..#svE6]] cat [mVE6(-svE6[i]): i in [1..#svE6]];
 rooO7:=[ExtractOrbitNumber(&+[Integers()!x[i]*bas7alt[i]: i in [1..6]]): x in liO7];
 assert #rooO7 eq 72;
 assert &+[1:x in rooO7| x eq 7] eq 36;
 assert &+[1:x in rooO7| x eq 4] eq 18;
 assert &+[1:x in rooO7| x eq 3] eq 12;
 assert &+[1:x in rooO7| x eq 1] eq 6;    
 

 print "Section 3 verified!"