/****************************************************
 * Reduction of ternary forms
 *
 * Michael Stoll, started 2002-02-21
 ****************************************************/

/****************************************************
 We associate to a nonsingular ternary form C(x,y,z)
 a positive definite quadratic form Q(x,y,z).
 
 This is done as follows. When C is a nonsingular cubic:
 Let H(x,y,z) = det(matrix of second partial derivatives of C)
 be the Hessian of C. Then the intersection of C and H
 consists of 9 distinct points, the flexes of C.
 Three of them are real, the others come in three complex
 conjugate pairs. There are twelve lines each containing
 three of teh flexes, coming in four triples of lines
 that do not meet in a flex. One of these triples has all
 three lines real, call them L_11, L_12, L_13. Another one
 has one line real, call it L_21, and two complex conjugate
 lines, call them L_22 and L_23. Then our form Q spans
 the one-dimensional intersection of the spaces spanned by
 L_11^2, L_12^2 and L_13^2, and by L_21^2 and L_22 L_23,
 respectively.
 
 This Q is positive definite and (as a Hermitian form) 
 a covariant of C (up to scaling); in fact the unique
 (real) Hermitian covariant of C.
 
 If C is a higher degree form, we still compute the flexes
 (3 d (d-2) of them, where d is the degree of C). Now we
 use them to find a covariant positive definite Hermitian
 (or real quadratic) form Q by some optimisation procedure.
 
 Q defines a lattice which we can Minkowsi- or LLL-reduce.
 Applying the corresponding transformations to C, we
 obtain a reduced curve.
 
 I.e., C is (LLL-)reduced if Q is reduced if the lattice
 defined by Q is Minkowski-reduced (LLL-reduced).
 ****************************************************/

declare verbose TernaryReduce, 3;
declare verbose Minimise, 3;

intrinsic PolMat(pol::RngMPolElt, mat::Mtrx) -> RngMPolElt
{pol^mat; workaround for bug in Magma.}
  P := Parent(pol);
  n := Rank(P);
  require n eq NumberOfRows(mat) and n eq NumberOfColumns(mat):
          "mat must be a square matrix of order equal to the number of variables in pol.";
  return Evaluate(pol, [&+[mat[i,j]*P.j : j in [1..n]] : i in [1..n]]);
end intrinsic;

// Generate random matrices in SL(n, Z)
intrinsic RandomSLnZ(n::RngIntElt, k::RngIntElt, l::RngIntElt) -> AlgMatElt
{A random element in SL(n,Z). k and l have an influence on the size
 of the entries.}
  eltmat := func<i,j,e | Matrix([[ii eq jj select 1
                                   else ii eq i and jj eq j select e
                                   else 0
                                   : jj in [1..n]] : ii in [1..n]])>;
  return &*[eltmat(Random(1,n), Random(1,n), Random(-k,k)) : j in [1..l]];
end intrinsic;

// Generate random forms
intrinsic RandomForm(P::RngMPol, deg::RngIntElt, k::RngIntElt) -> RngMPolElt
{Generates a random form of degree deg in P. k bounds the coefficients.}
  return &+[Random(-k,k)*m : m in MonomialsOfDegree(P, deg)];
end intrinsic;

intrinsic ComplexVariety(C::RngMPolElt, H::RngMPolElt) -> SeqEnum
{Computes the complex points of the zero-dimensional projective variety 
 given by the two ternary forms C and H over Q.}
  P := Parent(C);
  Q := CoefficientRing(P);
  require Type(Q) eq FldRat: "Base field must be the rationals.";
  require P eq Parent(H): "Both forms must have the same parent.";
  require Rank(P) eq 3: "Forms must be in three variables.";
  dC := TotalDegree(C);
  dH := TotalDegree(H);
  require forall{m : m in Monomials(C) | TotalDegree(m) eq dC}
           and forall{m : m in Monomials(H) | TotalDegree(m) eq dH}:
          "Polynomials must be forms (i.e., homogeneous).";
  require TotalDegree(GCD(C, H)) eq 0: "Curves have a component in common.";
  PC := PolynomialRing(ComplexField());
  Proj := ProjectiveSpace(Q, 2);
  PP := CoordinateRing(Proj);
  CC := Curve(Proj, PP!C);
  CH := Curve(Proj, PP!H);
  if Evaluate(C, [1,0,0]) eq 0 and Evaluate(H, [1,0,0]) eq 0 then
    pt := Proj![1,0,0];
    points := [<[1,0,0], IntersectionNumber(CC, CH, pt)>];
  else
    points := [];
  end if;
  PQ := PolynomialRing(Q);
  res := GCD(Evaluate(C, [PQ.1,1,0]), Evaluate(H, [PQ.1,1,0]));
  fact := Factorisation(res);
  fact := [<f[1], IntersectionNumber(BaseChange(CC,PrK), BaseChange(CH,PrK), 
                                     PrK![r,1,0])>
             where PrK := BaseChange(Proj, K)
             where r := Evaluate(res, K.1) eq 0 select K.1
                         else Roots(res, K)[1,1]
             where K := ext<Q | f[1]>
            : f in fact];
  points cat:= &cat[[<[r[1], 1, 0], f[2]> : r in Roots(f[1], ComplexField())]
                      : f in fact];
  res := Evaluate(Resultant(Evaluate(C, [P.1,P.2,1]), Evaluate(H, [P.1,P.2,1]),
                            P.2),
                  [PQ.1,0,0]);
  fact := Factorisation(res);
  for f in fact do
    f1 := f[1];
    if Degree(f1) eq 1 then
      r := -Coefficient(f1,0);
      res1 := GCD(Evaluate(C, [r, PQ.1, 1]), Evaluate(H, [r, PQ.1, 1]));
      fact1 := SquarefreeFactorization(res1);
      for ff in fact1 do
        points cat:= [<[r, rs[1], 1], ff[2]>
                        : rs in Roots(ff[1], ComplexField())];
      end for;
    elif Degree(f1) le 6 or f[2] ge 2 then
      // Do computations exactly as far as possible 
      // to get muliplicities right.
      K := ext<Q | f1>;
      if Evaluate(f1, K.1) eq 0 then
        r := K.1;
      else
        r := Roots(f1, K)[1,1];
      end if;
      PK := PolynomialRing(K);
      res1 := GCD(Evaluate(C, [r, PK.1, 1]), Evaluate(H, [r, PK.1, 1]));
      fact1 := SquarefreeFactorization(res1);
      for rp in Roots(f1, ComplexField()) do
        h := hom<K -> ComplexField() | rp[1]>;
        for ff in fact1 do
          points cat:= [<[rp[1], r[1], 1], ff[2]>
                          : r in Roots(PC![h(c) : c in Coefficients(ff[1])])];
        end for;
      end for;
    else
      // Try to get by without working in high degree number fields.
      // Result may be less accurate, though.
      _, prec := HasAttribute(FldPr, "Precision");
      oldprec := prec;
      repeat
        ptsloc := [];
        flag := true;
        for r in [r[1] : r in Roots(f1, ComplexField())] do
          pol1 := Evaluate(C, [r, PC.1, 1]);
          pol2 := Evaluate(H, [r, PC.1, 1]);
          // We know that pol1 and pol2 have a unique common root.
          // Do the euclidean algorithm until we are left with a linear factor.
          pol1 *:= 1/LeadingCoefficient(pol1);
          pol2 *:= 1/LeadingCoefficient(pol2);
          p1 := pol1; p2 := pol2;
          while Degree(pol2) gt 1 do
            temp := pol2 mod pol1;
            if temp eq 0 then break; end if;
            pol2 := pol1;
            if Degree(pol2) gt 1 then
              pol1 := 1/LeadingCoefficient(temp)*temp;
            end if;
          end while;
          r2 := -Coefficient(pol2, Degree(pol2)-1)/Degree(pol2);
          if Norm(Evaluate(p1, r2)) ge 1.0e-6
              or Norm(Evaluate(p2, r2)) ge 1.0e-6
          then
            prec +:= 50;
            AssertAttribute(FldPr, "Precision", prec);
            vprintf TernaryReduce, 2:
                    "ComplexVariety: Increasing precision to %o\n", prec;
            flag := false;
            break;
          end if;
          Append(~ptsloc, <[ComplexField()|r, r2, 1], 1>);
        end for;
      until flag;
      AssertAttribute(FldPr, "Precision", oldprec);
      points cat:= ptsloc;
    end if;
  end for;
  return points;
end intrinsic;

// Find real and complex flexes of a ternary form over Q
intrinsic RealAndComplexFlexes(C::RngMPolElt) -> SeqEnum, SeqEnum
{Finds the real and complex flexes of a ternary form over Q.
 First value is the sequence of real flexes, second value is a sequence
 containing one complex flex out of each conjugate pair.}
  P := Parent(C);
  deg := TotalDegree(C);
  require Rank(P) eq 3
          and forall{m : m in Monomials(C) | TotalDegree(m) eq deg}:
          "C must be a ternary form.";
  vprintf TernaryReduce, 1: "Flexes: C = %o\n", C;
  // compute Hessian
  H := Determinant(Matrix([[Derivative(Derivative(C,i),j) : j in [1..3]]
                            : i in [1..3]]));
  vprintf TernaryReduce, 2: "  Hessian H = %o\n", H;
  flexes := &cat[[v[1] : i in [1..v[2]]] : v in ComplexVariety(C, H)];
  require #flexes eq 3*deg*(deg-2):
          "Did not find the correct number of flexes.";
  flr := [fl : fl in flexes | IsReal(fl[1]) and IsReal(fl[2])];
  flc0 := [fl : fl in flexes | not IsReal(fl[1]) or not IsReal(fl[2])];
  // find representatives of conjugate pairs
  flc := [];
  for i := 1 to ExactQuotient(#flc0, 2) do
    new := flc0[#flc0];
    cnew := [ComplexConjugate(new[i]) : i in [1..3]];
    Append(~flc, new);
    min, pos := Minimum([&+[Norm(cnew[i]-flc0[j][i]) : i in [1..3]]
                          : j in [1..#flc0-1]]);
    assert min le 1.0e-5*&+[Norm(new[i]) : i in [1..3]];
    // assert pos ne #flc0;
    Remove(~flc0, pos);
    Prune(~flc0);
  end for;
  vprintf TernaryReduce, 3: "Flexes: %o\n", flexes;
  return flr, flc;
end intrinsic;

// Write my own Kernel function, since I got a ``division by zero in gdiv''
function MyKernel(mat)
  // mat = [[a,b], [c,d], [e,f]]
  // find linear combination of rows that vanishes
  // so operate on columns
  V := KSpace(CoefficientRing(mat), 3);
  // find largest entry in matrix
  max, pos := Max([Abs(c) : c in Eltseq(mat)]);
  if IsEven(pos) then // 2nd col
    SwapColumns(~mat, 1, 2);
    pos -:= 1;
  end if;
  r := ExactQuotient(pos+1, 2); // row number
  mat *:= 1/mat[r,1];
  AddColumn(~mat, -mat[r,2], 1, 2); // rth col is (1 0)
  // find largest entry in 2nd col
  max, pos := Max([Abs(mat[i,2]) : i in [1..3]]);
  if pos eq r or max eq 0 then
    return Eltseq([V.i - mat[i,1]*V.r : i in [1..3] | i ne r][1]);
  end if;
  MultiplyColumn(~mat, 1/mat[pos,2], 2);
  AddColumn(~mat, -mat[pos,1], 2, 1); // pos'th col is (0 1)
  t := 6-r-pos; // the third element of {1,2,3} = {r,pos,t}
  return Eltseq(V.t - mat[t,1]*V.r - mat[t,2]*V.pos);
end function;

// Given the flexes of a ternary cubic, find the relevant lines
function LinesFromFlexes(freal, fcomplex)
  // Given sequences containing the three real flexes and three non-conjugate 
  // complex flexes of a ternary cubic, respectively, finds the
  // line through the 3 real flexes, its complementary factor, and
  // the three real lines through a real and two complex flexes.
  P3 := PolynomialRing(RealField(), 3);
  normalise := function(l)
                 _, pos := Maximum([Norm(a) : a in l]);
                 norm := Sqrt(&+[Norm(a) : a in l])
                           * l[pos]/Modulus(l[pos]);
                 return [a/norm : a in l];
               end function;
  line := function(pt1, pt2)
            l := MyKernel(Matrix([[pt1[i], pt2[i]] : i in [1..3]]));
            return normalise(l);
          end function;
  lines := [line(freal[1],freal[2]),
            line(freal[2],freal[3]),
            line(freal[3],freal[1])];
  liner := normalise([&+[l[i] : l in lines] : i in [1..3]]);
  vprintf TernaryReduce, 3: "  Real line: %o\n", liner;
  linec := [[Real(a) : a in line(fl, [ComplexConjugate(c) : c in fl])]
             : fl in fcomplex];
  vprintf TernaryReduce, 3: "  Real triple of lines:\n%o\n", linec;
  linetest1 := line(fcomplex[1], fcomplex[2]);
  linetest2 := line(fcomplex[1], [ComplexConjugate(c) : c in fcomplex[2]]);
  vprintf TernaryReduce, 3:
          "  Possible complex lines:\n  %o\n  %o\n", linetest1, linetest2;
  m1 := Min(Norm(&+[linetest1[i]*fcomplex[3,i] : i in [1..3]]),
            Norm(&+[linetest1[i]*ComplexConjugate(fcomplex[3,i])
                     : i in [1..3]]));
  m2 := Min(Norm(&+[linetest2[i]*fcomplex[3,i] : i in [1..3]]),
            Norm(&+[linetest2[i]*ComplexConjugate(fcomplex[3,i])
                     : i in [1..3]]));
  error if Min(m1, m2) ge 1.0e-10, "Lines do not fit!";
  vprintf TernaryReduce, 3: 
          "  Values of two possible lines at third point: %o, %o\n", m1, m2;
  linet := m1 lt m2 select linetest1 else linetest2;
  compl := Norm(linet[1])*P3.1^2 + Norm(linet[2])*P3.2^2 + Norm(linet[3])*P3.3^2
            + 2*Real(linet[1]*ComplexConjugate(linet[2]))*P3.1*P3.2
            + 2*Real(linet[1]*ComplexConjugate(linet[3]))*P3.1*P3.3
            + 2*Real(linet[2]*ComplexConjugate(linet[3]))*P3.2*P3.3;
  vprintf TernaryReduce, 3: "  Complementary factor:\n%o\n", compl;
  return [&+[liner[i]*P3.i : i in [1..3]], compl]
           cat [&+[l[i]*P3.i : i in [1..3]] : l in linec];
end function;

function orthonorm(basis)
  for i := 1 to #basis do
    s := [&+[basis[j,k]*basis[i,k] : k in [1..#basis[i]]] : j in [1..i-1]];
    if i gt 1 then
      basis[i] := [basis[i,k] - &+[s[j]*basis[j,k] : j in [1..#s]]
                    : k in [1..#basis[i]]];
    end if;
    n := Sqrt(&+[a^2 : a in basis[i]]);
    basis[i] := [a/n : a in basis[i]];
  end for;
  return basis;
end function;

// Now find the form from the lines
function CovariantFromLines(lines)
  // Finds the covariant positive definite Hermitian form from the
  // given lines through the flexes of a ternary cubic.
  P := Universe(lines);
  // look for element in space <lines[1]^2, lines[2]> that is closest to
  // <lines[3]^2, lines[4]^2, lines[5]^2>
  mons := MonomialsOfDegree(P, 2);
  w1 := [ MonomialCoefficient(lines[1]^2, mon) : mon in mons ];
  w2 := [ MonomialCoefficient(lines[2], mon) : mon in mons ];
  // orthonormalise
  vprintf TernaryReduce, 3: "  First basis:\n  %o\n  %o\n", w1, w2;
  w := orthonorm([w1, w2]);
  vprintf TernaryReduce, 3: "  orthonormalised:\n  %o\n  %o\n", w[1], w[2];
  // basis of other space
  v1 := [ MonomialCoefficient(lines[3]^2, mon) : mon in mons ];
  v2 := [ MonomialCoefficient(lines[4]^2, mon) : mon in mons ];
  v3 := [ MonomialCoefficient(lines[5]^2, mon) : mon in mons ];
  // orthonormalise
  vprintf TernaryReduce, 3: "  Second basis:\n  %o\n  %o\n  %o\n", v1, v2, v3;
  v := orthonorm([v1, v2, v3]);
  vprintf TernaryReduce, 3: "  orthonormalised:\n  %o\n  %o\n  %o\n",
                            v[1], v[2], v[3];
  // look at cos(phi)*w[1] + sin(phi)*w[2] such that distance to <v>
  // is minimal. Leads to
  scmat := [[&+[wi[k]*vj[k] : k in [1..#mons]] : wi in w] : vj in v];
  vprintf TernaryReduce, 3: "  Sclar products matrix:\n%o\n", scmat;
  phi := 1/2*Arctan2(&+[a[1]*a[2] : a in scmat],
                     1/2*&+[a[1]^2 - a[2]^2 : a in scmat]);
  cphi := Cos(phi); sphi := Sin(phi);
  qs1 := [cphi*w[1,k] + sphi*w[2,k] : k in [1..#mons]];
  qs2 := [-sphi*w[1,k] + cphi*w[2,k] : k in [1..#mons]];
  ls1 := [cphi*a[1] + sphi*a[2] : a in scmat];
  ls2 := [-sphi*a[1] + cphi*a[2] : a in scmat];
  dist1 := 1 + &+[l^2 : l in ls1]
             - 2*&+[ls1[i]*(cphi*scmat[i,1] + sphi*scmat[i,2])
                     : i in [1..#ls1]];
  dist2 := 1 + &+[l^2 : l in ls2]
             - 2*&+[ls2[i]*(-sphi*scmat[i,1] + cphi*scmat[i,2])
                     : i in [1..#ls2]];
  if dist1 lt dist2 then
    qs := qs1; ls := ls1; dist := dist1;
  else
    qs := qs2; ls := ls2; dist := dist2;
  end if;
  vprintf TernaryReduce, 3: "  Nearest vecor in first space:\n  %o\n", qs;
  vv := [&+[ls[i]*v[i,k] : i in [1..#ls]] : k in [1..#mons]];
  vprintf TernaryReduce, 3: "  Nearest vector in second space:\n  %o\n", vv;
  Q := &+[qs[i]*mons[i] : i in [1..#mons]];
  if MonomialCoefficient(Q, Parent(Q).1^2) lt 0 then Q := -Q; end if;
  Qmat := Matrix([[(i eq j select 1 else 1/2)*MonomialCoefficient(Q, P.i*P.j)
                     : j in [1..3]] : i in [1..3]]);
  error if not IsPositiveDefinite(Qmat),
          "Covariant matrix is not positive definite -- try with higher precision.";
  Qmat *:= Determinant(Qmat)^(-1/3);
  return Qmat, dist;
end function;

intrinsic Covariant(C::RngMPolElt : Al := "Default", Eps := 1.0e-20)
  -> RngMPolElt
{Computes the positive definite Hermitian covariant of a ternary cubic
 as a real quadratic form.}
  require Rank(Parent(C)) eq 3: "C must be a ternary form.";
  deg := TotalDegree(C);
  require forall{m : m in Monomials(C) | TotalDegree(m) eq deg}:
          "Polynomial must be homogeneous.";
  require deg ge 3: "Degree must be at least 3.";
  require Al in {"Default", "SteepestDescent", "Newton"}:
          "Al must be one of \"Default\", \"SteepestDescent\" and \"Newton\"";
  if deg eq 3 and Al eq "Default" then
    flr, flc := RealAndComplexFlexes(C);
    return CovariantFromLines(LinesFromFlexes(flr, flc));
  else
    flr, flc := RealAndComplexFlexes(C);
    return Covariant(flr, flc : Al := Al, Eps := Eps);
  end if;
end intrinsic;

///////////////////////////////////////////////////////////////////////
// Ad-hoc reduction functions
///////////////////////////////////////////////////////////////////////

function siz(pol)
  return Max([Abs(c) : c in Coefficients(pol)]);
end function;

function redstep(cubic)
  P := Parent(cubic);
  z := [P.1, P.2, P.3];
  trans := IdentityMatrix(Integers(), 3);
  for i := 1 to 3 do
    mc := MonomialCoefficient(cubic, z[i]^3);
    if mc ne 0 then
      as := [j eq i select 1 
                    else -Round(MonomialCoefficient(cubic, z[i]^2*z[j])/(3*mc))
              : j in [1..3]];
      tr := Matrix([ii eq i select as
                            else [j eq ii select 1 else 0 : j in [1..3]]
                     : ii in [1..3]]);
      cubic := PolMat(cubic, tr);
      trans := trans*tr;
    end if;                             
  end for;
  return cubic, trans;
end function;

function redrep(C)
  trans := IdentityMatrix(Integers(), 3);
  s0 := siz(C);
  s1 := siz(C);
  C1 := C;
  tr1 := trans;
  repeat
    s := s1;
    C := C1;
    trans := tr1;
    C1, tr := redstep(C);
    tr1 := trans*tr;
    s1 := siz(C1);
  until s1 ge s;
  return C, trans;
end function;

function em(i,j,k) // elementary matrix
  return IdentityMatrix(Integers(),3)
          + Matrix(3, [(ii eq i and jj eq j) select k else 0
                          : ii,jj in [1..3]]);
end function;

// set up a list of (monoid) generators
gs1 := [em(i,j,1) : i,j in [1..3] | i ne j];
gs2 := [em(i,j,-1) : i,j in [1..3] | i ne j];
gens1 := gs1 cat gs2;
gens2 := gens1 cat [gs1[i]*gs1[j] : j in [i..#gs1], i in [1..#gs1]]
               cat [gs1[i]*gs2[j] : i,j in [1..#gs1] | i ne j]
               cat [gs2[i]*gs2[j] : j in [i..#gs2], i in [1..#gs2]];


function red0(pol, gens)
  pols := [PolMat(pol, g) : g in gens];
  _, pos := Min([siz(mm) : mm in pols]);
  return pols[pos], gens[pos];
end function;

function red0rep(C)
  trans := IdentityMatrix(Integers(), 3);
  s1 := siz(C);
  C1 := C;
  tr1 := trans;
  dig := Ceiling(Log(10, s1));
  repeat
    s := s1;
    C := C1;
    trans := tr1;
    C1, tr := red0(C, gens2);
    tr1 := trans*tr;
    s1 := siz(C1);
    dig1 := Ceiling(Log(10, s1));
    if dig1 lt dig then
      vprintf TernaryReduce, 2: "  red0rep: size now = %o digits\n",
                                dig1;
      dig := dig1;
    end if;
  until s1 ge s;
  return C, trans;
end function;

function AdHocRedCubic(C)
  trans := IdentityMatrix(Integers(), 3);
  s1 := siz(C);
  C1 := C;
  tr1 := trans;
  dig := Ceiling(Log(10, s1));
  count := 100;
  repeat
    count -:= 1;
    s := s1;
    C := C1;
    trans := tr1;
    C1, tr := redrep(C);
    tr1 := trans*tr;
    C1, tr := red0(C1, gens2);
    tr1 := tr1*tr;
    s1 := siz(C1);
    dig1 := Ceiling(Log(10, s1));
    vprintf TernaryReduce, 3: "   Ad-hoc: gain = %o\n", s-s1;
    if dig1 lt dig then
      vprintf TernaryReduce, 2: "  Ad-hoc cubic red: size now = %o digits\n",
                                dig1;
      dig := dig1;
    end if;
  until s1 ge s or count le 0;
  return C, trans;
end function;

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

intrinsic Reduce(C::RngMPolElt : Al := "Default", Steps := 3, AdHoc := true)
  -> RngMPolElt, AlgMatElt
{Computes an LLL-reduced form in the SL_3(Z)-equivalence class of C.}
  P := Parent(C);
  require Rank(P) eq 3: "C must be a ternary form.";
  deg := TotalDegree(C);
  require forall{m : m in Monomials(C) | TotalDegree(m) eq deg}:
          "Polynomial must be homogeneous.";
  require deg ge 3: "Degree must be at least 3.";
  C0 := C;
  vprintf TernaryReduce, 1: "Reduce: C = %o\n", C;
  P1 := PolynomialRing(BaseRing(P));
  PR := PolynomialRing(RealField());
  trans := IdentityMatrix(Integers(), 3);
  // do some preliminary steps
  if AdHoc then
    repeat
      if deg eq 3 then
        vprintf TernaryReduce, 2:
                "  Preliminary ad-hoc reduction for cubics...\n";
        C1, tr := AdHocRedCubic(C);
        trans := trans*tr;
      else
        vprintf TernaryReduce, 2:
                "  Preliminary ad-hoc reduction for general curves...\n";
        C1, tr := red0rep(C);
        trans := trans*tr;
      end if;
      ff := C1 eq C;
      C := C1;
      vprintf TernaryReduce, 2: "  Ad-hoc reduction...\n";
      deg1 := 3*deg*(deg-2);
      count := 10;
      repeat
        s := siz(C);
        H := TernaryHessian(C);
        vprintf TernaryReduce, 2: "   Computing resultants...\n";
        s1 := siz(C);
        BI1 := Resultant(C, H, P.3);
        BI2 := Resultant(C, H, P.2);
        BI3 := Resultant(C, H, P.1);
        pol12 := Evaluate(BI1, [P1.1,1,0]);
        pol21 := Evaluate(BI1, [1,P1.1,0]);
        pol13 := Evaluate(BI2, [P1.1,0,1]);
        pol31 := Evaluate(BI2, [1,0,P1.1]);
        pol23 := Evaluate(BI3, [0,P1.1,1]);
        pol32 := Evaluate(BI3, [0,1,P1.1]);
        pols := [pol12,pol21,pol13,pol31,pol23,pol32];
        pols := [pol/LeadingCoefficient(pol) : pol in pols];
        // see if complex roots are close together
        eps := 1/1000;
        rtsav := [Degree(pol) lt deg1 select 0 
                                      else -Coefficient(pol,deg1-1)/deg1
                   : pol in pols];
        rtsr := [Round(r) : r in rtsav];
        flags := [rtsr[i] ne 0
                   and forall{j : j in [0..deg1-1] 
                                | Abs(Coefficient(pp,j))
                                   le Binomial(deg1,j)*eps^(deg1-j)}
                       where pp := Evaluate(pols[i], P1.1+rtsav[i])
                   : i in [1..6]];
        mats := [];
        if flags[1] and flags[6] then // affine y = 1
          Append(~mats, Matrix([[1,rtsr[1],0], [0,1,0], [0,rtsr[6],1]]));
        end if;
        if flags[2] and flags[4] then // affine x = 1
          Append(~mats, Matrix([[1,0,0], [rtsr[2],1,0], [rtsr[4],0,1]]));
        end if;
        if flags[3] and flags[5] then // affine z =1
          Append(~mats, Matrix([[1,0,rtsr[3]], [0,1,rtsr[5]], [0,0,1]]));
        end if;
        if not IsEmpty(mats) then
          _, pos := Max([&+[Abs(c) : c in Eltseq(m)] : m in mats]);
         tr := mats[pos];
        elif &or flags then
          if flags[1] then
            Append(~mats, Matrix([[1,rtsr[1],0], [0,1,0], [0,0,1]]));
          end if;
          if flags[2] then
            Append(~mats, Matrix([[1,0,0], [rtsr[2],1,0], [0,0,1]]));
          end if;
          if flags[3] then
            Append(~mats, Matrix([[1,0,rtsr[3]], [0,1,0], [0,0,1]]));
          end if;
          if flags[4] then
            Append(~mats, Matrix([[1,0,0], [0,1,0], [rtsr[4],0,1]]));
          end if;
          if flags[5] then
            Append(~mats, Matrix([[1,0,0], [0,1,rtsr[5]], [0,0,1]]));
          end if;
          if flags[6] then
            Append(~mats, Matrix([[1,0,0], [0,1,0], [0,rtsr[6],1]]));
          end if;
          _, pos := Max([&+[Abs(c) : c in Eltseq(m)] : m in mats]);
         tr := mats[pos];
        else
          break;
        end if;
        ff := false;
        C := PolMat(C, tr);
        trans := trans*tr;
        s1 := siz(C);
        vprintf TernaryReduce, 2: "  Ad-hoc red: Size now %o digits\n",
                                  Ceiling(Log(10, s1));
        if s1 ge s then
          count -:= 1;
        else
          s := s1; count := 5;
        end if;
      until count eq 0;
    until ff;
  end if; //AdHoc
      
  assert C eq PolMat(C0, trans);
  vprintf TernaryReduce, 2: " After ad-hoc reduction, curve is\n%o\n", C;
  // now do the serious work
  for eps in [1.0e-5, 1.0e-10, 1.0e-20][1..Max(1,Min(Steps,3))] do
    vprintf TernaryReduce, 2: " Reduce: Using eps = %o\n", eps;
    cov := Covariant(C : Al := Al, Eps := eps);
    _, tr := LLLGram(cov);
    tr := Transpose(tr);
    vprintf TernaryReduce, 2: " Transformation:\n%o\n", tr;
    C := PolMat(C, tr);
    vprintf TernaryReduce, 2: " New C = %o\n", C;
    trans := trans*tr;
  end for;
  // reduce completely
  mins1 := [];
  c := 16;
  repeat
    mins := mins1;
    repeat
      c *:= 2;
      M := Matrix(3,[Round(c*e) : e in Eltseq(cov)]);
    until IsPositiveDefinite(M);
    _, mins1 := SuccessiveMinima(LatticeWithGram(M));
    mins1 := [Eltseq(v) : v in mins1];
  until mins1 eq mins;
  tr := Transpose(Matrix(mins));
  vprintf TernaryReduce, 2: " Final transformation:\n%o\n", tr;
  C := PolMat(C, tr);
  vprintf TernaryReduce, 2: " Final C = %o\n", C;
  return C, trans*tr;
end intrinsic;

// Now try general ternary forms.
// Use the flexes to find a covariant point geometrically

// Substitute for missing Exp(mat)
function SimpleMatExp(mat)
  n := NumberOfRows(mat);
  e := Parent(mat)!IdentityMatrix(RealField(), n);
  mat -:= 1/n*Trace(mat)*e;
  norm := Sqrt(&+[Norm(c) : c in Eltseq(mat)]);
  b := Ceiling(10*norm);
  mat *:= 1/b;
  m := mat;
  for j := 1 to 10 do
    e +:= m;
    m *:= 1/(j+1)*mat;
  end for;
  e +:= m;
  e := e^b;
  e *:= (Determinant(e))^(-1/n);
  return e;
end function;

intrinsic Covariant(ptseqr::SeqEnum, ptseqc::SeqEnum 
                     : Steps := 0, Eps := 1.0e-20, Al := "Newton")
   -> AlgMatElt, FldPrElt
{Given two sequences of real, resp. complex, points, represented as
 n-tuples of coordinates, not all equal to zero, finds the covariant
 positice definite Hermitian form as a symmetric real matrix. The collection
 of points (real + complex + conjugate complex) must be balanced.
 (I.e., less than d/n of the points may be in a d-dimensional subspace
 for d = 1, ..., n-1. This is not checked; if not satisfied, teh algorithm
 will not terminate.)}
  s := {#pt : pt in ptseqr cat ptseqc};
  require #s eq 1: "All points must have the same number of coordinates.";
  n := Rep(s);
  nr := #ptseqr; nc := #ptseqc;
  npts := nr + 2*nc;
  require nr + 2*nc gt n: "Need at least %o points.", n+1;
  require Al in {"Newton", "SteepestDescent", "Default"}:
          "Al must be one of \"Default\", \"Newton\" and \"SteepestDescent\"";
  vprintf TernaryReduce, 1: 
          "Covariant of points:\n Real: %o\n Complex: %o\n", ptseqr, ptseqc;
  R := RealField();
  // work at I
  trans := IdentityMatrix(R, n);
  M := Parent(trans);
  theta := func<ptseqr, ptseqc |
                &+[R | Log(&+[pt[i]^2 : i in [1..n]]) : pt in ptseqr]
                 + 2*&+[R | Log(&+[Norm(pt[i]) : i in [1..n]]) : pt in ptseqc]>;
  norm := func<pt | &+[Norm(pt[i]) : i in [1..n]]>;
  grad := func<ptseqr, ptseqc |
               &+[M | Matrix([[pt[i]*pt[j] : j in [1..n]] : i in [1..n]])
                       * 1/norm(pt) : pt in ptseqr]
                + 2*&+[M | Matrix([[Real(ComplexConjugate(pt[i])*pt[j])
                                    : j in [1..n]] : i in [1..n]]) 
                            * 1/norm(pt)
                            : pt in ptseqc]
               - npts/n * IdentityMatrix(R, n)>;
  matexp := func<mat | SimpleMatExp(mat)>;
  newseqs := func<ptseqr, ptseqc, mat |
                  [Eltseq(Vector([pt[i] : i in [1..n]])*mat) : pt in ptseqr],
                  [Eltseq(Vector([pt[i] : i in [1..n]])*matc) : pt in ptseqc]
                  where matc := ChangeRing(mat, ComplexField())>;
  // Do a first step if points are very close together
  ptseq := ptseqr cat ptseqc;
  dth := 0;
  if forall{pt : pt in ptseq | pt[3] eq 1} then
    mid := [(&+[pt[i] : pt in ptseqr] + 2*&+[Real(pt[i]) : pt in ptseqc])/npts
              : i in [1,2]];
    rad := Max([Norm(pt[1]-mid[1]) + Norm(pt[2]-mid[2]) : pt in ptseq]);
    if rad lt 0.1 then
      // shift mid to (0,0) and scale to make rad = 1
      vprintf TernaryReduce, 2:
              "  Points are very close to (%o, %o)\n", mid[1], mid[2];
      mrad := (&+[Norm(pt[1]-mid[1]) + Norm(pt[2]-mid[2]) : pt in ptseqr]
                + 2*&+[Norm(pt[1]-mid[1]) + Norm(pt[2]-mid[2]) : pt in ptseqc])
                / npts;
      r := mrad^(-1/6);
      trans := Matrix(R, [[1,0,0], [0,1,0], [-mid[1],-mid[2],1]])
                * DiagonalMatrix([r,r,r^-2]);
      vprintf TernaryReduce, 3: "  Transforming by\n%o\n", trans;
      ptseqr, ptseqc := newseqs(ptseqr, ptseqc, trans);
      ptseqr := [[r^2*e : e in pt] : pt in ptseqr];
      ptseqc := [[r^2*e : e in pt] : pt in ptseqc];
      dth := npts*4*Log(r);
      vprintf TernaryReduce, 3: "  Points are now\n%o\n%o\n", ptseqr, ptseqc;
    end if;
  end if;
  prec := Min([pr : c in &cat[[pt[1],pt[2],pt[3]] : pt in ptseqr cat ptseqc]
                | pr ge 0 where pr := Precision(c)]);
  _, oldprec := HasAttribute(FldPr, "Precision");
  Eps := Max(Eps, Exp(-1.5*Log(10)*prec));
  prec := Max(oldprec, prec);
  AssertAttribute(FldPr, "Precision", prec);
  vprintf TernaryReduce, 2: "  Computing with precision = %o\n", prec;
  // Set up machinery for Newton steps
  eltmat := func< i,j | Matrix(R, n, 
                               [ii eq i and jj eq j select 1/Sqrt(2) else 0
                                : jj, ii in [1..n]]) >;
  tr := Orthonormalize(Matrix(R, n-1,
                              [i eq j select 2 else 1 : i,j in [2..n]]));
  b1mat := tr^-1*Matrix(R, [[i eq 1 select 1 else i eq j select -1 else 0
                              : i in [1..n]] : j in [2..n]]);
  basis := [DiagonalMatrix([b1mat[i,j] : j in [1..n]]) : i in [1..n-1]]
            cat [eltmat(i,j) + eltmat(j,i) : j in [i+1..n], i in [1..n]];
  // basis is an orthonormal basis of the tangent space at I
  ipr := func<pt, mat | &+[pt[i]*mat[i,j]*pt[j] : i,j in [1..n]]>;
  ipc := func<pt, mat | Real(&+[ComplexConjugate(pt[i])*mat[i,j]*pt[j]
                                   : i,j in [1..n]])>;
  function newton(ptseqr, ptseqc)
    // 1st derivative
    der1 := [&+[ipr(pt, b)/norm(pt) : pt in ptseqr]
               + 2*&+[ipc(pt, b)/norm(pt) : pt in ptseqc]
              : b in basis];
    // 2nd derivative
    der2 := Matrix([[&+[ipr(pt,b1*b2)/norm(pt)
                        - (ipr(pt,b1)*ipr(pt,b2))/norm(pt)^2 : pt in ptseqr]
                      + 2*&+[ipc(pt,b1*b2)/norm(pt)
                              - (ipc(pt,b1)*ipc(pt,b2))/norm(pt)^2
                              : pt in ptseqc]
                     : b2 in basis] : b1 in basis]);
    // Newton step
    step := -1/2*Vector(der1)*der2^-1;
    gr := &+[der1[i]*basis[i] : i in [1..#basis]]; // gradient
    return gr, &+[step[i]*basis[i] : i in [1..#basis]];
  end function;
  // Try Newton steps first, do steepest descent if that doesn't reduce theta.
  c := Exp(0); // get 1.0 to current precision
  oldtheta := theta(ptseqr, ptseqc);
  vprintf TernaryReduce, 2: "  Initial theta = %o\n", oldtheta-dth;
  repeat
    // Newton first
    flag := false;
    if Al eq "Newton" then
      gr, step := newton(ptseqr, ptseqc);
      gr := -1/2*(gr + Transpose(gr));
      vprintf TernaryReduce, 2: "  Gradient = \n%o\n", gr;
      step := 1/2*(step + Transpose(step));
      vprintf TernaryReduce, 3: "  Newton step =\n%o\n", step;
      if &+[Norm(c) : c in Eltseq(step)] gt 10 then
        vprintf TernaryReduce, 3: "  Step too large -- try steepest descent.\n";
        flag := true;
      else
        mat := matexp(step);
        mat := 1/2*(mat + Transpose(mat));
        vprintf TernaryReduce, 3: "  Exponentiated =\n%o\n", mat;
        psr1, psc1 := newseqs(ptseqr, ptseqc, mat);
        th1 := theta(psr1, psc1);
        if th1 lt oldtheta then
          vprintf TernaryReduce, 3: "  Newton step successful.\n";
        else
          vprintf TernaryReduce, 3:
                  "  Newton step unsuccessful -- try steepest descent.\n";
          flag := true;
        end if;
      end if;
    else
      gr := -grad(ptseqr, ptseqc);
      vprintf TernaryReduce, 2: "  Gradient = \n%o\n", gr;
      flag := true;
    end if;
    if flag then
      mat := matexp(c*gr);
      mat := 1/2*(mat + Transpose(mat));
      psr1, psc1 := newseqs(ptseqr, ptseqc, mat);
      th1 := theta(psr1, psc1);
      while th1 gt oldtheta do
        c /:= 2;
        vprintf TernaryReduce, 3: "  th1 > th0 : reduce c to %o\n", c;
        mat := matexp(c*gr);
        mat := 1/2*(mat + Transpose(mat));
        psr1, psc1 := newseqs(ptseqr, ptseqc, mat);
        th1 := theta(psr1, psc1);
      end while;
      mat2 := matexp(c/2*gr);
      mat2 := 1/2*(mat2 + Transpose(mat2));
      psr2, psc2 := newseqs(ptseqr, ptseqc, mat2);
      th2 := theta(psr2, psc2);
      if th2 lt th1 then
        th1 := th2;
        c /:= 2;
        psr1 := psr2; psc1 := psc2;
        mat := mat2;
      else
        mat2 := matexp(c*2*gr);
        mat2 := 1/2*(mat2 + Transpose(mat2));
        psr2, psc2 := newseqs(ptseqr, ptseqc, mat2);
        th2 := theta(psr2, psc2);
        if th2 lt th1 then
          th1 := th2;
          c *:= 2;
          psr1 := psr2; psc1 := psc2;
          mat := mat2;
        end if;
      end if;
      vprintf TernaryReduce, 3: "  c now = %o\n", c;
    end if;
    oldtheta := th1;
    vprintf TernaryReduce, 2: "  Theta now = %o\n", th1 - dth;
    ptseqr := psr1; ptseqc := psc1;
    vprintf TernaryReduce, 3: "  Transform points by:\n%o\n", mat;
    trans := trans*mat;
    vprintf TernaryReduce, 3: "  New total transformation:\n%o\n", trans;
  until &+[c^2 : c in Eltseq(gr)] lt Eps or c^2 lt Eps;
  AssertAttribute(FldPr, "Precision", oldprec);
  herm := trans*Transpose(trans);
  vprintf TernaryReduce, 2: "Covariant is\n%o\n", herm;
  vprintf TernaryReduce, 2: "Theta = %o\n", oldtheta - dth;
  return herm, oldtheta - dth;
end intrinsic;

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

// Ueberschiebung(F,G,k)
//  = c*det([[d/dx1 d/dx2] [d/dy1 d/dy2]])^k
//       F(x1,y1) G(x2,y2)|(x1=x2=x, y1=y2=y)
// Can generalise this to three ternary forms:
// TernaryUeberschiebung(F_1,F_2,F_3, k)
//  = c*det((d/dx_ij))^k F_1(x_11, x_12, x_13) F_2() F_3()|(x_ij = x_j)

intrinsic TernaryUeberschiebung(F1::RngMPolElt, F2::RngMPolElt, F3::RngMPolElt, k::RngIntElt) -> RngMPolElt
{Computes the k-th Ueberschiebung of three ternary forms.}
  P3 := Parent(F1);
  require Parent(F2) eq P3 and Parent(F3) eq P3:
          "TernaryUeberschiebung(",F1,", ",F2,", ",F3,", ",k,"): ",
          "Forms must live in the same polynomial ring.";
  require Rank(P3) eq 3:
          "TernaryUeberschiebung(",F1,", ",F2,", ",F3,", ",k,"): ",
          "Forms must be ternary.";
  x0 := [P3.1, P3.2, P3.3];
  n1 := TotalDegree(F1); n2 := TotalDegree(F2); n3 := TotalDegree(F3);
  require forall{m : m in Monomials(F1) | TotalDegree(m) eq n1}
           and forall{m : m in Monomials(F2) | TotalDegree(m) eq n2}
           and forall{m : m in Monomials(F3) | TotalDegree(m) eq n3}:
          "TernaryUeberschiebung(",F1,", ",F2,", ",F3,", ",k,"): ",
          "Polynomials must be homogeneous.";
  if n1 lt k or n2 lt k or n3 lt k then return P3!0; end if;
  P9<[x]> := PolynomialRing(CoefficientRing(P3), 9);
  td := func<pol, i,j,k | Derivative(Derivative(Derivative(pol, i), j), k)>;
  op := func<pol | td(pol, 1,5,9) + td(pol, 2,6,7) + td(pol, 3,4,8)
                    - td(pol, 3,5,7) - td(pol, 2,4,9) - td(pol, 1,6,8)>;
  res := Evaluate(F1, x[1..3])*Evaluate(F2, x[4..6])*Evaluate(F3, x[7..9]);
  for i := 1 to k do
    res := op(res);
  end for;
  coef := 1/(&*[n : n in [n1-k+1..n1]]
              * &*[n : n in [n2-k+1..n2]]
              * &*[n : n in [n3-k+1..n3]]);
  return (CoefficientRing(P3)!coef)*Evaluate(res, x0 cat x0 cat x0);
end intrinsic;

intrinsic TernaryHessian(F::RngMPolElt) -> RngMPolElt
{Compte the Hessian of a ternary form.}
  P3 := Parent(F);
  require Rank(P3) eq 3:
          "TernaryHessian(",F,"): Form must be ternary.";
  n := TotalDegree(F);
  require forall{m : m in Monomials(F) | TotalDegree(m) eq n}:
          "TernaryHessian(",F,"): Polynomial must be homogeneous.";
  dd := func<pol, i,j | Derivative(Derivative(pol, i), j)>;
  dF := [dd(F,1,1), dd(F,2,2), dd(F,3,3),
           dd(F,2,3), dd(F,3,1), dd(F,1,2)];
  return dF[1]*dF[2]*dF[3] + 2*dF[4]*dF[5]*dF[6]
          - dF[1]*dF[4]^2 - dF[2]*dF[5]^2 - dF[3]*dF[6]^2;
end intrinsic;

P10<[a]> := PolynomialRing(Integers(), 10);

S := -(
   -144*a[1]*a[4]*a[8]*a[10] + 48*a[1]*a[4]*a[9]^2 + 216*a[1]*a[5]*a[7]*a[10] - 
    24*a[1]*a[5]*a[8]*a[9] - 144*a[1]*a[6]*a[7]*a[9] + 48*a[1]*a[6]*a[8]^2 + 
    48*a[2]^2*a[8]*a[10] - 16*a[2]^2*a[9]^2 - 144*a[2]*a[3]*a[7]*a[10] + 
    16*a[2]*a[3]*a[8]*a[9] - 24*a[2]*a[4]*a[5]*a[10] + 16*a[2]*a[4]*a[6]*a[9] + 
    8*a[2]*a[5]^2*a[9] - 24*a[2]*a[5]*a[6]*a[8] + 48*a[2]*a[6]^2*a[7] + 
    48*a[3]^2*a[7]*a[9] - 16*a[3]^2*a[8]^2 + 48*a[3]*a[4]^2*a[10] - 
    24*a[3]*a[4]*a[5]*a[9] + 16*a[3]*a[4]*a[6]*a[8] + 8*a[3]*a[5]^2*a[8] - 
    24*a[3]*a[5]*a[6]*a[7] - 16*a[4]^2*a[6]^2 + 8*a[4]*a[5]^2*a[6] - a[5]^4);

T :=
    5832*a[1]^2*a[7]^2*a[10]^2 - 3888*a[1]^2*a[7]*a[8]*a[9]*a[10] + 
    864*a[1]^2*a[7]*a[9]^3 + 864*a[1]^2*a[8]^3*a[10] - 216*a[1]^2*a[8]^2*a[9]^2 
    - 3888*a[1]*a[2]*a[4]*a[7]*a[10]^2 + 1296*a[1]*a[2]*a[4]*a[8]*a[9]*a[10] - 
    288*a[1]*a[2]*a[4]*a[9]^3 + 1296*a[1]*a[2]*a[5]*a[7]*a[9]*a[10] - 
    864*a[1]*a[2]*a[5]*a[8]^2*a[10] + 144*a[1]*a[2]*a[5]*a[8]*a[9]^2 + 
    1296*a[1]*a[2]*a[6]*a[7]*a[8]*a[10] - 864*a[1]*a[2]*a[6]*a[7]*a[9]^2 + 
    144*a[1]*a[2]*a[6]*a[8]^2*a[9] + 1296*a[1]*a[3]*a[4]*a[7]*a[9]*a[10] - 
    864*a[1]*a[3]*a[4]*a[8]^2*a[10] + 144*a[1]*a[3]*a[4]*a[8]*a[9]^2 + 
    1296*a[1]*a[3]*a[5]*a[7]*a[8]*a[10] - 864*a[1]*a[3]*a[5]*a[7]*a[9]^2 + 
    144*a[1]*a[3]*a[5]*a[8]^2*a[9] - 3888*a[1]*a[3]*a[6]*a[7]^2*a[10] + 
    1296*a[1]*a[3]*a[6]*a[7]*a[8]*a[9] - 288*a[1]*a[3]*a[6]*a[8]^3 + 
    864*a[1]*a[4]^3*a[10]^2 - 864*a[1]*a[4]^2*a[5]*a[9]*a[10] - 
    864*a[1]*a[4]^2*a[6]*a[8]*a[10] + 576*a[1]*a[4]^2*a[6]*a[9]^2 + 
    648*a[1]*a[4]*a[5]^2*a[8]*a[10] + 72*a[1]*a[4]*a[5]^2*a[9]^2 + 
    1296*a[1]*a[4]*a[5]*a[6]*a[7]*a[10] - 720*a[1]*a[4]*a[5]*a[6]*a[8]*a[9] - 
    864*a[1]*a[4]*a[6]^2*a[7]*a[9] + 576*a[1]*a[4]*a[6]^2*a[8]^2 - 
    540*a[1]*a[5]^3*a[7]*a[10] - 36*a[1]*a[5]^3*a[8]*a[9] + 
    648*a[1]*a[5]^2*a[6]*a[7]*a[9] + 72*a[1]*a[5]^2*a[6]*a[8]^2 - 
    864*a[1]*a[5]*a[6]^2*a[7]*a[8] + 864*a[1]*a[6]^3*a[7]^2 + 
    864*a[2]^3*a[7]*a[10]^2 - 288*a[2]^3*a[8]*a[9]*a[10] + 64*a[2]^3*a[9]^3 - 
    864*a[2]^2*a[3]*a[7]*a[9]*a[10] + 576*a[2]^2*a[3]*a[8]^2*a[10] - 
    96*a[2]^2*a[3]*a[8]*a[9]^2 - 216*a[2]^2*a[4]^2*a[10]^2 + 
    144*a[2]^2*a[4]*a[5]*a[9]*a[10] + 144*a[2]^2*a[4]*a[6]*a[8]*a[10] - 
    96*a[2]^2*a[4]*a[6]*a[9]^2 + 72*a[2]^2*a[5]^2*a[8]*a[10] - 
    48*a[2]^2*a[5]^2*a[9]^2 - 864*a[2]^2*a[5]*a[6]*a[7]*a[10] + 
    144*a[2]^2*a[5]*a[6]*a[8]*a[9] + 576*a[2]^2*a[6]^2*a[7]*a[9] - 
    216*a[2]^2*a[6]^2*a[8]^2 - 864*a[2]*a[3]^2*a[7]*a[8]*a[10] + 
    576*a[2]*a[3]^2*a[7]*a[9]^2 - 96*a[2]*a[3]^2*a[8]^2*a[9] + 
    144*a[2]*a[3]*a[4]^2*a[9]*a[10] - 720*a[2]*a[3]*a[4]*a[5]*a[8]*a[10] + 
    144*a[2]*a[3]*a[4]*a[5]*a[9]^2 + 1296*a[2]*a[3]*a[4]*a[6]*a[7]*a[10] - 
    48*a[2]*a[3]*a[4]*a[6]*a[8]*a[9] + 648*a[2]*a[3]*a[5]^2*a[7]*a[10] - 
    24*a[2]*a[3]*a[5]^2*a[8]*a[9] - 720*a[2]*a[3]*a[5]*a[6]*a[7]*a[9] + 
    144*a[2]*a[3]*a[5]*a[6]*a[8]^2 + 144*a[2]*a[3]*a[6]^2*a[7]*a[8] + 
    144*a[2]*a[4]^2*a[5]*a[6]*a[10] - 96*a[2]*a[4]^2*a[6]^2*a[9] - 
    36*a[2]*a[4]*a[5]^3*a[10] - 24*a[2]*a[4]*a[5]^2*a[6]*a[9] + 
    144*a[2]*a[4]*a[5]*a[6]^2*a[8] - 288*a[2]*a[4]*a[6]^3*a[7] + 
    12*a[2]*a[5]^4*a[9] - 36*a[2]*a[5]^3*a[6]*a[8] + 72*a[2]*a[5]^2*a[6]^2*a[7] 
    + 864*a[3]^3*a[7]^2*a[10] - 288*a[3]^3*a[7]*a[8]*a[9] + 64*a[3]^3*a[8]^3 + 
    576*a[3]^2*a[4]^2*a[8]*a[10] - 216*a[3]^2*a[4]^2*a[9]^2 - 
    864*a[3]^2*a[4]*a[5]*a[7]*a[10] + 144*a[3]^2*a[4]*a[5]*a[8]*a[9] + 
    144*a[3]^2*a[4]*a[6]*a[7]*a[9] - 96*a[3]^2*a[4]*a[6]*a[8]^2 + 
    72*a[3]^2*a[5]^2*a[7]*a[9] - 48*a[3]^2*a[5]^2*a[8]^2 + 
    144*a[3]^2*a[5]*a[6]*a[7]*a[8] - 216*a[3]^2*a[6]^2*a[7]^2 - 
    288*a[3]*a[4]^3*a[6]*a[10] + 72*a[3]*a[4]^2*a[5]^2*a[10] + 
    144*a[3]*a[4]^2*a[5]*a[6]*a[9] - 96*a[3]*a[4]^2*a[6]^2*a[8] - 
    36*a[3]*a[4]*a[5]^3*a[9] - 24*a[3]*a[4]*a[5]^2*a[6]*a[8] + 
    144*a[3]*a[4]*a[5]*a[6]^2*a[7] + 12*a[3]*a[5]^4*a[8] - 
    36*a[3]*a[5]^3*a[6]*a[7] + 64*a[4]^3*a[6]^3 - 48*a[4]^2*a[5]^2*a[6]^2 + 
    12*a[4]*a[5]^4*a[6] - a[5]^6;

Delta := ExactQuotient(S^3 - T^2, 1728);

intrinsic T3Invariants(F::RngMPolElt) -> RngElt, RngElt, RngElt
{Compute the invariants S, T and Delta of a ternary cubic.}
  P3 := Parent(F);
  require Rank(P3) eq 3:
          "T3Invariants(",F,"): Form must be ternary.";
  require forall{m : m in Monomials(F) | TotalDegree(m) eq 3}:
          "T3Invariants(",F,"): Polynomial must be a cubic form.";
  mons := [P3.1^i*P3.2^j*P3.3^(3-i-j) : j in [3-i..0 by -1], i in [3..0 by -1]];
  cofs := [MonomialCoefficient(F, m) : m in mons];
  SF := Evaluate(S, cofs); TF := Evaluate(T, cofs);
  if Characteristic(P3) notin {2, 3} then
    DF := ExactQuotient(SF^3 - TF^2, Parent(SF)!1728);
  else
    DF := Evaluate(Delta, cofs);
  end if;
  return SF, TF, DF;
end intrinsic;

/////////////////////////////////////////////////////////////////////
// Minimisation of cubics
/////////////////////////////////////////////////////////////////////

intrinsic Minimise(C::RngMPolElt) -> RngMPolElt, AlgMatElt
{Minimise the integral ternary cubic form C.
 Returns the minimised form and the transformation.}
  P3 := Parent(C);
  require Rank(P3) eq 3: "Form must be ternary.";
  require forall{m : m in Monomials(C) | TotalDegree(m) eq 3}:
          "Polynomial must be a cubic form.";
  require BaseRing(P3) cmpeq Rationals():
          "Form must live in polynomial ring over the rationals.";
  require forall{c : c in Coefficients(C) | Denominator(c) eq 1}:
          "Form must have integral coefficients.";
  vprintf Minimise, 1: "Minimise C =\n%o\n", C;
  trans := IdentityMatrix(Rationals(), 3);
  S, T, D := T3Invariants(C);
  g := GCD([Integers() | S, T, D]);
  vprintf Minimise, 2: "  Factoring GCD of invariants ...\n";
  plist := [p : p in PrimeDivisors(g)
              | (S eq 0 or Valuation(S, p) ge 4)
                 and (T eq 0 or Valuation(T, p) ge 6)
                 and (D eq 0 or Valuation(D, p) ge 12)];
  vprintf Minimise, 2: "  primes where model is nonminimal:\n    %o\n",
                           plist;
  for p in plist do
    C, tr := pMinimise(C, p);
    trans := trans*tr;
  end for;
  vprintf Minimise, 1: "Minimise -- done.\n";
  return C, trans;
end intrinsic;

intrinsic pMinimise(C::RngMPolElt, p::RngIntElt) -> RngMPolElt, AlgMatElt
{Minimise the integral ternary cubic Form C at the prime p.
 Returns the minimised form and the transformation.}
  P3 := Parent(C);
  require Rank(P3) eq 3: "Form must be ternary.";
  require forall{m : m in Monomials(C) | TotalDegree(m) eq 3}:
          "Polynomial must be a cubic form.";
  require BaseRing(P3) cmpeq Rationals():
          "Form must live in polynomial ring over the rationals.";
  require forall{c : c in Coefficients(C) | Denominator(c) eq 1}:
          "Form must have integral coefficients.";
  vprintf Minimise, 1: "pMinimise at p = %o ...\n", p;
  C0 := C; // save for later check
  trans := IdentityMatrix(Rationals(), 3);
  indices := [<i, j, 3-i-j> : j in [0..3-i], i in [0..3]];
  mons := [P3.1^(a[1])*P3.2^(a[2])*P3.3^(a[3]) : a in indices];
  MC := MonomialCoefficient;
  F := GF(p);
  PF := PolynomialRing(F, 3);
  S, T, D := T3Invariants(C);
  steps := Min([Floor(Valuation(pair[1],p)/pair[2]) 
                 : pair in [<S,4>, <T,6>, <D,12>]
                 | pair[1] ne 0]);
  if p eq 3 then // check Kraus's conditions
    Tred := Integers()!(T/3^(6*steps));
    if Tred mod 27 in {9, 18} then
      steps -:= 1;
    end if; 
  elif p eq 2 then // check Kraus's conditions
    Tred := Integers()!(T/2^(6*steps));
    if Tred mod 4 ne 3
        and (Valuation(S, 2) lt 4*steps+4 or Tred mod 32 notin {0, 8})
    then
      steps -:= 1;
    end if;
  end if;
  // now steps is the number of minimisation steps to make
  vprintf Minimise, 2: "  %o minimisation steps to make.\n", steps;
  while steps gt 0 do
    vprintf Minimise, 3: "    remaining steps: %o\n", steps;
    CF := PF!C;
    if CF eq 0 then
      C := 1/p*C;
      vprintf Minimise, 2:
              "    C mod %o is zero. Dividing C by %o.\n", p, p;
      steps -:= 1;
    else
      fact := Factorisation(CF);
      if forall(f){f : f in fact | f[2] eq 1} then
        vprintf Minimise, 3: "    C mod %o is squarefree.\n", p;
        I := ideal<PF | CF, [Derivative(CF, i) : i in [1..3]]>;
        assert Dimension(I) eq 1; // isolated singularities
        // we know that there is one singular point; move it to (0:0:1)
        V1 := Variety(ideal<PF | I, PF.3-1>);
        if IsEmpty(V1) then
          V2 := Variety(ideal<PF | I, PF.2-1>);
          if IsEmpty(V2) then
            V3 := Variety(ideal<PF | I, PF.1-1>);
            assert not IsEmpty(V3);
            sing := V3[1]; // = (1:0:0)
            C := Evaluate(C, [P3.3, P3.2, -P3.1]);
            trans := trans*Matrix(Rationals(), [[0,0,1], [0,1,0], [-1,0,0]]);
          else // sing = (x:1:0)
            sing := V2[1];
            C := Evaluate(C, [P3.1 + (Integers()!sing[1])*P3.3, P3.3, -P3.2]);
            trans := trans*Matrix(Rationals(),
                                [[1,0,Integers()!sing[1]], [0,0,1], [0,-1,0]]);
          end if;
        else
          sing := V1[1];
          C := Evaluate(C, [P3.1 + (Integers()!sing[1])*P3.3,
                            P3.2 + (Integers()!sing[2])*P3.3, 
                            P3.3]);
          trans := trans*Matrix(Rationals(),
                                [[1,0,Integers()!sing[1]],
                                 [0,1,Integers()!sing[2]],
                                 [0,0,1]]);
        end if;
        CF := PF!C;
        // look at degree 2 part in X,Y
        deg2 := Evaluate(CF, [PF.1, PF.2, 1]) - Evaluate(CF, [PF.1, PF.2, 0]);
        if deg2 eq 0 then
          // CF only depends on X,Y, so describes three lines through one point
          // since CF is squarefree, the lines are distinct
          // if minimisation is possible, then by z --> z/p
          vprintf Minimise, 3:
                  "    C mod %o is union of three distinct lines.\n", p;
          if forall{i : i in [1..#indices]
                      | Valuation(MC(C, mons[i]), p) ge indices[i,3]}
          then
            C := Evaluate(C, [P3.1, P3.2, 1/p*P3.3]);
            trans := trans*Matrix(Rationals(), [[1,0,0], [0,1,0], [0,0,1/p]]);
            steps -:= 1;
          else
            break;
          end if;
        else // have nonzero degree 2 part
          // degree 2 part must be a square
          vprintf Minimise, 3:
                  "    C mod %o is cuspidal cubic or line touching conic.\n", p;
          a := MC(deg2, PF.1^2);
          b := MC(deg2, PF.1*PF.2);
          c := MC(deg2, PF.2^2);
          assert 4*a*c eq b^2;
          // transform to Y^2
          if c ne 0 then
            if p eq 2 then
              s := a/c;
            else
              s := b/(2*c);
            end if;
            C := Evaluate(C, [P3.1, P3.2 - (Integers()!s)*P3.1, P3.3]);
            trans := trans*Matrix(Rationals(),
                                  [[1,0,0], [-Integers()!s,1,0], [0,0,1]]);
          else // c eq 0
            C := Evaluate(C, [P3.2, -P3.1, P3.3]);
            trans := trans*Matrix(Rationals(), [[0,1,0], [-1,0,0], [0,0,1]]);
          end if; // c ne 0
          CF := PF!C;
          // now CF = c*Y^2*Z + F(X,Y) must be either a cuspidal cubic
          //  or a line touching a conic; in both cases we must have
          //  p^2 | a003, a102
          assert Valuation(MC(C, P3.3^3), p) ge 2
                  and Valuation(MC(C, P3.1*P3.3^2), p) ge 2;
          vprintf Minimise, 3:
                  "      transform to another model of same level.\n";
          C := p*Evaluate(C, [P3.1, P3.2, 1/p*P3.3]);
          trans := trans*Matrix(Rationals(), [[1,0,0], [0,1,0], [0,0,1/p]]);
          // this leads to CF being a product of linear factors
        end if; // deg2 eq 0
      else // a multiple line
        // f = <line, 2 or 3>
        // move multiple line to Y
        line := f[1];
        a,b,c := Explode([MC(line, PF.i) : i in [1..3]]);
        if b ne 0 then
          C := Evaluate(C, [P3.1, 
                            P3.2 - (Integers()!(a/b))*P3.1
                                 - (Integers()!(c/b))*P3.3,
                            P3.3]);
          trans := trans*Matrix(Rationals(),
                                [[1,0,0],
                                 [-Integers()!(a/b),1,-Integers()!(c/b)],
                                 [0,0,1]]);
        elif c ne 0 then // line = aX + cZ
          C := Evaluate(C, [P3.1, -P3.3,
                            P3.2 - (Integers()!(a/c))*P3.1]);
          trans := trans*Matrix(Rationals(),
                                [[1,0,0], [0,0,-1], [-Integers()!(a/c),1,0]]);
        else // line = aX
          C := Evaluate(C, [P3.2, -P3.1, P3.3]);
          trans := trans*Matrix(Rationals(), [[0,1,0], [-1,0,0], [0,0,1]]);
        end if;
        CF := PF!C;
        if f[2] eq 2 then // a double line and a simple one
          // move to Y^2 Z
          vprintf Minimise, 3:
             "    C mod %o is the union of a double line and a simple line.\n",
             p;
          line := ExactQuotient(CF, PF.2^2);
          a,b,c := Explode([MC(line, PF.i) : i in [1..3]]);
          if c ne 0 then
            C := Evaluate(C, [P3.1, P3.2,
                              P3.3 - (Integers()!(a/c))*P3.1
                                   - (Integers()!(b/c))*P3.2]);
            trans := trans*Matrix(Rationals(), 
                                  [[1,0,0], [0,1,0],
                                   [-Integers()!(a/c),-Integers()!(b/c),1]]);
          else
            assert a ne 0;
            C := Evaluate(C, [P3.3, P3.2, -P3.1]);
            trans := trans*Matrix(Rationals(), [[0,0,1], [0,1,0], [-1,0,0]]);
          end if;
          CF1 := PF!(1/p*&+[MC(C, P3.1^i*P3.3^(3-i))*P3.1^i*P3.3^(3-i)
                             : i in [0..3]]);
          if MC(CF1, PF.1^3) eq 0 then
            if forall{i : i in [1..#indices]
                        | Valuation(MC(C, mons[i]), p) ge indices[i,1]}
            then
              C := Evaluate(C, [1/p*P3.1, P3.2, P3.3]);
              trans := trans*Matrix(Rationals(), [[1/p,0,0], [0,1,0], [0,0,1]]);
              steps -:= 1;
            elif forall{i : i in [1..#indices]
                          | Valuation(MC(C, mons[i]), p) 
                             ge indices[i,1]+indices[i,3]-1}
            then
              C := p*Evaluate(C, [1/p*P3.1, P3.2, 1/p*P3.3]);
              trans := trans*Matrix(Rationals(), 
                                    [[1/p,0,0], [0,1,0], [0,0,1/p]]);
              steps -:= 1;
            else
              break;
            end if;
          else // MC(CF1, PF.1^3) ne 0
            // CF1 must be a cube
            fact := Factorisation(CF1);
            if fact[1,2] ne 3 then break; end if;
            r := -MC(fact[1,1], PF.3)/MC(fact[1,1], PF.1); // triple root
            C := Evaluate(C, [P3.1 + (Integers()!r)*P3.3, P3.2, P3.3]);
            trans := trans*Matrix(Rationals(), 
                                  [[1,0,Integers()!r], [0,1,0], [0,0,1]]);
            if forall{i : i in [1..#indices]
                        | Valuation(MC(C, mons[i]), p) ge 1-indices[i,2]}
            then
              vprintf Minimise, 3:
                      "      transform to another model of same level.\n";
              C := 1/p*Evaluate(C, [P3.1, p*P3.2, P3.3]);
              trans := trans*Matrix(Rationals(), [[1,0,0], [0,p,0], [0,0,1]]);
              // makes CF a triple line
            else
              break;
            end if;
          end if; // MC(CF1, PF.1^3) eq 0
        else // a triple line
          // CF = Y^3
          vprintf Minimise, 3:
                  "    C mod %o is a triple line.\n", p;
          CF1 := PF!(1/p*&+[MC(C, P3.1^i*P3.3^(3-i))*P3.1^i*P3.3^(3-i)
                             : i in [0..3]]);
          // this must have a multiple root
          if CF1 eq 0 then
            C := 1/p^2*Evaluate(C, [P3.1, p*P3.2, P3.3]);
            trans := trans*Matrix(Rationals(), [[1,0,0], [0,p,0], [0,0,1]]);
            steps -:= 1;
          else
            fact := Factorisation(CF1);
            if not exists(f){f : f in fact | f[2] ge 2} then break; end if;
            lin := f[1];
            a := MC(lin, PF.1);
            b := MC(lin, PF.3);
            if a ne 0 then
              C := Evaluate(C, [P3.1 - (Integers()!(b/a))*P3.3, P3.2, P3.3]);
              trans := trans*Matrix(Rationals(),
                                   [[1,0,-Integers()!(b/a)], [0,1,0], [0,0,1]]);
            else
              C := Evaluate(C, [-P3.3, P3.2, P3.1]);
              trans := trans*Matrix(Rationals(), [[0,0,-1], [0,1,0], [1,0,0]]);
            end if;
            if forall{i : i in [1..#indices]
                        | Valuation(MC(C, mons[i]), p) ge indices[i,3]}
            then
              C := Evaluate(C, [P3.1, P3.2, 1/p*P3.3]);
              trans := trans*Matrix(Rationals(), [[1,0,0], [0,1,0], [0,0,1/p]]);
              steps -:= 1;
            else
              break;
            end if;
          end if;
        end if; // if f[2] eq 2
      end if; // forall(f){f : f in fact | f[2] eq 1}
    end if;
    assert IsDivisibleBy(Evaluate(C0, [&+[trans[i,j]*P3.j : j in [1..3]]
                                        : i in [1..3]]),
                         C);
  end while;
  if steps gt 0 then
    vprintf Minimise, 1:
            "pMinimise at p = %o: not completely minimised.\n", p;
  end if;
  vprintf Minimise, 1: "pMinimise at p = %o -- done.\n", p;
  l := LCM([Denominator(c) : c in Eltseq(trans)]);
  trans *:= l;
  g := GCD([Numerator(c) : c in Eltseq(trans)]);
  trans *:= 1/g;
  return C, trans;
end intrinsic;

// Convert (E, pt) --> (C, pt') where E is an elliptic curve/Q with
// a rational point pt and C is the corresponding plane cubic, with a
// rational point pt' that maps to pt on E.
intrinsic CubicFromPoint(E::CrvEll, P::PtEll) -> RngMPolElt, SeqEnum
{Given an elliptic curve and a point on it, find the corresponding
 3-covering as a ternary cubic, together with a point on it that
 maps to the given point on E.}
  // Bring (E, P) into the form
  // E : y(y-b) = (x-a)(x^2+c), P = (a,0)
  R := CoefficientRing(E);
  require Characteristic(R) ne 2:
          "This does not work in characteristic 2.";
  a1,a2,a3,a4,a6 := Explode(aInvariants(E));
  require P[3] ne 0: "Point may not be at infinity.";
  px := P[1]/P[3]; py := P[2]/P[3];
  P<X,Y> := PolynomialRing(R, 2);
  pol := Y^2 + a1*X*Y + a3*Y - (X^3 + a2*X^2 + a4*X + a6);
  // shift point to (0,0)
  pol := Evaluate(pol, [X+px, Y+py]);
  assert Evaluate(pol, [0,0]) eq 0;
  // get rid of the xy term
  pol := Evaluate(pol, [X, Y-1/(R!2)*MonomialCoefficient(pol,X*Y)*X]);
  assert MonomialCoefficient(pol, X*Y) eq 0 and Evaluate(pol, [0,0]) eq 0;
  // shift x to get RHS into form (x-a)(x^2+c) (from x(x^2+a2*x+a4))
  a := -1/(R!2)*MonomialCoefficient(pol, X^2);
  pol := Evaluate(pol, [X-a, Y]);
  assert MonomialCoefficient(ExactQuotient(Coefficient(pol,Y,0),(X-a)),X) eq 0;
  assert Evaluate(pol, [a,0]) eq 0;
  b := -MonomialCoefficient(pol,Y);
  c := -MonomialCoefficient(pol,X);
  P3<x,y,z> := PolynomialRing(R, 3);
  // cubic and point
  C := y*z^2 + x^2*z + a*x^2*y + b*x*y^2 + c*y^3;
  ps := [R | 0,0,1];
  // if R = rationals, then minimise and reduce
  if R cmpeq Rationals() then
    l := LCM([Denominator(a), Denominator(b), Denominator(c)]);
    C *:= l;
    C, tr := Minimise(C);
    ps := Eltseq(tr^-1*Matrix(R, [[c] : c in ps]));
    C, tr := Reduce(C);
    ps := Eltseq(ChangeRing(tr^-1,R)*Matrix(R, [[c] : c in ps]));
    l := LCM([Denominator(c) : c in ps]);
    ps := [l*c : c in ps];
    g := GCD([Numerator(c) : c in ps]);
    ps := [c/g : c in ps];
    assert Evaluate(C, ps) eq 0;
  end if;
  return C, ps;
end intrinsic;


/*
// Idea: Given that we know by how much we have to minimise the given
//   cubic (it is N = (disc(cubic)/minimal_disc(Jac(cubic)))^(1/12)),
//   we just enter into the computation (excepting p = 2 and 3). If
//   something goes wrong, then because we have an element mod N that
//   is neither zero nor invertible and thus will split N. In this case,
//   split N and repeat with the factors.

forward extendBD;

function split(a,b)
  // printf "split(%o,%o)\n", a, b;
  g := GCD(a,b);
  if g eq 1 then return [c : c in [a,b] | c ne 1]; end if;
  a := ExactQuotient(a,g);
  b := ExactQuotient(b,g);
  return extendBD([a,b], g);
end function;

function extendBD(l, n)
  // printf "extendBD(%o,%o)\n", l, n;
  gl := [GCD(a, n) : a in l];
  l1 := [l[i] : i in [1..#l] | gl[i] eq 1];
  l2 := [l[i] : i in [1..#l] | gl[i] ne 1];
  if IsEmpty(l2) then return [c : c in l cat [n] | c ne 1]; end if;
  s := split(l2[1], n);
  res := l2[2..#l2];
  for a in s do
    res := extendBD(res, a);
  end for;
  return l1 cat res;
end function;

function BD(l)
  res := [];
  for a in l do
    res := extendBD(res, a);
  end for;
  return res;
end function;

function BreakUpWith(n, l)
  // use the numbers in the sequence l to split n into factors
  bl := BD([n] cat [GCD(n,a) : a in l]);
  function val(m, p)
    e := 0;
    while IsDivisibleBy(m, p) do
      e +:= 1;
      m := m div p;
    end while;
    return m, e;
  end function;
  res := [];
  for b in bl do
    n, e := val(n, b);
    Append(~res, <b, e>);
  end for;
  return res;
end function;

intrinsic Minimise(C::RngMPolElt, l::SeqEnum : Check := true)
  -> RngMPolElt, Mtrx
{Minimise C according to the factorisation sequence l (the factors need
 not to be prime).}
  vprintf Minimise, 1: "Minimise(l = %o)\n", l;
  tr := IdentityMatrix(Rationals(), 3);
  for ll in l do
    for i := 1 to ll[2] do
      C, tr1 := Minimise(C, ll[1] : Check := Check);
      tr := tr*tr1;
    end for;
  end for;
  return C, tr;
end intrinsic;

//--------------------------------------------------------------

function val(m, p)
  if m eq 0 then return Infinity(); end if;
  m := Integers()!m;
  p := Integers()!p;
  e := 0;
  while IsDivisibleBy(m ,p) do
    e +:= 1;
    m := m div p;
  end while;
  return e;
end function;

function valmax(m, p, max)
  m := Integers()!m;
  p := Integers()!p;
  e := 0;
  while IsDivisibleBy(m ,p) do
    e +:= 1;
    m := m div p;
    if e ge max then return e; end if;
  end while;
  return e;
end function;

intrinsic MinimiseQuick(C::RngMPolElt) -> RngMPolElt
{Do all minimisation that is possible by scalings onnly.}
  P3 := Parent(C);
  require Rank(P3) eq 3: "Form must be ternary.";
  require forall{m : m in Monomials(C) | TotalDegree(m) eq 3}:
          "Polynomial must be a cubic form.";
  require BaseRing(P3) cmpeq Rationals():
          "Form must live in polynomial ring over the rationals.";
  require forall{c : c in Coefficients(C) | Denominator(c) eq 1}:
          "Form must have integral coefficients.";
  C := 1/GCD(ChangeUniverse(Coefficients(C), Integers()))*C;
  S, T, D := T3Invariants(C);
  S := Integers()!S; T := Integers()!T; D := Integers()!D;
  list := [l[1] : l in BreakUpWith(GCD([S, T, D]), 
                                   ChangeUniverse(Coefficients(C), 
                                                  Integers()))];
  Sort(~list); Reverse(~list);
  mons := MonomialsOfDegree(P3, 3);
  for n in list do
    repeat
      vals := [valmax(MonomialCoefficient(C, m), n, 10) : m in mons];
      inds := [<Degree(m,P3.1), Degree(m,P3.2), Degree(m,P3.3)> : m in  mons];
      // find best possible scaling
      test := func<sy,sz | Min([vals[i] + sy*inds[i,2] + sz*inds[i,3]
                                 : i in [1..#vals]])
                            - (sy+sz)>;
      test1 := function(sy)
                 sz := 0;
                 t0 := test(sy,0);
                 tp := test(sy,1);
                 tm := test(sy,-1);
                 if tp gt t0 then
                   sz := 1;
                   repeat
                     sz +:= 1;
                     t0 := tp;
                     tp := test(sy,sz);
                   until tp le t0;
                   return sz-1, t0;
                 elif tm gt t0 then
                   sz := -1;
                   repeat
                     sz -:= 1;
                     t0 := tm;
                     tm := test(sy,sz);
                   until tm le t0;
                   return sz+1, t0;
                 else
                   return 0, t0;
                 end if;
               end function;
      sy := 0;
      sz,t0 := test1(0);
      szp,tp := test1(1);
      szm,tm := test1(-1);
      if tp gt t0 then
        sy := 1;
        repeat
          sy +:= 1;
          t0 := tp; sz := szp;
          szp; tp := test1(sy);
        until tp le t0;
        sy -:= 1;
      elif tm gt t0 then
        sy := -1;
        repeat
          sy -:= 1;
          t0 := tm; sz := szm;
          szm, tm := test1(sy);
        until tm le t0;
        sy +:= 1;
      end if;
      C := n^-(t0+sy+sz)*Evaluate(C, [P3.1, n^sy*P3.2, n^sz*P3.3]);
      assert forall{c : c in Coefficients(C) | Denominator(c) eq 1};
    until t0 eq 0;
  end for;
  return C;
end intrinsic;

//--------------------------------------------------------------

function Mini300(C, N)
  // here C mod N is const*x^3
  printf "Not yet implemented (Mini300) -- sorry.\n";
  return C;
end function;

function Mini210(C, N)
  // here C mod N is const*x^2*y
  vprintf Minimise, 2: "Mini210(N = %o)\n", N;
  
  return C;
end function;

function Mini2var(C, N)
  // here C mod N does not involve z
  F := Integers(N);
  PF := PolynomialRing(F,3);
  Cred := PF!C;
  vprintf Minimise, 2: "Mini2Var(Cred= %o, N = %o)\n", Cred, N;
  lc := MonomialCoefficient(Cred, PF.1^3);
  if lc eq 0 then // C = y (a x^2 + b x y + c y^2) mod N
    lc := MonomialCoefficient(Cred, PF.1^2*PF.2);
    if lc eq 0 then // C = y^2 (a x + b y) mod N
      lc := MonomialCoefficient(Cred, PF.1*PF.2^2);
      if lc eq 0 then // C = c y^3 mod N
        lc := MonomialCoefficient(Cred, PF.2^3);
        assert IsUnit(lc);
        tr1 := Matrix([[0,-1,0], [1,0,0], [0,0,1]]);
        return Mini300(PolMat(C, tr1), N);
      elif IsUnit(lc) then
        // x <- -y, y <- x - b/a y
        b := MonomialCoefficient(Cred, PF.2^3);
        tr1 := Matrix([[0,-1,0], [1,-Integers()!(b/lc),0], [0,0,1]]);
        return Mini210(PolMat(C, tr1), N);
      else
        return Minimise(C, BreakUpWith(N, [Integers() | lc]) : Check := false);
      end if;
    elif IsUnit(lc) then
      // check discriminant
      b := MonomialCoefficient(Cred, PF.1*PF.2^2);
      c := MonomialCoefficient(Cred, PF.2^3);
      disc := b^2 - 4*lc*c;
      if disc eq 0 then // C = y (...)^2 mod N
        // find linear factor whose square we know
        if IsEven(N) and N ne 2 then
          return Minimise(C, BreakUpWith(N, [Integers() | 2]) : Check := false);
        elif N eq 2 then // b = 0 mod 2, factor is lc*x + c*y = x + c*y
          tr1 := Matrix([[1,-Integers()!(c/lc),0], [0,1,0], [0,0,1]]);
          return Mini210(PolMat(C, tr1), N);
        else // N is odd -> get x + b/(2*lc)*y
          tr1 := Matrix([[1,-Integers()!(b/(2*lc)),0], [0,1,0], [0,0,1]]);
          return Mini210(PolMat(C, tr1), N);
        end if;
      elif IsUnit(disc) then // three distinct roots
        Cnew := PolMat(C, DiagonalMatrix([1,1,1/N]));
        if forall{c : c in Coefficients(Cnew) | Denominator(c) eq 1} then
          return Cnew;
        else
          printf "Minimisation at N = %o not possible.\n", N;
          return C;
        end if;
      else
        return Minimise(C, BreakUpWith(N, [Integers() | disc])
                         : Check := false);
      end if;
    else
      return Minimise(C, BreakUpWith(N, [Integers() | lc]) : Check := false);
    end if;
  elif IsUnit(lc) then
    P1 := PolynomialRing(F);
    C1 := Evaluate(Cred, [P1.1,1,0]);
    a := Coefficient(C1, 2);
    b := Coefficient(C1, 1);
    c := Coefficient(C1, 0);
    disc := Discriminant(C1);
    if disc eq 0 then // multiple root
      DC1 := Derivative(C1);
      if DC1 eq 0 then // must have N eq 3
        assert N eq 3;
        tr1 := Matrix([[1,-Integers()!(c/lc),0], [0,1,0], [0,0,1]]);
        return Mini300(PolMat(C, tr1), N);
      end if;
      disc1 := Discriminant(DC1);
      if disc1 eq 0 then // triple root
        if N eq 3 then // -> (x + c/lc*y)^3
          tr1 := Matrix([[1,-Integers()!(c/lc),0], [0,1,0], [0,0,1]]);
          return Mini300(PolMat(C, tr1), N);
        elif IsDivisibleBy(N, 3) then
          return Minimise(C, BreakUpWith(N, [3]) : Check := false);
        else // N prime to 3 -> (x + a/(3*lc)*y)^3
          tr1 := Matrix([[1,-Integers()!(a/(3*lc)),0], [0,1,0], [0,0,1]]);
          return Mini300(PolMat(C, tr1), N);
        end if;
      elif IsUnit(disc1) then
        if N eq 3 then // DC1 = 2*a*u + b
          assert IsUnit(a);
          dr := -b/(2*a); // double root
          sr := -(a/lc + 2*dr); // single root
          tr1 := Matrix([[1,Integers()!dr,0], [0,1,0], [0,0,1]]);
          sr := sr-dr;
          assert IsUnit(sr);
          tr1 := tr1*Matrix([[1,0,0], [Integers()!(1/sr),1,0], [0,0,1]]);
          return Mini210(PolMat(C, tr1), N);
        elif IsDivisibleBy(N, 3) then
          return Minimise(C, BreakUpWith(N, [3]) : Check := false);        
        else
          DC1 := 1/LeadingCoefficient(DC1)*DC1;
          lin := C1 - lc*P1.1*DC1;
          lin := lin - Coefficient(lin, 2)*DC1;
          ll := Coefficient(lin ,1);
          assert IsUnit(ll);
          dr := -Coefficient(lin,0)/ll; // double root
          sr := -(a/lc + 2*dr); // single root
          tr1 := Matrix([[1,Integers()!dr,0], [0,1,0], [0,0,1]]);
          sr := sr-dr;
          assert IsUnit(sr);
          tr1 := tr1*Matrix([[1,0,0], [Integers()!(1/sr),1,0], [0,0,1]]);
          return Mini210(PolMat(C, tr1), N);
        end if;
      else
        return Minimise(C, BreakUpWith(N, [Integers() | disc1])
                         : Check := false);
      end if;
    elif IsUnit(disc) then // three distinct roots
      Cnew := PolMat(C, DiagonalMatrix([1,1,1/N]));
      if forall{c : c in Coefficients(Cnew) | Denominator(c) eq 1} then
        return Cnew;
      else
        printf "Minimisation at N = %o not possible.\n", N;
        return C;
      end if;
    else
      return Minimise(C, BreakUpWith(N, [Integers() | disc]) : Check := false);
    end if;
  else
    return Minimise(C, BreakUpWith(N, [Integers() | lc]) : Check := false);
  end if;
  return C;
end function;

function MinLine(C, N, line)
  vprintf Minimise, 2: " MinLine(N = %o, line = %o)\n", N, line;
  F := Integers(N);
  PF := PolynomialRing(F,3);
  if line[1] eq 1 then
    line[2] := -line[2]; line[3] :=-line[3];
    tr1 := Matrix(Rationals(), 
                  [ChangeUniverse(line, Integers()), [0,1,0], [0,0,1]]);
  elif line[2] eq 1 then
    line[1] := -line[1]; line[3] :=-line[3];
    tr1 := Matrix(Rationals(), 
                  [[1,0,0], ChangeUniverse(line, Integers()), [0,0,1]]);
    tr1 := tr1*Matrix(Rationals(), [[0,1,0], [-1,0,0], [0,0,1]]);
  elif line[3] eq 1 then
    line[1] := -line[1]; line[2] :=-line[2];
    tr1 := Matrix(Rationals(), 
                  [[1,0,0], [0,1,0], ChangeUniverse(line, Integers())]);
    tr1 := tr1*Matrix(Rationals(), [[0,0,1], [0,1,0], [-1,0,0]]);
  else
    error "line equation must have one coefficient = 1.";
  end if;
  C := PolMat(C, tr1);
  Cred := PF!C;
  assert Setseq(Monomials(Cred)) subset {PF.1^2*PF.i : i in [1..3]};
  fact2 := &+[PF.i*MonomialCoefficient(Cred, PF.1^2*PF.i) : i in [1..3]]; 
  a,b,c := Explode([MonomialCoefficient(fact2, PF.i) : i in [1..3]]);
  if b eq 0 and c eq 0 then
    return Mini300(C, N);
  elif IsUnit(b) then
    tr1 := Matrix(Rationals(),
                  [[1,0,0], [Integers() | -a/b, 1, -c/b], [0,0,1]]);
    return Mini210(PolMat(C, tr1), N);
  elif IsUnit(c) then
    tr1 := Matrix(Rationals(),
                  [[1,0,0], [0,1,0], [Integers() | -a/c, -b/c, 1]]);
    return Mini210(PolMat(C, tr1), N);
  else
    return Minimise(C, BreakUpWith(N, [Integers() | b,c]) : Check := false);
  end if;
end function;

function MinSing1(C, N)
  // here C = c x^2 z + F(x,y) mod N, c invertible
  vprintf Minimise, 2: " MinSing1(N = %o)\n", N;
  F := Integers(N);
  PF := PolynomialRing(F,3);
  Cred := PF!C;
  g := GCD([Integers() | MonomialCoefficient(Cred, PF.1^i*PF.2^(3-i))
                       : i in [0..3]]);
  if g eq 1 then
    P3 := Parent(C);
    Cnew := N*Evaluate(C, [P3.1, P3.2, 1/N*P3.3]);
    if forall{c : c in Coefficients(Cnew) | Denominator(c) eq 1} then
      tr1 := Matrix([[1,0,0], [0,0,-1], [0,1,0]]);
      return Mini2var(PolMat(Cnew, tr1), N);
    else
      printf "Minimisation not possible.\n";
      return C;
    end if;
  elif F!g eq 0 then
    tr1 := Matrix([[1,0,0], [0,0,-1], [0,1,0]]);
    return Mini210(PolMat(C, tr1), N);
  else
    return Minimise(C, BreakUpWith(N, [g]) : Check := false);
  end if;
end function;

function MinSing(C, N, sing)
  // here sing are coordinates of a singular point mod N
  vprintf Minimise, 2: " MinSing(N = %o, sing = %o)\n", N, sing;
  F := Integers(N);
  PF := PolynomialRing(F,3);
  // move sing to (0 : 0 : 1)
  if sing[1] eq 1 then
    tr1 := Matrix([[1,0,0],
                   [Integers() | sing[2],1,0],
                   [Integers() | sing[3],0,1]]);
    tr1 := tr1*Matrix([[0,0,-1], [0,1,0], [1,0,0]]);
  elif sing[2] eq 1 then
    tr1 := Matrix([[Integers() | 1,sing[1],0],
                   [0,1,0],
                   [Integers() | 0,sing[3],1]]);
    tr1 := tr1*Matrix([[1,0,0], [0,0,-1], [0,1,0]]);
  elif sing[3] eq 1 then
    tr1 := Matrix([[Integers() | 1,0,sing[1]],
                   [Integers() | 0,1,sing[2]],
                   [0,0,1]]);
  else
    error "One of the coordinates of sing must be 1.";
  end if;
  C := PolMat(C, tr1);
  Cred := PF!C;
  assert forall{m : m in {PF.3^3, PF.1*PF.3^2, PF.2*PF.3^2}
                  | MonomialCoefficient(Cred, m) eq 0};
  // now must have that degree-2 part in (x,y) is a square
  a, b, c := Explode([MonomialCoefficient(Cred, PF.1^(2-i)*PF.2^i*PF.3)
                       : i in [0..2]]);
  if b^2 ne 4*a*c then
    printf "Minimisation not possible.\n";
    return C;
  end if;
  if a eq 0 then
    if b ne 0 then
      return Minimise(C, BreakUpWith(N, [Integers() | b]) : Check := false);
    elif c eq 0 then // no F_2(x,y)*z
      return Mini2var(C, N);
    elif IsUnit(c) then
      tr1 := Matrix([[0,-1,0], [1,0,0], [0,0,1]]);
      return MinSing1(PolMat(C, tr1), N);
    else
      return Minimise(C, BreakUpWith(N, [Integers() | c]) : Check := false);
    end if;
  elif IsUnit(a) then
    if N eq 2 then
      tr1 := Matrix([[1,-Integers()!(c/a),0], [0,1,0], [0,0,1]]);
      return MinSing1(PolMat(C, tr1), N);
    elif IsDivisibleBy(N, 2) then
      return Minimise(C, BreakUpWith(N, [2]) : Check := false);
    else // N is odd
      tr1 := Matrix([[1,-Integers()!(b/(2*a)),0], [0,1,0], [0,0,1]]);
      return MinSing1(PolMat(C, tr1), N);
    end if;
  else
    return Minimise(C, BreakUpWith(N, [Integers() | a]) : Check := false);
  end if;
  return C;
end function;

intrinsic Minimise(C::RngMPolElt, N::RngIntElt : Check := true)
  -> RngMPolElt, Mtrx
{Minimise C by N.}
  vprintf Minimise, 1: "Minimise(N = %o),\n Curve = %o\n", N, C;
  N := Abs(N);
  P3 := Parent(C);
  require Rank(P3) eq 3: "Form must be ternary.";
  require forall{m : m in Monomials(C) | TotalDegree(m) eq 3}:
          "Polynomial must be a cubic form.";
  require BaseRing(P3) cmpeq Rationals():
          "Form must live in polynomial ring over the rationals.";
  require forall{c : c in Coefficients(C) | Denominator(c) eq 1}:
          "Form must have integral coefficients.";
  C := 1/GCD(ChangeUniverse(Coefficients(C), Integers()))*C;
  if Check then
    S,T := T3Invariants(C);
    S := Integers()!S;
    T := Integers()!T;
    require IsDivisibleBy(S, N^4) and IsDivisibleBy(T, N^6):
            "N^4, N^6 do not divide the invariants.";
  end if;
  Corig := C;
  if N eq 1 then return C; end if;
  F := Integers(N);
  PF := PolynomialRing(F, 3);
  Cred := PF!C;
  derivs := [Derivative(Cred, i) : i in [1..3]];
  if derivs[1] eq 0 then // no x
    return Mini2var(C, N); 
  elif derivs[2] eq 0 then // no y
    tr1 := Matrix(Rationals(), [[0,1,0], [-1,0,0], [0,0,1]]);
    return Mini2var(PolMat(C, tr1), N);
  elif derivs[3] eq 0 then // no z
    tr1 := Matrix(Rationals(), [[0,0,1], [0,1,0], [-1,0,0]]);
    return Mini2var(PolMat(C, tr1), N);
  end if;
  gl := [GCD(N, GCD(ChangeUniverse(Coefficients(Fd), Integers()))) 
          : Fd in derivs];
  if exists{g : g in gl | g ne 1} then
    return Minimise(Corig, BreakUpWith(N, gl) : Check := false);
  end if;
  // here, all the derivatives are of content prime to N
  mat := Matrix([[MonomialCoefficient(Fd, PF.1*PF.j) : j in [1..3]]
                   : Fd in derivs]);
  det := Determinant(mat);
  if det eq 0 then
    // distinguish further depending on rank of mat
    minors := Minors(mat, 2);
    if forall{m : m in minors | m eq 0} then
      // rank is at most 1
      if mat eq 0 then
        // rank is 0
        // derivatives do not involve x -- singularities along a line
        // derivatives must have common factor
        lcs := [MonomialCoefficient(Fd, PF.2^2) : Fd in derivs];
        _, pos := Position([IsUnit(l) : l in lcs], true);
        if forall{c : c in lcs | c eq 0} then
          // z is common factor, line is y = 0
          return MinLine(C, N, [F | 0, 1, 0]);
        elif pos ne 0 then
          q1 := derivs[pos];
          qs := [derivs[i] : i in [1..3] | i ne pos];
          q2 := qs[1] - MonomialCoefficient(qs[1], PF.2^2)*q1;
          q2 := &+[MonomialCoefficient(q2, PF.2*m)*m : m in [PF.2,PF.3]];
          q3 := qs[2] - MonomialCoefficient(qs[2], PF.2^2)*q1;
          q3 := &+[MonomialCoefficient(q3, PF.2*m)*m : m in [PF.2,PF.3]];
          lq2 := MonomialCoefficient(q2, PF.2);
          lq3 := MonomialCoefficient(q3, PF.2);
          assert lq2 ne 0 and lq3 ne 0;
          if IsUnit(lq2) then
            return MinLine(C, N, [0, 1, MonomialCoefficient(q2, PF.3)/lq2]);
          elif IsUnit(lq3) then
            return MinLine(C, N, [0, 1, MonomialCoefficient(q3, PF.3)/lq3]);
          else
            return Minimise(Corig, BreakUpWith(N, [Integers() | lq2, lq3])
                            : Check := false);
          end if;
        else
          return Minimise(Corig, BreakUpWith(N, ChangeUniverse(lcs, Integers()))
                          : Check := false);
        end if; // all lcs = 0
      elif exists{c : c in Eltseq(mat) | IsUnit(c)} then
        // rank is 1
        _, pos := Position([IsUnit(c) : c in Eltseq(mat)], true);
        ii := 1 + ((pos-1) div 3);
        jj := 1 + ((pos-1) mod 3);
        p1 := 1/mat[ii,jj]*derivs[ii];
        ps := [derivs[i] - MonomialCoefficient(derivs[i], PF.1*PF.jj)*p1
                : i in [1..3] | i ne ii];
        assert forall{p : p in ps | Seqset(Monomials(p))
                                     subset {PF.2^i*PF.3^(2-i) : i in [0..2]}};
        // two possibilities:
        // - either ps[1] = ps[2] = 0, and p1 is a square
        // - or ps[1] and ps[2] have a common root
        if ps[1] eq 0 and ps[2] eq 0 then
          if IsEven(N) and N ne 2 then
            return Minimise(Corig, BreakUpWith(N, [2]) : Check := false);
          end if;
          if N eq 2 then
            assert forall{m : m in [PF.1*PF.2, PF.2*PF.3, PF.3*PF.1] 
                            | MonomialCoefficient(p1, m) eq 0};
            return MinLine(C, N, 
                           [MonomialCoefficient(p1, PF.i^2) : i in [1..3]]);
          else // N is odd
            assert jj eq 1;
            // so p1 = x^2 + ...
            assert Coefficient(p1, PF.1, 1)^2 eq 4*Coefficient(p1, PF.1, 0);
            return MinLine(C, N,
                           [F | 1, MonomialCoefficient(p1, PF.1*PF.2)/2,
                                   MonomialCoefficient(p1, PF.1*PF.3)/2]);
          end if;
        end if;
        q2 := ps[1];
        q3 := ps[2];
        lq2 := MonomialCoefficient(q2, PF.2);
        lq3 := MonomialCoefficient(q3, PF.2);
        if lq2 eq 0 and lq3 eq 0 then
          // z is common factor
          p1 := Evaluate(p1, [PF.1, PF.2, 0]);
          if p1 eq 0 then
            return MinLine(C, N, [F | 0, 0, 1]);
          end if;
          g := GCD(N, GCD(ChangeUniverse(Coefficients(p1), Integers())));
          if g ne 1 then
            return Minimise(Corig, BreakUpWith(N, [g]) : Check := false);
          else
            // p1 must be a square
            lp1 := MonomialCoefficient(p1, PF.1^2);
            if IsEven(N) and N ne 2 then
              return Minimise(Corig, BreakUpWith(N, [2]) : Check := false);
            elif IsUnit(lp1) then
              if N eq 2 then
                assert MonomialCoefficient(p1, PF.1*PF.2) eq 0;
                return MinSing(C, N,
                               [1, MonomialCoefficient(p1, PF.2^2), 0]);
              else
                a, b, c := Explode([MonomialCoefficient(p1, PF.1^(2-i)*PF.2^i)
                                     : i in [0..2]]);
                assert b^2 eq 4*a*c;
                return MinSing(C, N, [1, -b/(2*a), 0]);
              end if;
            elif lp1 eq 0 then
              assert MonomialCoefficient(p1, PF.1*PF.2) eq 0;
              return MinSing(C, N, [F | 1, 0, 0]);
            else
              return Minimise(Corig, BreakUpWith(N, [Integers()!lp1])
                              : Check := false);
            end if;
          end if;
        elif IsUnit(lq2) then
          // ???
          printf "Not yet implemented (IsUnit(lq2)) -- sorry.\n";
          return C;
        elif IsUnit(lq3) then
          // ???
          printf "Not yet implemented (IsUnit(lq3)) -- sorry.\n";
          return C;
        else
          return Minimise(Corig, BreakUpWith(N, [Integers() | lq2, lq3])
                          : Check := false);
        end if;
      else
        // split N
        return Minimise(Corig, 
                        BreakUpWith(N, ChangeUniverse(Eltseq(mat), Integers()))
                        : Check := false);
      end if; // mat eq 0
    elif exists{m : m in minors | IsUnit(m)} then
      // rank is 2
      // ???
      printf "Not yet implemented (rank(mat) = 2) -- sorry.\n";
      return C;
    else
      // split N
      return Minimise(Corig, BreakUpWith(N, ChangeUniverse(minors, Integers()))
                      : Check := false);
    end if; // all minors = 0
  elif not IsUnit(det) then
    return Minimise(Corig, BreakUpWith(N, [Integers()!det]) : Check := false);
  end if; // det eq 0
  // here mat is invertible
  mati := mat^-1;
  pols := [&+[mati[i,j]*derivs[j] : j in [1..3]] : i in [1..3]];
  qpols := [pols[i] - PF.1*PF.i : i in [1..3]];
  assert forall{q : q in qpols | Derivative(q, 1) eq 0};
  p1 := PF.3*qpols[2] - PF.2*qpols[3];
  p2 := PF.2*PF.3*qpols[1] + qpols[2]*qpols[3];
  if p1 eq 0 and p2 eq 0 then
    // singular line given by x = -Q2(y,z)/y = -Q3(y,z)/z
    assert Evaluate(qpols[2], [0, 0, PF.3]) eq 0;
    Q2 := qpols[2];
    return MinLine(C, N, [1, MonomialCoefficient(Q2, PF.2^2),
                             MonomialCoefficient(Q2, PF.2*PF.3)]);
  end if;
  gl := [GCD(N, GCD(ChangeUniverse(Coefficients(Fd), Integers()))) 
          : Fd in [p1, p2]];
  if exists{g : g in gl | g ne 1} then
    return Minimise(Corig, BreakUpWith(N, gl) : Check := false);
  end if;
  // p1 and p2 are primitive. Find common root.
  lc1 := MonomialCoefficient(p1, PF.2^3);
  lc2 := MonomialCoefficient(p2, PF.2^4);
  if lc1 eq 0 and lc2 eq 0 then
    // singularity at (-Q3(0,1) : 0 : 1)
    return MinSing(C, N, [-Evaluate(qpols[3], [0,0,1]), 0, 1]);
  end if;
  if not IsUnit(lc1) or not IsUnit(lc2) then
      return Minimise(Corig, BreakUpWith(N, [Integers() | lc1,lc2])
                       : Check := false);
  end if; 
  p1 := 1/lc1*p1;
  p2 := 1/lc2*p2;
  ptemp := p2 - PF.2*p1;
  ptemp := ptemp - MonomialCoefficient(ptemp, PF.2^3*PF.3)*PF.2*p1;
  p3 := &+[mon*MonomialCoefficient(ptemp, PF.2^2*mon)
           : mon in [PF.2^i*PF.3^(2-i) : i in [0..2]]];
  lc3 := MonomialCoefficient(p3, PF.2^2);
  if lc3 eq 0 then
    // lc3 is linear, so must be the GCD
    lc3 := MonomialCoefficient(p3, PF.2*PF.3);
    assert lc3 ne 0;
    if not IsUnit(lc3) then
    return Minimise(Corig, BreakUpWith(N, [Integers() | lc3])
                       : Check := false);
    end if;
    zsing := -MonomialCoefficient(p3, PF.3^2)/lc3;
    assert Evaluate(p1, [0,1,zsing]) eq 0 and Evaluate(p2, [0,1,zsing]) eq 0;
    return MinSing(C, N, [-Evaluate(qpols[2], [0,1,zsing]), 1, zsing]);
  elif not IsUnit(lc3) then
    return Minimise(Corig, BreakUpWith(N, [Integers() | lc3])
                       : Check := false);
  else
    p3 := 1/lc3*p3;
    ptemp := p1 - PF.2*p3;
    ptemp := ptemp - MonomialCoefficient(ptemp, PF.2^2*PF.3)*PF.2*p3;
    p4 := &+[mon*MonomialCoefficient(ptemp, PF.2^2*mon) : mon in [PF.2,PF.3]];
    lc4 := MonomialCoefficient(p4,PF.2);
    assert lc4 ne 0;
    if not IsUnit(lc4) then
    return Minimise(Corig, BreakUpWith(N, [Integers() | lc4])
                       : Check := false);
    end if;
    zsing := -MonomialCoefficient(p4, PF.3);
    assert Evaluate(p1, [0,1,zsing]) eq 0 and Evaluate(p2, [0,1,zsing]) eq 0;
    return MinSing(C, N, [-Evaluate(qpols[2], [0,1,zsing]), 1, zsing]);
  end if; // lc3 eq 0
end intrinsic;
*/
