
// File for generic test routines.
// WARNING: The code is very slow, and meant to test the geometry over a (small) finite field.

ZZ := Integers();
P6 := ProjectiveSpace(ZZ,6);
c := ProjectiveSpace(ZZ,3);
_<t> := PolynomialRing(CoordinateRing(c));
g := &+[c.i*t^(i-1) : i in [1 .. 4]];

square_map := map< c->P6 | Coefficients(g^2)>;
VarSqrsZZ := Image(square_map);


// Computes the tritangents to a generic space sextic curve C
// The code fails on a small percentage of inputs. I have no
// idea why, but it is not a priority to figure this out.
function ComputeTritangents(C)

    k := BaseField(C);
    VarSqrs := ChangeRing(VarSqrsZZ,k);
    P3<x0,x1,x2,x3> := AmbientSpace(C);
    de := DefiningEquations(C);
    q := [e : e in de | Degree(e) eq 2][1];
    f := [e : e in de | Degree(e) eq 3][1];

//    P6k := ChangeRing(P6,k);
    P6k := AmbientSpace(VarSqrs);

    // Assert the quadric is normalized
    assert q eq x0*x2-x1*x3;

    // Solve for the tritangents
    R1<s,t,u,v> := PolynomialRing(k,4);
    phi1 := hom< CoordinateRing(P3) ->  R1 | [s*u, t*u, t*v, s*v]>;

    R2<a,b,c,d> := PolynomialRing(k,4);
    P3dual := Proj(R2);
    R3<S,U> := PolynomialRing(R2,2);
    phi2 := hom<R1 -> R3 | [S,1,U,1]>;

    F := phi2(phi1(f));
    h := a*S*U + b*U + c + d*S;

    Fs := Resultant(F,h,U);
    maptoS := map<P3dual -> P6k | [Coefficient(Fs, S, i) : i in [0 .. 6]]>; 

    Fu := Resultant(F,h,S);
    maptoU := map<P3dual -> P6k | [Coefficient(Fu, U, i) : i in [0 .. 6]]>; 

    Tgt_scheme := (VarSqrs @@ maptoS) meet (VarSqrs @@ maptoU);
    _, Tgt_scheme := IsCluster(Tgt_scheme);

    // Sometimes the base field causes problems. I'm not sure why.
    assert Degree(Tgt_scheme) eq 120;
    Tgt_pts := PointsOverSplittingField(Tgt_scheme);

    // **********************************************************************
    // Base field, curve, and ambient projective space redefined.
    // This is to be consitent with other code that has been written.
    // It may not be a bad idea to discuss the design of the final code.

    // Modify the base field of objects.
    k := Parent(Random(Tgt_pts)[1]);
    P3<x0,x1,x2,x3> := ChangeRing(P3, k);
    C := Curve(Scheme(P3, [f,q]));
    Tgt_list := [ p[1]*x0 + p[2]*x1 + p[3]*x2 + p[4]*x3 : p in Tgt_pts];

    return k, P3, C, Tgt_list;
end function;

// A test function to ensure that the tritangents have been computed
// correctly. Assumes that the curve and tritangents share the same
// field of definition.
function testTgtNoLabels(crv, Tgt_list)
    P3 := AmbientSpace(crv);
    L := [];
    for tgt in Tgt_list do
        L cat:= [Degree(ReducedSubscheme(Scheme(P3,tgt) meet crv))];
    end for;
    return L;
end function;

