/* ****************** LinesCalculated ********************* */

/* This MAGMA script gives the equations of the 240 lines on a Del Pezzo
surface of degree 1. It assumes that the file CoeffEqs has been
loaded in MAGMA. More specifically, we compute the lines of

	w^2 = x^6 + y^6 + z^3.

Theory shows they are of the form z = a*x^2 + b*x*y + c^y^2, w = \pm r*x^3 + 
s*x^2*y + t*x*y^2 + u*y^3. This file computes arrays Zvalues and Wvalues,
and the equation of any line on the surface is of the form z = 
Zvalues[i], w = Wvalues[i], where i = 1,...,240.. */

/* For each of the families of coefficients computed by the script "lines", 
the array w(i)zs below stores the quadratic forms that we get for z. The 
numbering below does not correspond to the numbering in "lines". However, 
they are identified as follows:

      family # in lines2           in this script
      1,2,3,4,5,6,7,8        =     (51,52),4,2,5,3,6,7,8

Note family 1 in lines2 had to be split up into two families here (to make 
life easier). */

w1zs := [ (r1[j])*x^2 + (p1[i])*x*y + (q1[j])*y^2 
                : j in [(i-1)*#p1 + 1..i*#p1], i in [1..#p1] ]; 
w2zs := [ (r2[i])*x^2 + (p2[i])*x*y + (q2[i])*y^2 
                : i in [1..#p2] ]; 
w3zs := [ (r3[j])*x^2 + (p3[i])*x*y + (q3[j])*y^2 
                : j in [(i-1)*#p3 + 1..i*#p3], i in [1..#p3] ]; 
w4zs := [ (r4[i])*x^2 + (p4[i])*x*y + (q4[i])*y^2 
                : i in [1..#p4] ]; 
w51zs := [ (r5[j])*x^2 + (q5[i])*y^2 
                : j in [(i-1)*3 + 1..i*3], i in [1..4] ]; 
w52zs := [ (r5[j])*x^2 + (q5[i])*y^2 
                : j in [(i-2)*4 + 1..(i-1)*4], i in [5..7] ];
w6zs := [ (r6[(i-1)*#q6 + k])*x^2 + (p6[i])*x*y + (q6[k])*y^2 
                : k in [1..#q6], i in [1..#p6] ]; 
w7zs := [ (r7[j])*x^2 + (p7[i])*x*y + (q7[j])*y^2 
                : j in [(i-1)*3 + 1..i*3], i in [1..#p7] ]; 
w8zs := [ (r8[j])*x^2 + (p8[i])*x*y + (q8[j])*y^2 
                : j in [(i-1)*3 + 1..i*3], i in [1..#p8] ];

/* Now put all possible quadratic forms for z together in one array indexed 
by family */

zvalues := [ w1zs, w2zs, w3zs, w4zs, w51zs, w52zs, w6zs, w7zs, w8zs ];
wvalues := [];

/* The following for loop computes the equation w = cubic form for each 
quadratic form in zvalues. The result is stored in the array wvalues, which
is indexed the same way zvalues is indexed. */

for j in [1..#zvalues] do
	partwvalues := [];
	for i in [1..#zvalues[j]] do
		f, u := Factorization(x^6 + y^6 + (zvalues[j][i])^3);
       		A := Factorization(x^2 - u);
       		unit := -Reductum(A[1][1]);
		if f[1][2] ne 2 then
			Append(~partwvalues, unit*(&*[factr[1] : factr in f])^3);
		else
	       		Append(~partwvalues, unit*(&*[factr[1] : factr in f]));
		end if;
	end for;
	Append(~wvalues, partwvalues);
end for;

/* The "if" in the middle of the for-loop is designed for cases when f = 
<x,6> or f = <y,6>. In that case we want w = x^3 or w = y^3 and not w = x or 
w = y. */

// If you want to store the results in a file you can do this:
/* 
FL := Open("LinesFile", "w");
Put(FL, "zvalues := ");
Write("LinesFile", zvalues);
Puts(FL, ";");
Put(FL, "wvalues := ");
Write("LinesFile", wvalues);
Puts(FL, ";");
*/

/* Now we check that we do indeed have the equations of a line */
/*
Iamthechampion := { };
for i in [1..#zvalues] do
	for j in [1..#zvalues[i]] do
		ans := wvalues[i][j]^2 - x^6 - y^6 - zvalues[i][j]^3;
		if ans ne 0 then
			print "We have a Problem...";
		end if;
		Include(~Iamthechampion, ans);
	end for;
end for;
print "I am the Champion.";
print Iamthechampion;
*/

// Now we put all of zvalues into a single array with 120 elements.
// Ditto for wvalues.

pZvalues := [ zvalues[i][j] : j in [1..#zvalues[i]], i in [1..#zvalues] ];
pWvalues := [ wvalues[i][j] : j in [1..#wvalues[i]], i in [1..#wvalues] ];

// Finally, we duplicate pZvalues to get 240 elements.
// We also put pWvalues and -pWvalues into an array of 240 elements.

Zvalues := [ pZvalues[i] : i in [1..#pZvalues], j in [1..2] ];
Wvalues := [ ((-1)^(k+1))*pWvalues[i] : i in [1..#pWvalues], k in [1..2] ];

// Finally, just for kicks, we assemble the lines:

LINES := [ [Zvalues[i], Wvalues[i] ] : i in [1..240] ];

// Write("LinesFile", LINES);

/* ****************** Intersection ********************* */

/* The following script calculates the 240x240 matrix of intersections of the 
exceptional curves on the surface w^2 = z^3 + x^6 + y^6. It assumes the file 
LinesCalculated has been loaded. */

P1<xx,yy> := ProjectiveSpace(Q,1);

IntEqs := function(f,g)
	f1 := Evaluate(f,[xx,yy,0,0,0,0,0]);
	g1 := Evaluate(g,[xx,yy,0,0,0,0,0]);
	S := Scheme(P1,[f1,g1]);
	return Degree(S);
end function;

Intersection := function(i,j)
        if i eq j then
                return -1;
        end if;
	f := Zvalues[i] - Zvalues[j];
	g := Wvalues[i] - Wvalues[j];
	return IntEqs(f,g);
end function;

Mat := MatrixRing(IntegerRing(), 240);
time IntersectionMatrix := Mat![ Intersection(i,j) : j in [1..240], i in [1..240]];
// print IntersectionMatrix;
// Write("Matrix", IntersectionMatrix);