//////////////////////////////////////////////////////////////////////
//
// ellfib.mag
// GBD, DJR 19 June 2008
//
// Online at:
//	http://www.kent.ac.uk/ims/personal/gdb/ellfib.mag
// A file of examples is at
//	http://www.kent.ac.uk/ims/personal/gdb/examples.ellfib
// 
// 
//	EllFib package - elliptic fibrations on cubic surfaces
//	======================================================
//
// 
// Suppose $EF is the path to the directory containing this file.
// Load the EllFib package in a magma session with the line:
// 
// Attach("$EF/ellfib.mag");
// 
//////////////////////////////////////////////////////////////////////

/*
List of intrinsics
==================

This file has the following intrinsics (other functions are
auxilliary to these):

GeiserInvolution(X::Sch,p::Pt) -> MapSch
  {Let X be a cubic surface in P3 and p a nonsingular point of X.
  Return the Geiser involution X -> X centred in p.
  (X is not required to be minimal - the usual requirement - but
  in that case the Geiser involution may not be defined
  and an error will be reported.)}

BertiniInvolution(X::Sch,Z::Sch) -> MapAutSch
  {Let X be a cubic surface in P3 and Z a finite subscheme of X of length 2.
  Return the Bertini involution X -> X centred in Z.}

ImposeBasepoint(X::Sch,L::LinearSys,p::Pt,m::RngIntElt) -> LinearSys
ImposeBasepoint(X::Sch,L::LinearSys,Z::Sch,m::RngIntElt) -> LinearSys
  {Let X be a surface in projective 3-space P3, L a linear system on P3
  and p a nonsingular point of X (or Z a subscheme of degree 2).
  Return the sublinear system of L of forms that vanish at p (or Z) to
  order at least m.}

Two supporting intrinsics:

Degree2SplittingField(Z::Sch) -> Rng
  {Let Z be a zero-dimensional scheme of degree 2 in P^n over a field k.
  Return an extension of k (of deg 1 or 2) over which Z splits geometrically.}

VanishingConditions(X::Sch,p::Pt,d::RngIntElt,m::RngIntElt) -> Mtrx
  {Let X be a cubic surface in projective 3-space P3 and p be a nonsingular
  point of X.
  Return the matrix of coefficients of linear forms of degree d on P3
  that vanish to order m at p when restricted to X.}
*/


//////////////////////////////////////////////////////////////////////
//
//		Main vanishing calculation
//
//////////////////////////////////////////////////////////////////////


// Input: X a surface in P3, p a nonsingular rational point of X.
// Output: a blowup map A3 -> P3 where A3 = A3<x,y,z> and
//      E_big = (z = 0), E_X = X~ \cap E_big.
// Also, to first order (i.e. the tangent space)
//      E_big = (z = 0), E_X = (y=0 in E_big) = (x-axis).
//
function good_blowup_patch(X,p)
        assert IsNonsingular(X!p);
        k := BaseField(X);
        P3 := Ambient(X);
        // First move p to (0001).
        phi1 := Translation(X,p);
        // following 3 lines replace the natural line
        // X1 := phi1(X);
        // which failed for a while in magma v2.14
        f :=DefiningEquation(X);
        am_phi1 := AlgebraMap(phi1^-1);
        X1 := Scheme(Ambient(X),f @ am_phi1);
        p1 := phi1(p);
        f1 := Equation(X1);     // phi1(p) = (0:0:0:1) in (f1=0).
        T1 := TangentSpace(X1,p1);
        t1 := Equation(T1);
        coeffs1 := [ MonomialCoefficient(t1,P3.i) : i in [1..4] ];
        // Second move T1 to y=0 at p1=(0001).
        // Please do not laugh at us.
        M := Transpose(Matrix(4,[1,0,0,0] cat coeffs1 cat [0,0,1,0,0,0,0,1]));
        if Determinant(M) eq 0 then
                M[2,1] := 1;
        end if;
        if Determinant(M) eq 0 then
                M[2,1] := 0;
                M[2,3] := 1;
        end if;
        assert Determinant(M) ne 0;
        // following 6 lines replace the natural lines
        // phi2 := Automorphism(P3,M);
        // assert phi2(T1) eq Scheme(P3,P3.2);
        // which failed for a while in magma v2.14
        eqns := [ &+[ M[j,i] * P3.j : j in [1..4]] : i in [1..4]];
        phi2 := map< P3 -> P3 | eqns >;
        am_phi2 := AlgebraMap(phi2^-1);
        f := DefiningEquation(T1);
        T2 := Scheme(Ambient(T1),f @ am_phi2);
        assert T2 eq Scheme(P3,P3.2);

        // Third blowup using the z patch.  i.e. using xz,yz,z.
        A3<xx,yy,zz> := AffineSpace(k,3);
        xx := A3.1; yy := A3.2; zz := A3.3;
        phi3 := map< A3 -> P3 | [xx*zz, yy*zz, zz, 1] >;
        return phi3 * Inverse(phi2) * Inverse(phi1);
end function;

// Input: X a surface in P3, p in X a nonsingular rational point,
//      m a positive integer.
// Output: In the previous notation, compute y as a function in k(x)[[z]].
// NB. Of course m will be a multiplicity in the end, but here its
// role is only to decide what precision to expand a power series out to.
//
function local_along_E(X,p,m)
        // blowup X at p (in lovely previous coords)
        phi := good_blowup_patch(X,p);
        A3 := Domain(phi);
        R3 := CoordinateRing(A3);
        Xtotalblowup := X @@ phi;
        f_tot := DefiningEquation(Xtotalblowup);
        bool,f_bir := IsDivisibleBy(f_tot, A3.3);
        assert bool and not IsDivisibleBy(f_tot, A3.3^2);
        // make the magic ring k(x)[[z]][y].
        k := BaseField(X);
        K := RationalFunctionField(k);  xxx := K.1;
        P := PowerSeriesRing(K: Precision:=m);  zzz := P.1;
        R := PolynomialRing(P); yyy := R.1;
        // map our poly ring k[x,y,z] to k(x)[[z]][y].
        psi := hom< R3 -> R | f :-> Evaluate(f,[xxx,yyy,zzz]) >;
        // solve, finding a root with zero constant term.
        a := ImplicitFunction(psi(f_bir),1,m);  // Expand to order m.
        return a, psi, phi;
end function;

intrinsic VanishingConditions(X::Sch,p::Pt,d::RngIntElt,m::RngIntElt) -> Mtrx
{Let X be a cubic surface in projective 3-space P3 and p be a nonsingular
point of X.
Return the matrix of coefficients of linear forms of degree d on P3
that vanish to order m at p when restricted to X.}
        k := BaseField(X);
        P3 := Ambient(X);
        P3coords := CoordinateRing(P3);
        // Find y as implicit function of z at generic point along x-axis.
        // Notation: Y is y(x).  Also we get blowup patch phi here.
        Y,psi,phi := local_along_E(X,p,m);
        A3 := Domain(phi);
        A3coords := CoordinateRing(A3);
        assert Domain(psi) eq A3coords;
        phi_hash := hom< P3coords -> A3coords | DefiningEquations(phi) >;

        // Now build a generic d-tic pulled back to the generic point
        // of the exceptional curve of phi_{|X}, localised there.
        N := Binomial(4+d-1,d);
        monos_on_P3 := MonomialsOfDegree(P3coords,d);
        monos_lifted_to_local := [ psi(phi_hash(M)) : M in monos_on_P3 ];
        local_ring := Codomain(psi);
        local_ring_unknowns<[a]> := PolynomialRing(local_ring,N);
        // pull generic d-tic back, evaluating y at Y as we go.
        F := &+ [ a[i] * Evaluate(monos_lifted_to_local[i],Y) : i in [1..N] ];
        // get coefficients of F(x,Y,z) up to order z^(m-1).
        coeffs_of_z := [
        &+[ a[j] * Coefficient(Coefficient(MonomialCoefficient(F,a[j]),0),i) : j in [1..N]]
                : i in [0..m-1] ];

        // build matrix of conditions.
        rows_of_a_matrix := [];
        for e in [0..m-1] do
                coeff_to_kill := coeffs_of_z[e+1];
                new_row :=
                [ MonomialCoefficient(coeff_to_kill,a[i]) : i in [1..N] ];
                Append(~rows_of_a_matrix,new_row);
        end for;
        // we know these coeffs are elts of the basering, so coerce
        // back down to that.
        all_entries := &cat rows_of_a_matrix;
        // these entries will all lie in k(x).
        H<s> := RationalFunctionField(k);
        ChangeUniverse(~all_entries, H);
        Moverkx := Matrix(N,all_entries);
        // cool.
        // Each row may have x's in it; split this into separate
        // conditions, each degree of x at a time.
        entries_over_k := [ k | ];
        biggest_power_of_x := Maximum([Degree(p): p in Eltseq(Moverkx)]);
        for i in [1..NumberOfRows(Moverkx)] do
                row := Eltseq(Moverkx[i]);
                // We believe magma has rigged the rational functions in x
                // to have denom=1 - what we need, in fact, is only that they
                // have precisely the same denominator. We make the following
                // sufficient assert, although could do better if it failed:
                assert &and [ Denominator(q) eq 1 : q in row ];
                rowdata := [ [Coefficient(Numerator(q),i) :
                                i in [0..biggest_power_of_x]] : q in row];
                entries_over_k cat:= &cat[ [ rowdata[i][j] :
                                i in [1..N] ] : j in [1..#rowdata[1]] ];
        end for;
        Moverk := Matrix(N,entries_over_k);
        // remove zeros
        V := Parent(Moverk[1]);
        zero := Zero(V);
        Moverk := Matrix([ Eltseq(Moverk[i]) : i in [1..NumberOfRows(Moverk)]
                        | Moverk[i] ne zero ]);
        return Moverk;
end intrinsic;


//////////////////////////////////////////////////////////////////////
//
//	Imposing a basepoint on a linear system
//
//////////////////////////////////////////////////////////////////////

forward impose1, impose2;

intrinsic ImposeBasepoint(X::Sch,L::LinearSys,p::Pt,m::RngIntElt) -> LinearSys
{}
    require Dimension(X) eq 2: "X must be a surface";
    require IsNonsingular(X,p): "p must be a nonsingular point of X";
    require Ambient(L) eq Ambient(X): "L and X must lie in the same ambient space";
    d := Degree(L);
    return L meet LinearSystem(Ambient(L),impose1(X,d,p,m));
end intrinsic;

intrinsic ImposeBasepoint(X::Sch,L::LinearSys,Z::Sch,m::RngIntElt) -> LinearSys
{Let X be a surface in projective 3-space P3, L a linear system on P3
and p a nonsingular point of X (or Z a subscheme of degree 2).
Return the sublinear system of L of forms that vanish at p (or Z) to
order at least m.}
    require Dimension(X) eq 2: "X must be a surface";
    require Ambient(L) eq Ambient(X): "L and X must lie in the same ambient space";
    deg := Degree(Z);
    require deg in {1,2}: "The proposed basepoint Z must have degree 1 or 2 over the base field";
    if deg eq 1 then
        p := Support(Z);
        return ImposeBasepoint(X,L,p,m);
    elif deg eq 2 then
        return L meet LinearSystem(Ambient(L),impose2(X,Degree(L),Z,m));
    end if;
end intrinsic;

//////////////////////////////////////////////////////////////////////

function impose1(X,d,p,m)
    conds_mx := VanishingConditions(X,p,d,m);
    P3 := Ambient(X);
    monos_on_P3 := MonomialsOfDegree(CoordinateRing(P3),d);
    N := #monos_on_P3;
    conds := Kernel(Transpose(conds_mx));
    polys := [ &+[ Eltseq(conds.i)[j] * monos_on_P3[j] : j in [1..N] ]
                        : i in [1..Dimension(conds)] ];
    return polys;
end function;

// This is the case where p is a point of degree 2 - Type = Sch.
function impose2(X,d,p,m)
    P3 := Ambient(X);
    k1 := Degree2SplittingField(p);
    P3_k1 := BaseChange(P3,k1);
    X_k1 := BaseChange(X,P3_k1);
    Z_k1 := BaseChange(p,P3_k1);
    p_k1 := Representative(Support(Z_k1));
    conds_mx := VanishingConditions(X_k1,p_k1,d,m);
      // These conditions are linear over k1. Resolve them into twice as
      // many conditions over the original field k (see e.g. below)
      // by splitting into 'real and imaginary' parts.
    entries_over_k1 := Eltseq(Transpose(conds_mx));
    entries_over_k := &cat[ Eltseq(k1!z) : z in entries_over_k1 ];
    conds_mx2 := Transpose(Matrix(2*NumberOfRows(conds_mx),entries_over_k));
    conds := Kernel(Transpose(conds_mx2));
      // OK: calculation is now back over the original field k.
    monos_on_P3 := MonomialsOfDegree(CoordinateRing(P3),d);
    N := #monos_on_P3;
    polys := [ &+[ Eltseq(conds.i)[j] * monos_on_P3[j] :
                j in [1..N] ] : i in [1..Dimension(conds)] ];
    return polys;
end function;

/*
sample to explain the "split real and imag parts of a matrix" line.
> k := Rationals();
> U:=PolynomialRing(k);
> k1<a> := ext<k|U.1^2+2>;
> Ma := Matrix(2,[1+2*a,3+4*a,5+6*a,7+8*a,9+10*a,11+12*a]);
> Mk := Transpose(Matrix(2*NumberOfRows(Ma),&cat[ Eltseq(k1!z) : z in Eltseq(Transpose(Ma))]));
> Ma;
[  2*a + 1   4*a + 3]
[  6*a + 5   8*a + 7]
[ 10*a + 9 12*a + 11]
> Mk;
[ 1  3]
[ 2  4]
[ 5  7]
[ 6  8]
[ 9 11]
[10 12]
*/


//////////////////////////////////////////////////////////////////////
//                      Geiser
//
// Input:
//              X       cubic surface in P^3=P(1,1,1,1) over a field
//              p in X  nonsingular rational point
//
// Output:
//              phi     quadratic map P3 - -> P3 that takes X->X.
//
// Notes:
// 1. Don't know (or test) whether X needs to be smooth and of
// Picard rank 1.
// 2. Don't know whether we've used char=0 or not.
//
// Convention: treat all points as lying in P3 when it makes sense.
//
//////////////////////////////////////////////////////////////////////

forward LineParametrisation, simplex_map;

intrinsic GeiserInvolution(X::Sch,p::Pt) -> MapSch
{Let X be a cubic surface in P3 and p a nonsingular point of X.
Return the Geiser involution X -> X centred in p. (X is not required
to be minimal - the usual requirement - but
in that case the Geiser involution may not be defined
and an error will be reported.)}
    // Check that X is cubic surface and p in X is nonsingular, etc.
    P3 := Ambient(X);
    require (Type(P3) eq Prj and Gradings(P3) eq [[1,1,1,1]]):
        "X must lie in (standard unweighted) projective 3-space";
    require (IsProjective(X) and Dimension(X) eq 2 and IsIrreducible(X)):
        "X must be an irreducible, projective surface";
    bool,p := IsCoercible(X,p);
    require bool: "p must lie on X";
    require IsNonsingular(X,p): "X must be nonsingular at p";
    p := P3 ! p;

    // If Tp \cap X splits as line+conic, then X isn't minimal:
    // an assert at the end will fail in this case.

    Tp := TangentSpace(X,p);
    Xp := Intersection(X,Tp);
        // Either Xp is irreducible over kbar, or
        // Xp splits over k, or p is an Eckardt point.
    /*
    if not IsIrreducible(Xp) then
        print "Plane section splits over k:",IrreducibleComponents(Xp);
        error if true, "X does not have Picard rank 1";
    end if;
    */

    k := BaseRing(P3);
    // we must work over an exact field (the rationals is checked separately)
    require (IsField(k) and (Type(k) eq FldRat or IsExact(k))):
        "X must be defined over an exact field";
    // Note: we need to make quadratic extensions of k, and magma
    // may impose additional conditions on us to be able to do that.

    // Get the nonlinear equations for the map from VanishingConditions
    // function and make a map G0 : P3 -> P3 using them.
    // Note that this won't necessarily take X -> X because we have
    // not taken care to choose the right basis of H^0(X,2A-3p).

    L := ImposeBasepoint(X,LinearSystem(P3,2),p,3);
    require #Sections(L) eq 4:
        "h^0(X,2A - 3p) is not 4; in particular, X is not minimal";
    G0 := map< P3 -> P3 | Sections(L) >;

    // Find the linear change of coordinates H : P3 -> P3 to arrange
    // that the Geiser G0*H (G0 first) really maps X -> X in the right way.
    // Method: compute >= 5 points that span P3 and whose images
    // under the final Geiser involution G we know. Then rig up the
    // linear map to match these correct images with those images under G0.

    // First find some good lines through p that have 3 distinct
    // intersection points.
    // The Geiser will interchange the two points that are not p
    // (although we may have to make a field extension to see them).
    Ip := Ideal(Cluster(p));
    B := Basis(Ip);
    Ptpairs := [];
    L1 := Scheme(P3,[B[1],B[2]]);
    L2 := Scheme(P3,[B[1],B[3]]);
    L3 := Scheme(P3,[B[2],B[3]]);
    lines := [ L1,L2,L3 ];
    Li := [];
    for L in lines do
        // Note: L cannot lie inside X (by the Tp condition at the start),
        // but we keep the next line for testing daft examples.
        if L subset X then
            continue L;
        end if;
        // intersect L with X and remove the known point p
        pair := Scheme(P3,ColonIdeal(Ideal(Intersection(X,L)),Ip));
        if (not(pair in Ptpairs)) and IsReduced(pair) then
            Append(~Li, L);
            Append(~Ptpairs, pair);
        end if;
    end for;
    easy_pairs := Ptpairs;

repeat
    // We need at least 5 pairs of points, so generate some
    // new painfully random ones if we haven't succeeded yet.
    Ptpairs := easy_pairs;      // we may repeat this step, so initialise.
    while #Ptpairs lt 5 do
        rands := [ Random([1..10]) : i in [1..6] ];
        L := Scheme(P3,[rands[1]*B[1]+rands[2]*B[2]+rands[3]*B[3],
                rands[4]*B[1]+rands[5]*B[2]+rands[5]*B[3]]);
        // We avoid cases where L lies in Tp(X) (so Eckardt points work e.g.)
        if L subset Tp then
            continue;
        end if;
        pair := Scheme(P3,ColonIdeal(Ideal(Intersection(X,L)),Ip));
        if (not(L in Li)) and IsReduced(pair) then
            Append(~Li, L);
            Append(~Ptpairs, pair);
        end if;
    end while;

    // Extend the field to see the 10 points.
    k1 := k;    // k1 is our working field extension
    for i in [1..5] do
        L := Li[i];
        psi,PP1<u,v> := LineParametrisation(L);
        Z := Ptpairs[i];
        pullback_of_Z := [ AlgebraMap(psi)(f) : f in Basis(Ideal(Z))];
        F := GroebnerBasis(pullback_of_Z)[1];
        // dehomogenise this quadratic poly (no extension is nec'y
        // if the poly is not quadratic).
        RR := PolynomialRing(k1);
        F1 := Evaluate(F,[RR.1,1]);
        if Degree(F1) eq 2 and IsIrreducible(F1) then
                // we must extend k1 by a root of F1
                k1 := ext< k1 | F1 >;
        end if;
    end for;
    five_pairs := [ SetToSequence(Support(Z,k1)) : Z in Ptpairs ];
    five_ps := [ Z[1] : Z in five_pairs ];
    five_qs := [ Z[2] : Z in five_pairs ];
    // The Geiser involution should take the five_ps to the five_qs, resp'y.
    // Tiny check: we can't allow any of the five_ps to be the original pt p;
    // it's fine for the five_qs though.
    for i in [1..5] do
        if five_ps[i] eq p then
            // switch this p and q.
            five_ps[i] := five_qs[i];
            five_qs[i] := p;
            require five_ps[i] ne p: "Unexpected error; please report (X,p): p lies triply on some line";
        end if;
    end for;

    // Test that the sets of points form a simplex; don't need to check qs too.
    M := Matrix([Eltseq(p):p in five_ps]);
until &and[ m ne 0 : m in Minors(M,4)];

    // We know that the Geiser involution must do:
    //      p1 <-> q1, p2 <-> q2, etc.
    // (maybe fewer than 5 of these would be enough... that would speed it up.)
    // (We don't use that it must usually also take _any_ point on Xp to p
    // because it's not true at Eckardt points, for instance.)
    // So build a linear H : P3 -> P3 (defined over the initial base
    // field) that takes G0(p1) -> q1 etc.
    // This should work, because we've got 5 lin indep points in P3.
    // But, of course, we work over the bigger field k1 at first.

    P3_ext := BaseChange(P3,k1);
    H1 := simplex_map(five_qs, P3_ext);
    H2 := simplex_map([G0(pp) : pp in five_ps], P3_ext);
    H := H2^(-1) * H1;

    // Now define the same map H over the smaller field k - call it H0.
    // (Point: often there's a factor in the equations from k1 which makes
    // H not def'd over k - but we only went to k1 to concoct H.)
    // Then we'll do G0 * H0 to get the answer.

    aaa := DefiningEquations(H);        // Note: aaa cannot be zero.
    for j in [1..4] do
        for i in [1..#aaa] do
            a := MonomialCoefficient(aaa[i],P3_ext.j);
            if a ne 0 then
                break j;
            end if;
        end for;
    end for;
    bbb := [ f/a : f in aaa];
    H0 := Automorphism(P3, [CoordinateRing(P3) ! f : f in bbb]);

    // finally, the Geiser trans is G = first G0 then H0:
    G := G0 * H0;
    assert five_qs eq G(five_ps);
    return G;   // use Expand(G) to multiply out the composition, if you like
end intrinsic;


//////////////////////////////////////////////////////////////////////
//                      Bertini
//
// Input:
//              X       cubic surface in P^3=P(1,1,1,1) over a field
//              Z in X  subscheme of length 2 (nonsingular points over kbar)
//
// Output:
//              phi     quadratic map P3 - -> P3 that takes X->X.
//
// Notes:
// 1. Don't know (or test) whether X needs to be smooth and of
// Picard rank 1.
// 2. Don't know whether we've used char=0 or not.
//
// Convention: treat all points as lying in P3 when it makes sense.
//
//////////////////////////////////////////////////////////////////////


forward remove_comp, param, bertini_data;

intrinsic BertiniInvolution(X::Sch,Z::Sch) -> MapAutSch
{Let X be a cubic surface in P3 and Z a finite subscheme of X of length 2.
Return the Bertini involution X -> X centred in Z.}
    k := BaseRing(X);
    // We must work over an exact field (the rationals is checked separately)
    require (IsField(k) and (Type(k) eq FldRat or IsExact(k))):
        "X must be defined over an exact field";
    // Note: we need to make quadratic extensions of k, and magma
    // may impose additional conditions on us to be able to do that.

    P3 := Ambient(X);
    require (Type(P3) eq Prj and Gradings(P3) eq [[1,1,1,1]]):
        "X must lie in (standard unweighted) projective 3-space";
    // Check that X is cubic surface and p in X is nonsingular, etc.
    require (IsProjective(X) and Dimension(X) eq 2 and IsIrreducible(X)):
        "X must be an irreducible, projective surface";
    require TotalDegree(DefiningEquation(X)) eq 3: "X must be a cubic surface";
    bool := Z subset X;
    require bool: "Z must lie on X";
    // Nonsingularity of X at Z tested in ImposeBasepoint algorithm.

    L0 := ImposeBasepoint(X,LinearSystem(P3,5),Z,6);
    L := Complement(L0,X);
    require #Sections(L) eq 4:
        "h^0(X,5A - 6Z) is not 4; in particular, X is not minimal";
    B0 := map< P3 -> P3 | Sections(L) >;

    // Note that B0 won't necessarily take X -> X because we have
    // not taken care to choose the right basis of H^0(X,5A-6p).
    // We change basis as follows. Take 5 points of X - these are called Ps.
    // This determines two other sets of five points, Qs and Rs:
    // Qs = five images under B0; Rs = what we know the five images
    // must be under the Bertini involution.
    // So all we have to do is use a linear automorphism to shift Qs to Rs.

repeat
    // First get those points, Ps and Rs:
    invol_pairs := bertini_data(X,Z); // = < [P1,-P1],[P2,-P2],...>

    // Second, map the Qs to the Rs.
    coeffs := [];
    for k in [1..5] do
        P := invol_pairs[k][1];
        R := invol_pairs[k][2];
        Q := B0(P);

        RR<[u]> := PolynomialRing(Parent(Q[1]),16);
        MQ := Matrix(4,[RR|a:a in Eltseq(Q)]);
        b1,b2,b3,b4 := Explode([RR|a:a in Eltseq(R)]);
        // We will try to solve Q * (4x4 matrix in u1..u16) = scalar mult of R.
        // This is the same as Q * (4x4...) * ker R = (0,0,0).
        // So we compute ker R = 4x3 matrix of rank 3 with kernel-like columns.
        // In the end, a trick with trace gets this condition back to k.
        if b1 ne 0 then
            MR := Transpose(Matrix(4,[RR| b2,-b1,0,0, b3,0,-b1,0, b4,0,0,-b1]));
        elif b2 ne 0 then
            MR := Transpose(Matrix(4,[RR| -b2,b1,0,0, 0,b3,-b2,0, 0,b4,0,-b2]));
        elif b3 ne 0 then
            MR := Transpose(Matrix(4,[RR| -b3,0,b1,0, 0,-b3,b2,0, 0,0,b4,-b3]));
        else // b4 must be nonzero
            MR := Transpose(Matrix(4,[RR| -b4,0,0,b1, 0,-b4,0,b2, 0,0,-b4,b3]));
        end if;
        // to believe it's the right way round, compare with:
        // Transpose(Matrix(4,[RR| 1,-1,0,0, 2,0,-2,0, 3,0,0,-3]));
        Mu := Matrix(4,u);
        N := MQ * Mu * MR;
        coeffs cat:=
            [ [Trace(MonomialCoefficient(N[1,j],u))
                : u in [RR.i:i in [1..Rank(RR)]]]
                    : j in [1..3] ];
    end for;
    KK := Kernel(Transpose(Matrix(coeffs)));
until Dimension(KK) eq 1;
    // this dim = 1 condition could fail if the five first elts of invol_pairs
    // fail to span a simplex in P3 - rather than checking that, we simply
    // run the whole routine again (a bit daft).
    phi := Automorphism(P3,Matrix(4,Eltseq(Basis(KK)[1])));
    return B0 * phi;
end intrinsic;


function bertini_data(X,Z)
    P3 := Ambient(Z);
    // Step 1: find the line L through Z.
    L1 := LinearSystem(P3,1);
    lineqs := Sections(LinearSystem(L1,Z)); // the linear eqns in I_Z.
    L := Scheme(P3,lineqs);
    // Step 2: find the 3rd point of intersection q---must be k-rational.
    q := Representative(Support(remove_comp(Z,Intersection(X,L))));
    x_pts := <>;
    E_so_far := [];
    for i in [1..5] do
        // Step 3: generate a random plane Pi containing L; check
        // that E=X int Pi is nonsingular.
        repeat
            F := &+[ Random([1..100]) * l : l in lineqs ];
            E := Scheme(X,F);
        until E notin E_so_far and IsNonsingular(E);
        Append(~E_so_far,E);
        Pi := Scheme(P3,F);
        para := param(Pi);
        // Step 4: realise E as an elliptic curve magma type.
        Epar := Curve(E @@ para);
        Zpar := Z @@ para;
        assert Zpar subset Epar;
        qpar := Representative(Support(q@@para));
        assert qpar in Epar;
        Epec,fec := EllipticCurve(Epar,qpar);
        Zpar := Intersection(Zpar,Epar); // nuts, but asserts Zpar subset Epar.
        // Zec := Zpar @@ (fec^-1); //why doesn't this work? do by hand instead:
        fec_poly := AlgebraMap(fec^-1);
        Zec := Scheme(Epec, [ fec_poly(f) : f in Basis(Ideal(Zpar)) ]);
        // Step 5: use E,q as EllCrv to find point pairs (P,-P) on E
        // Note: could use TwoTorsion(E) just in case it has any / k.
        P2ec := Ambient(Epec); U := P2ec.1; W := P2ec.3;
        repeat
            line := Scheme(P2ec,&+[ Random([1..100]) * var : var in [U,W] ]);
            twopoints :=
                remove_comp(Cluster(P2ec![0,1,0]),Intersection(Epec,line));
            k2 := Degree2SplittingField(twopoints);
            supp := Support(twopoints,k2);
        until not (supp subset Zec);    // ensuring that supp is not Z.
        // Step 6: return the pair (P,-P) to X.
        // maps so far: E <-para- Epar -fec-> Epec
        Append(~x_pts, SetToSequence(((fec^-1) * para) (supp)));
    end for;
    // Remark: this function should really return 5 simplex pts; we should
    // check that condition here and continue searching if it fails.
    // The calling code is smart and will ask again if some necessary
    // condition fails - so it's not a math disaster, but it's not very
    // efficient if this gets called twice.
    return x_pts;
end function;


//////////////////////////////////////////////////////////////////////
//
//			Auxilliary functions
//
//////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////
// Input: zero-dimensional scheme Z of degree 2 in P^n.
// Output: a field extension over which all the points of Z are defined.
//////////////////////////////////////////////////////////////////////

intrinsic Degree2SplittingField(Z::Sch) -> Rng
{Let Z be a zero-dimensional scheme of degree 2 in P^n over a field k.
Return an extension of k (of degree 1 or 2) over which Z splits geometrically.}
    require Dimension(Z) eq 0: "Z must be a zero-dimensional scheme";
    require Degree(Z) eq 2: "Z must be a scheme of degree 2";
    Pn := Ambient(Z);
    i := 0;
    repeat
        // run through affine patches looking for Z
        i +:= 1;
        ZA := AffinePatch(Z,i);
        bad := IsEmpty(ZA);
    until not bad;
    // in the affine patch ZA, Z must be defined by linear eqns and 1 quadratic
    RA := CoordinateRing(Ambient(ZA));
    IZA := Ideal(ZA);
    coord_polys := [ Basis(EliminationIdeal(IZA,{j}))[1] :
                                        j in [1..Dimension(Pn)]];
    found_irred := false;
    for f in coord_polys do
        if TotalDegree(f) ge 2 then
            if IsIrreducible(f) then
                found_irred := true;
                goodf := f;
                break f;
            end if;
        end if;
    end for;
    if not found_irred then
            // The two points are already defined over the basefield.
            return BaseField(Z);
    else
        k := BaseField(Z);
        U := PolynomialRing(k);
                // stupid way to turn f into univariate poly:
        fU := Evaluate(goodf, [U.1:i in [1..Dimension(Pn)]]);
        k1 := ext< k | fU >;
        return k1;
    end if;
end intrinsic;


//////////////////////////////////////////////////////////////////////
// Input: X in P3 a nonsingular surface, f : X -> P^1.
// Output: true iff there is a maxl centre, together with the
//      first basepoint (of deg 1 or 2---checking deg 1 first)
//      on which f has mult > deg(f).
//////////////////////////////////////////////////////////////////////

forward test_basepoint;
intrinsic HasMaximalCentre(f::MapSch,X::Sch) -> BoolElt,Sch,RngIntElt
{Let X be a cubic surface in P3 and f : P3 -> P^1 a map that restricts to
an  elliptic pencil on X, i.e. a map defined by two homogeneous functions
(f:g) with genus 1 fibres.
Return true iff f has a basepoint on X of multiplicity at least 2, in which
case also return the basepoint as a scheme and its multiplicity.}
    f1 := DefiningEquations(f)[1];
    f2 := DefiningEquations(f)[2];
    degf := TotalDegree(f1);
    P3 := Ambient(X);
    monos := MonomialsOfDegree(CoordinateRing(P3),degf);
    coeffs := Transpose(Matrix(#monos,
         [ MonomialCoefficient(f1,m) : m in monos ] cat
            [ MonomialCoefficient(f2,m) : m in monos ]));

    // Step 1: retrieve set-theoretic base locus and split it into irreducibles.
    Ef := DefiningEquations(f);
    BL := ReducedSubscheme(Scheme(X,Ef));
    comp := IrreducibleComponents(BL);

    // Step 2: compute multiplicities of basepoints of deg <= 2 only.
    points := [];
    curves := [];
    for Z in comp do
        if Dimension(Z) ne 0 then
            Append(~curves,Z);
            continue Z;
        else
            Append(~points,Z);
        end if;
        bool,deg := test_basepoint(X,Z,degf,coeffs);
        if bool then
            return bool,Z,deg;
        end if;
    end for;

    // Step 3: haven't found anything yet; look for genuine base locus
    // embedded in the codim 1 part.
    // First remove any point components we've already found.
    I := ideal< CoordinateRing(P3) | DefiningEquation(X), Ef >;
    for P in points do
        I := Saturation(I,Ideal(P));
    end for;
    // Second remove just the right amount of curve components,
    // but not so much that you remove embedded points.
    Ired := Radical(I);
    for C in curves do
        IC := Ideal(C);
        I1 := I;
        repeat
            I1 := ColonIdeal(I1,IC);
        until not (I1 subset Ired);
        I := I1;
    end for;
    // Third, split the remaining ideal.
    Irad := Radical(I);
    if Irad eq Generic(I) then
        comp := [];
    else
        comp := IrreducibleComponents(Scheme(P3,Radical(I)));
    end if;
    for Z in comp do
        if Dimension(Z) ne 0 then
            Append(~curves,Z);
            continue Z;
        else
            Append(~points,Z);
        end if;
        degZ := Degree(Z);
        bool,deg := test_basepoint(X,Z,degf,coeffs);
        if bool then
            return bool,Z,deg;
        end if;
    end for;

    return false,_,_;
end intrinsic;

// return true (and the degree of Z) iff the linear system (given implicitly by
// the sequence 'coeffs') of polys of degree d vanishes to order >= d+1 at p.
function test_basepoint(X,Z,d,coeffs)
    degZ := Degree(Z);
    if degZ eq 1 then
        pt := Representative(Support(Z));
        conds_mx := VanishingConditions(X,pt,d,d+1);
        if IsZero(conds_mx*coeffs) then
            // this means that the lin sys vanishes too much at Z,
            // so Z is a centre that could be untwisted.
            return true,1;
        end if;
    elif degZ le 2 then
        // Note: the condition is the same computed at either point.
        k2 := Degree2SplittingField(Z);
        P := Ambient(X);
        P2 := BaseChange(P,k2);
        X2 := BaseChange(X,P2);
        Z2 := BaseChange(Z,P2);
        pt := Representative(Support(Z2));
        conds_mx := VanishingConditions(X2,pt,d,d+1);
        coeffs := ChangeRing(coeffs,k2);
        if IsZero(conds_mx*coeffs) then
            return true,2;
        end if;
    end if;
    return false,_;
end function;


//////////////////////////////////////////////////
// Input: a sequence of 5 points in a P3
// Output: an automorphism of P3 taking the standard simplex
// to the given 5 points.

function simplex_map(Q,P3)
        M1 := Matrix(4, &cat [Eltseq(p) : p in Q[1..4]]);
        lambdas := Solution(M1,Matrix(4,Eltseq(Q[5])));
        for i in [1..4] do
                M1[i] *:= lambdas[1,i];
        end for;
        M1 := ChangeRing(M1,BaseRing(P3));
        phi := Automorphism(P3,M1);
        return phi;
end function;
//////////////////////////////////////////////////


//////////////////////////////////////////////////
// Input: L a line in P3 given by two linear forms f,g
// Output: map P1 -> P3 parametrising L.

function line_param(P1,P3,f,g)
    k := BaseField(P3);
    fcoeffs := [k| MonomialCoefficient(f,P3.i) : i in [1..4]];
    gcoeffs := [k| MonomialCoefficient(g,P3.i) : i in [1..4]];
    M := Matrix([fcoeffs,gcoeffs]);
    K := Kernel(Transpose(M));
    u := P1.1;
    v := P1.2;
    map_eqns := [ (K.1)[i]*u + (K.2)[i]*v : i in [1..4] ];
    return map< P1 -> P3 | map_eqns >;
end function;

function LineParametrisation(L)
    P3 := Ambient(L);
    k := BaseField(P3);
    P1 := ProjectiveSpace(k,1);
    f,g := Explode(GroebnerBasis(L));
    return line_param(P1,P3,f,g),P1;
end function;

function plane_param(P2,P3,f)
    k := BaseField(P3);
    fcoeffs := [k| MonomialCoefficient(f,P3.i) : i in [1..4]];
    M := Matrix([fcoeffs]);
    K := Kernel(Transpose(M));
    u := P2.1;
    v := P2.2;
    w := P2.3;
    map_eqns := [ (K.1)[i]*u + (K.2)[i]*v + (K.3)[i]*w : i in [1..4] ];
    return map< P2 -> P3 | map_eqns >;
end function;

function PlaneParametrisation(L)
    P3 := Ambient(L);
    k := BaseField(P3);
    P2 := ProjectiveSpace(k,2);
    f := Explode(GroebnerBasis(L));
    return plane_param(P2,P3,f),P2;
end function;
//////////////////////////////////////////////////


//////////////////////////////////////////////////
// X a cubic surface, p a k-rational point on it.
// Return true iff p is an Eckardt point, i.e iff
// TpX intersects X in 3 lines.

function is_eckardt_point(X,p)
    Tp := TangentSpace(X,p);
    Xp := Intersection(X,Tp);
    phi := PlaneParametrisation(Tp);
    C := Curve(Xp @@ phi);
    pp := Representative(Support(p @@ phi));
    return TangentCone(C,pp) eq C;
end function;
//////////////////////////////////////////////////


//////////////////////////////////////////////////

function trace(p,k)
    return Scheme(p) ! [ Trace(a,k) : a in Eltseq(p) ];
end function;

// imagine a P2 in a P3: return map< P2 - P3 > with that image.
function param(plane)
    Pn := Ambient(plane);
    R := CoordinateRing(Pn);
    F := DefiningEquations(plane);
    // should check codim?
    f := F[1];
    M := Matrix(1,[ MonomialCoefficient(f,R.i) : i in [1..Length(Pn)] ]);
    BM := Basis(Kernel(M));
    P2 := ProjectiveSpace(BaseField(Pn),2);
    S := CoordinateRing(P2);
    map_polys := [ &+[ BM[j][i] * S.j : j in [1..#BM] ] :
                   i in [1..Length(Pn)] ];
    return map< P2 -> Pn | map_polys >;
end function;

// A,B two schemes; assume A subset B; return scheme of colonideal(B,A).
function remove_comp(A,B)
    IA := Ideal(A);
    IB := Ideal(B);
    assert IB subset IA;
    return Scheme(Ambient(B),ColonIdeal(IB,IA));
end function;