// Computing geometric invariants of genus one curves

// Author :   T.A.Fisher
// Version :  February 2008
// Requires Magma version 2.14.

MC := MonomialCoefficient;
MD := MonomialsOfDegree;
Coeff := Coefficient;

intrinsic RationalGCD(S::Setq[RngElt]) -> RngElt
{The greatest common divisor of the elements of the set/sequence S.}
   require #S gt 0 : "Illegal null set/sequence";
   B := Universe(S);
   require ISA(Type(B),FldRat) or ISA(Type(B),RngInt) :
     "Universe of set/sequence must be Z or Q.";  
   d := LCM([Denominator(Rationals()!x):x in S| x ne 0] cat [1]);
   return B!(GCD([IntegerRing()!(x*d): x in S])/d);
end intrinsic;

function PrimitiveModel(phi)
// The GCD of the coefficients of the model phi 
  B := CoefficientRing(phi);
  R := PolynomialRing(phi);
  n := Degree(phi);
  assert n in {3,4,5};
  if B cmpeq Rationals() then 
    coeffs := Eltseq(phi);
    dd := RationalGCD(coeffs);
    coeffs := [x/dd : x in coeffs];
return GenusOneModel(B,n,coeffs:PolyRing:=R);
  else 
    return phi;
  end if;
end function;
 
intrinsic InvariantDifferential(phi::ModelG1) -> AlgMatElt
{The invariant differential defined by the genus one model phi (of degree 3,4 or 5).}
  n := Degree(phi);
  require n in {3,4,5} : 
    "Input must be a genus one model of degree 3,4 or 5.";
  R := PolynomialRing(phi);
  B := CoefficientRing(phi);
  case n :
    when 3 :
      f := Equation(phi);
      Omega := ZeroMatrix(R,4,4);
      for pi in SymmetricGroup(3) do
        r,s,i := Explode([a^pi: a in [1..3]]); 
        Omega[r,s] := Sign(pi)*Derivative(f,i); 
      end for;
    when 4 :
      P := Equations(phi);
      Q := [-P[2],P[1]];
      Omega := ZeroMatrix(R,4,4);
      for pi in SymmetricGroup(4) do
        r,s,i,j := Explode([a^pi: a in [1..4]]); 
        Omega[r,s] +:= Sign(pi)*Derivative(P[1],i)*Derivative(P[2],j);
      end for;
    when 5 :
      P := Equations(phi);
      PHI := Matrix(phi);
      Jmat := Matrix(R,5,5,[Derivative(P[i],R.j): i,j in [1..5]]);
      M1 := [Matrix(B,5,5,[MC(PHI[i,j],R.k):i,j in [1..5]]): k in [1..5]];
      M2 := [Transpose(Jmat)*M*Jmat: M in M1];
      Omega := ZeroMatrix(R,5,5);
      for pi in SymmetricGroup(5) do
        r,s,i,j,k := Explode([a^pi: a in [1..5]]);
        if i lt j and j lt k then 
          Omega[r,s] := Sign(pi)*M2[j][i,k];
        end if;
      end for;
    end case;
  return Omega;
end intrinsic;

intrinsic InvariantDifferential(C::Crv) -> AlgMatElt
{An invariant differential on the genus one normal curve C, computed by linear algebra.}
  R := CoordinateRing(Ambient(C));
  K := BaseRing(R);
  n := Rank(R);  
  require n ge 4 : "C must have degree at least 4.";
  qq := Equations(C);
  mons := MD(R,2);
  V := VectorSpace(K,#mons);
  U := sub<V|[V![MC(q,mon):mon in mons] : q in qq]>;
  W,quomap := quo<V|U>;
  quads := [&+[x[i]*mons[i] : i in [1..#mons]] 
    where x is Eltseq(W.j @@ quomap) : j in [1..Dimension(W)]];
  bb := [Matrix(R,n,n,[<i,j,1>,<j,i,-1>]): i,j in [1..n] | i lt j];
  BB := [q*b: q in quads,b in bb];
  mons := MD(R,3);
  V:= VectorSpace(K,#mons);
  U := sub<V|[V![MC(f*q,mon):mon in mons] : q in qq,f in MD(R,1)]>;
  W,quomap := quo<V|U>;
  function polytoseq(f)
    vec := V![MC(f,mon):mon in mons];
    return Eltseq(quomap(vec));
  end function;
  Jmat := Matrix(#qq,n,[Derivative(q,R.i): i in [1..n],q in qq]);
  mat := Matrix(K,#BB,n*#qq*Dimension(W),
     [&cat[polytoseq(f): f in Eltseq(Jmat*B)]: B in BB]);
  km := KernelMatrix(mat);
  assert Nrows(km) eq 1;
  return &+[km[1,i]*BB[i]: i in [1..#BB]];
end intrinsic;

function ordinalstring(r) 
  str := "th";
  if r mod 10 eq 1 and r mod 100 ne 11 then str:= "st"; end if;
  if r mod 10 eq 2 and r mod 100 ne 12 then str:= "nd"; end if;
  if r mod 10 eq 3 and r mod 100 ne 13 then str:= "rd"; end if;
  return IntegerToString(r) cat str;
end function;

intrinsic WronskianPFunction(C::Crv,Omega::AlgMatElt)
     -> FldFunFracSchElt, FldFunFracSchElt
{The Weierstrass P function and its derivative, computed using Anderson\'s method.}
  Rx := CoordinateRing(Ambient(C));
  n := Rank(Rx);
  K := FunctionField(C);
  Ry := PreimageRing(Parent(Numerator(K.1*K.2)));
  assert forall{i : i in [1..n-1] | K!(Rx.i/Rx.n) eq K!Ry.i};
  dd := [K|(Omega[i,n]/Rx.n^2) : i in [1..n-1]];
  function deriv1(x)
    y := Ry!x;
    return &+[K!Derivative(y,Ry.i)*dd[i]: i in [1..n-1]];
  end function;
  function deriv(x)
    a := Ry!Numerator(x);
    b := Ry!Denominator(x);
    return K!(((K!b)*deriv1(a)-(K!a)*deriv1(b))/(K!b)^2);
  end function;
  b := [K|Rx.i/Rx.n : i in [1..n]];
  bb := [b];
  for r in [1..n+1] do
    printf "  Computing %o derivatives\n",ordinalstring(r);
    time bnew := [deriv(x)/r: x in b];
    b := bnew;
    bb cat:= [b];
  end for;
  W := Determinant(Matrix(n,n,[bb[i]: i in [1..n]]));
  W1 := Determinant(Matrix(n,n,[bb[i]: i in [1..n-1] cat [n+1]]));
  W2 := Determinant(Matrix(n,n,[bb[i]: i in [1..n-1] cat [n+2]]));
  W11 := Determinant(Matrix(n,n,[bb[i]: i in [1..n-2] cat [n,n+1]]));
  printf "  Computing P\n";
  time P := (W1/W)^2-(W2/W)-(W11/W);
  printf "  Computing P'\n";
  time PP := deriv(P)/n;
  return P,PP;
end intrinsic;

function NextElement(m,seq,y) 
  F1 := Domain(m);
  A := Codomain(m);
  P<x> := PolynomialRing(F1);
  poly := P!MinimalPolynomial(y);
  ff := [f[1]: f in Factorization(poly)];
  assert exists(f){f : f in ff |
    &+[m(Coeff(f,i))*y^i: i in [0..Degree(f)]] eq 0};
  F2 := ext<F1|f:DoLinearExtension>;
  F3 := AbsoluteField(F2);
  seq := [F3!x: x in seq] cat [F3!F2.1];
  f := Eltseq(F2!F3.1);
  elt := &+[m(f[i])*y^(i-1): i in [1..#f]];
  m := hom<F3->A|elt>;
  return m,seq;
end function;
  
function PassToField(seq)
// Computes the field generated by a sequence of elements 
// of an algebraically closed field, and returns the sequence 
// as elements of the new field.
  A := Universe(seq);
  m := hom<Rationals()->A|>;
  ans := [];
  for y in seq do
    m,ans := NextElement(m,ans,y);
  end for;
  return ans;
end function;

function IntersectWithHyperplane(C)
// Finds a point on C (over a suitable field extension)
// by intersecting with a hyperplane. 
  P := Ambient(C);
  R := CoordinateRing(P);
  X := Scheme(P,Equations(C) cat [R.1]);
  pt := PointsOverSplittingField(X)[1];
  pt := PassToField(Eltseq(pt));
  K<u> := Universe(pt);
  return C(K)!pt;
end function;

function MovePoint(C,Omega,P)
// Given a point P on the curve C, changes co-ordinates
// so that P = (0:0: ... :0:1). The matrix Omega, describing
// an invariant differential on C is also updated.
  K := Ring(Parent(P));
  C := BaseChange(C,K);
  R := CoordinateRing(Ambient(C));
  n := Rank(R);
  pt := Eltseq(P);
  assert exists(i){i : i in [n..1 by -1]|pt[i] ne 0};
  I := IdentityMatrix(K,n);
  II := [Eltseq(I[j]) : j in [1..n] | j ne i];
  M := Transpose(Matrix(K,n,n,II cat [pt]));
  eqns := Equations(C);
  eqns := [f^M : f in eqns];
  Omega := Matrix(R,n,n,[(R!Omega[i,j])^M: i,j in [1..n]]);
  M1 := M^(-1);
  Omega := M1*Omega*Transpose(M1);
  assert forall{f : f in eqns | MC(f,R.n^2) eq 0};
  return Scheme(Proj(R),eqns),Omega;
// Using "Curve" instead of "Scheme" prompts Magma
// to do some unnecessary and time consuming checks.
end function;

function ProjectFromPoint(C,Omega)
// Let C be a genus one normal curve of degree n>=4, that passes
// through the point P = (0:0: ... :0:1). Projecting away from P
// we obtain a genus one normal curve of degree n-1. The matrix Omega, 
// describing an invariant differential on C is updated, and 
// the image of P is also returned.
  eqns := Equations(C);
  R := Universe(eqns);
  K := BaseRing(R);
  n := Rank(R);
  assert #eqns eq n*(n-3)/2;
  assert forall{f : f in eqns | TotalDegree(f) eq 2};
  assert forall{f : f in eqns | MC(f,R.n^2) eq 0};
  mat := Matrix(K,#eqns,n,[MC(f,R.i*R.n):i in [1..n],f in eqns]);
  mat1 := Matrix(K,n-1,#eqns,[mat[i,j]:i in [1..#eqns],j in [1..n-1]]);
  function RemoveLast(f)
    vec := Vector(K,n,[MC(f,R.i*R.n): i in [1..n]]);
    sol := Solution(mat,vec);
    g := &+[sol[i]*eqns[i]: i in [1..#eqns]];
    return f-g;
  end function;
  for i,j in [1..n-1] do
    Omega[i,j] := RemoveLast(Omega[i,j]);
  end for;
  m := Integers()!((n-1)*(n-4)/2);
  km := RMatrixSpace(R,#eqns,m)!Transpose(KernelMatrix(mat));
  if n eq 4 then 
    ll := [&+[MC(f,R.i*R.4)*R.i: i in [1..3]] : f in eqns];
    eqns := [ll[1]*eqns[2]-ll[2]*eqns[1]];
  else 
    eqns := Eltseq(Vector(eqns)*km);
  end if;
  S := PolynomialRing(K,n-1);
  xx := [S.i: i in [1..n-1]] cat [0];
  eqns := [S|Evaluate(f,xx):f in eqns];
//  C := Curve(Proj(S),eqns);
  C := Scheme(Proj(S),eqns);
  Omega := Matrix(S,n-1,n-1,[Evaluate(Omega[i,j],xx):i,j in [1..n-1]]);
  km := KernelMatrix(mat1);
  P := C(K)!Eltseq(km[1]);
  return C,Omega,P;
end function;

function LinearDependence(seq)
// Finds a linear dependence between a sequence of 
// elements of a function field. 
  K := Universe(seq);
  n := #seq;
  for i in [1..n] do
    d := K!Denominator(seq[i]);
    seq := [d*x: x in seq];
  end for;
  A := Parent(Numerator(K.1));
  seq := ChangeUniverse(seq,A);
  R := PreimageRing(A);
  polys := ChangeUniverse(seq,R);
  mons := &join[SequenceToSet(Monomials(f)): f in polys];
  mons := SetToSequence(mons);
  mat := Matrix(#polys,#mons,[MC(f,mon): mon in mons,f in polys]);
  km := KernelMatrix(mat);
  return km[1];
end function;

function SchemeToGenusOneModel(C) 
  K := BaseRing(C);
  assert IsField(K);
  P := AmbientSpace(C);
  R := CoordinateRing(P);
  n := Rank(R);
  assert n in {3,4,5};
  eqns := MinimalBasis(Ideal(C));
  return GenusOneModel(eqns);
end function;

intrinsic GeometricInvariants(C::Sch,Omega::AlgMatElt
  : Method:="Invariants",nproj:=1,P:=0) -> SeqEnum 
{Computes the geometric invariants c4 and c6 of the pair (C,omega), 
where C is a genus one normal curve, and omega is an invariant
differential on C specified by a matrix of quadrics Omega by the rule 
omega = x_j^2 d(x_i/x_j) / Omega_\{ij\} for all i<>j. The methods available
are "Invariants", "Projection" and "Wronskians". The parameter nproj
specifies the number of projections used. A point P may be specified for
the projection method.}
  require Method in {"Invariants","Wronskians","Projection"} : 
  "Method should be either \"Invariants\", \"Wronskians\" or \"Projection\"";
  n := Dimension(Ambient(C))+1;
  if n eq 3 and Method eq "Projection" then
    Method := "Invariants"; 
  end if;
  printf "n = %o. Using the %o method\n",n,Method;
  case Method :
    when "Projection" :
      if P cmpeq 0 then 
        P := IntersectWithHyperplane(C);
      end if;
      C1,Omega1 := MovePoint(C,Omega,P);
      C2,Omega2,P := ProjectFromPoint(C1,Omega1);
      mm := nproj eq 1 select "Invariants" else "Projection";
      cc := GeometricInvariants(C2,Omega2:Method:=mm,nproj:=nproj-1,P:=P);
    when "Invariants" :
      require n le 5 : "Curve must have degree at most 5";
      phi := SchemeToGenusOneModel(C);
      phi := PrimitiveModel(phi);
      Omega1 := InvariantDifferential(phi);
      nf1 := NormalForm(Omega[1,2],Ideal(C));
      nf2 := NormalForm(Omega1[1,2],Ideal(C));
      u := nf1/nf2;
      c4,c6 := Explode(cInvariants(phi));
      cc := [u^4*c4,u^6*c6];
    when "Wronskians" :
      P,PP := WronskianPFunction(C,Omega);
      km := LinearDependence([PP^2-4*P^3,P,1]);
      cc := [12*km[2]/km[1],216*km[3]/km[1]];
  end case;
  return ChangeUniverse(cc,BaseRing(C));
end intrinsic;