
// Requires excecution of "main.m"

//************************************************************************
//       Utility functions.
//************************************************************************

// Takes a collection of cubics L defined over a (base extension of) P3.
// Returns a subspace of the space of cubics.
baseCSvectors := { Vector([MonomialCoefficient(QC*P3.i ,mon) : mon in Bcubics]) 
		 : i in [1 .. 4]};
function CayleySubspace(L)
    assert #L eq 4;
    R := Universe(L);
    K := BaseRing(R);
    BcubicsK := [ R ! mon : mon in Bcubics];
    VcubicsK := VectorSpace(K, 20);

    // The field of definition of each Cayley subspace is different.
    // We change the ring of the part defined by the quadric.
    csvectors := ChangeUniverse(baseCSvectors, VectorSpace(K,20));
    for f in L do
	Include(~csvectors, Vector([MonomialCoefficient(f,mon) 
				    : mon in BcubicsK]));
    end for;
    return sub<VcubicsK | csvectors>;
end function;

//************************************************************************
// New idea. Attempt to recover Cayley cubic via intersections
//************************************************************************

// DEPRECIATED. 
// Takes a *list* of Cayley spaces *over different ground fields*.
// Returns a cubic of some kind defining the curve. 
/*
function IntersectCayleySpaces2(cayleyList)

    // The correct subspace is known to be defined over k, so we
    // only look for k-subspaces contained in the Cayley spaces.
    K := BaseRing(cayleyList[1]);
    V := cayleyList[1];

    if not Type(K) eq FldFin then
	error "Not implemented for spaces over", Type(K);
    end if;

    for cs in cayleyList do
	K := CommonOverfield(K,BaseRing(cs));
	V := ChangeRing(V,K) meet ChangeRing(cs,K);
    end for;

    return V;
end function;
*/

// Takes a *list* of Cayley spaces *over different ground fields*.
// Returns a cubic of some kind defining the curve. 
function IntersectCayleySpaces(cayleyList)

    // The correct subspace is known to be defined over k, so we
    // only look for k-subspaces contained in the Cayley spaces.
    V := VectorSpace(k,20);
    for cs in cayleyList do
	V := V meet RestrictField(cs,k);
    end for;
    return V;
end function;



// **********************************************************************

// The dual to the projective space containing the space sextic, and
// a dummy factor for the product space.
P3dual<a0,a1,a2,a3> := ProjectiveSpace(k,3);
P3Y<y0,y1,y2,y3> := ProjectiveSpace(k,3);

// Product space and multiplication map.
// Multiplication takes two degree 1 forms and produces the degree 2 form 
// with the variables Xi identified with Yi.
KCxKC<X0,X1,X2,X3,Y0,Y1,Y2,Y3>, pi := DirectProduct(P3dual, P3Y);
bimons := 
[X0*Y0, 
X0*Y1 + X1*Y0, 
X1*Y1,
X0*Y2 + X2*Y0, 
X1*Y2 + X2*Y1,
X2*Y2,
X0*Y3 + X3*Y0,
X1*Y3 + X3*Y1,
X2*Y3 + X3*Y2,
X3*Y3];

// The space of quadrics.
VC := ProjectiveSpace(k,9);
mult := map< KCxKC->VC | bimons>;

// Given a subspace W of a 10 dimensional vector space
// return W as a subscheme of VC
function SubspaceToSubscheme(W)
	B := Matrix(Basis(W));
	ker := Basis(NullspaceOfTranspose(B));
	eqns := [&+[ VC.i*ker[j][i] : i in [1 .. 10]] : j in [1 ..#ker]];
	
	return Scheme(VC, eqns);
end function;


//************************************************************************
//       Main procedure.
//************************************************************************

curveList := [];
cayleyList := [* *];

for W in steinerSubspaces do

	Wsch := SubspaceToSubscheme(W);
	Xatatmult := Wsch @@ mult;
	Z := Xatatmult; 	            // Alias for easier typing
	
	// Ensure that the system is non-degenerate.
	assert Dimension(Xatatmult) eq 3;

	// Compute the fibres of the projection onto the first coordinate
	// The Scheme of points lying in fibres isomorphic to P2 is exactly the
	// Scheme of (nonsingular?) points such that the relative normal bundle has rank 1.
	// Computationally this is much easier done than said.

	JMat := JacobianMatrix(Z);
	JMaty:= ColumnSubmatrix(JMat, 5, 4);

	// Measures if projection fails to be invertible here.
	// In fact, we measure if a 2-dimensional subspace of the tangent bundle is 
	// killed by projection.
	eqns := DefiningEquations(Z) cat Minors(JMaty,2); 

	ExPlanes := Scheme(KCxKC, eqns);
	assert Dimension(ExPlanes) eq 2;

	// Project the image of this scheme to get the four hyperplanes through the 
	// four nodes of the Cayley cubic system.
	fourplanes := pi[1](ExPlanes);
	assert Dimension(fourplanes) eq 0;
	assert Degree(fourplanes) eq 4;
	assert IsReduced(fourplanes);

	/////////////////////////////////////////////////////////////////////////////////
	// Based on the various assert statements, I assume that everything has worked 
	// correctly up to this point. 
	/////////////////////////////////////////////////////////////////////////////////

	// We now compute the Cayley cubic system. 
	_, clstr := IsCluster(fourplanes);
	fourplanespts := PointsOverSplittingField(clstr);
	fldofdef := Ring(Universe(fourplanespts));

	// Abstract generators for the Cayley subspace
	MyRng<T0,T1,T2,T3> := PolynomialRing(fldofdef,4);
	genericCSgens := [ MyRng ! (T0*T1*T2*T3/MyRng.i) : i in [1 .. 4]];

	function DualPointToPlane(P)
		L := Eltseq(P);
		return &+[ MyRng.i*(fldofdef ! L[i]) : i in [1 .. 4]];
	end function;	

	h := [DualPointToPlane(p) : p in fourplanespts];
	csGens := [Evaluate(g, h) : g in genericCSgens];        

	// Add the Cayley subspace to the list;
	print "Subspace", #cayleyList + 1, "computed.";
	Append(~cayleyList, CayleySubspace(csGens));

end for;
