
load "genericTritangents.m";
load "steinerSolver.m";

// **********************************************************************
// The code assumes that the space sextic lies on the standard quadric
// for ease of computing the tritangents. However, the recovery
// functions should work generically.

// **********************************************************************
// Definition of ground projective space and declared basis elements.
// Definition of the standard quadric is also included here.

kbase := GF(7);
P3base<x0,x1,x2,x3> := ProjectiveSpace(kbase,3);

// **********************************************************************
// Curve defined.

// The smooth quadric (DO NOT CHANGE)
q := x0*x2-x1*x3;

// A Random cubic.
// f := x0^3 + x1^3 + x2^3 - 101*x2*x1*x0 + 432*x3^2*(x0+x1) + x1*x2^2;

// A transformed Cayley cubic.
f := x0*x1*x2 + x0*x1*x3 + x0*x2*x3 + x1*x2*x3;
f := 1/4*x0*x1*x2 + x0*x1*x3 + x0*x2*x3 + x1*x2*x3;
f := Evaluate(f, [x0-x1, x0-2*x2+x1, x2-x3 + 5*x0, x3-4*x1 + 11*x0]);

// The curve.
C := Curve(Scheme(P3base, [q,f]));
assert IsNonSingular(C);


k, P3, C, Tgt_list := ComputeTritangents(C);
steiner := ComputeSteinerSetsBySyzygy(C, Tgt_list);

// The other option, which takes much longer.
// time steiner := ComputeSteinerSetsByLinearEquivalence(C, Tgt_list);

// **********************************************************************
// The ambient projective space has been defined, so we define new
// Utility functions

// CONSTANTS ************************************************************
Vquadrics := VectorSpace(k,10);
Bquadrics := Monomials(CoordinateRing(P3) ! (x0+x1+x2+x3)^2);
Bcubics := Monomials(CoordinateRing(P3) ! (x0+x1+x2+x3)^3);

// FUNCTIONS  ***********************************************************

// Compute the Steiner subspaces
function SteinerSubspace(ss)
  quadraticMonomials := [CoordinateRing(P3) ! mon : mon in Bquadrics];
  spvectors := {};
  for sp in ss do
	  spQuadric := &*sp;
	  Include(~spvectors, Vector([MonomialCoefficient(spQuadric,mon)
				                        : mon in quadraticMonomials]));
  end for;
  return sub<Vquadrics | spvectors>;
end function;

// Given a vector, return a quadric equation of P3
function VectorToQuadric(v)
  R := CoordinateRing(P3);
  return &+[ (R ! v[i])*Bquadrics[i] : i in [1 .. #Bquadrics]];
end function;

// Given a vector, return a cubic equation of P3
function VectorToCubic(v)
  R := CoordinateRing(P3);
  return &+[ (R ! v[i])*Bcubics[i] : i in [1 .. #Bcubics]];
end function;

// Tests if the list of linear forms ``Tgt_list'' are the equations
// defining the 120 tritangents of the curve ``crv''.
//
// INPUTS:
// crv -- A space sextic in P3
// Tgt_list -- A list of 120 (unlabelled) linear forms.
//
// OUTPUTS:
// returns ``true'' if Tgt_list are 120 planes which are tangent to
// ``crv'' at each point of contact.
// Also return a list of 120 booleans, one for each form in Tgt_list.
function testTgt(crv, Tgt_list)
  L := [];
  for tgt in Tgt_list do
	  Z := Scheme(P3,tgt) meet crv;
	  Sing := SingularSubscheme(Z);
	  Red  := ReducedSubscheme(Z);

	  // We check if every point in the intersection is singular.
	  // We then check to make sure that (crv meet tgt) does not
	  // consist of 2 points, each of multiplicity 3.
	  if not IsSubscheme(Red, Sing) then
	    Append(~L, false);
	    continue;
	  elif Degree(Red) ge 3 or Degree(Red) eq 1 then
	    Append(~L, true);
	    continue;
    end if;

	  comps := IrreducibleComponents(Sing);
	  // The only intersection types which can occur at this point
	  // are (2,4) or (3,3). The individual points may not be
	  // defined over the ground field, hence the indirect checks.
	  // (Note: We use the singular subscheme to count multiplicities)
	  if Degree(Red) eq 2 and #comps eq 1 then
	    Append(~L, false);
    else
      Append(~L, &and[IsOdd(Degree(co)) : co in comps]);
	  end if;
  end for;
  return &and L, L;
end function;

// **********************************************************************
// Resume main procedure.
// **********************************************************************

steinerSubspaces := [SteinerSubspace(ss) : ss in steiner];

print "intersecting steiner subspaces";
Vquadric := &meet steinerSubspaces;
assert Dimension(Vquadric) eq 1;

print "basis of quadrics:";
print Bquadrics;
print "intersection of steiner subspaces:";
print Vquadric;

// Quadric equation defined.
QC := VectorToQuadric(Basis(Vquadric)[1]);

print "Quadric equation:";
print QC;

//************************************************************************
// Recover the cayley cubic.
//************************************************************************
load "cayleySpaceCode.m";

//************************************************************************
// Recover the curve.
//************************************************************************

Vcurve:= IntersectCayleySpaces(cayleyList);
defeq := [ VectorToCubic(v) : v in Basis(Vcurve)];
Ctest := Curve(Scheme(P3, defeq));

print "Check if the reconstruction was successful:",  C eq Ctest;
