
 //////////////////////////////////////////////////////
 //                                                  //
 //   omega.m                                        //
 //                                                  //
 //   Computations accompanying the article          //
 //       "A formula for the Jacobian of             //
 //    a genus one curve of arbitrary degree"        //
 //                                                  //
 //   Author: T.A. Fisher                            // 
 //   Version: June 2018                             //
 //                                                  //
 //////////////////////////////////////////////////////

MD := MonomialsOfDegree;
MC := MonomialCoefficient;
Deriv := Derivative;

load "~/public_html/papers/kleinmatrix.m";
// available from 
// http://www.dpmms.cam.ac.uk/~taf1000/papers/kleinmatrix.m

// The invariants c4(Omega) and c6(Omega) as defined in Section 1

function OmegaMatrixInvariants(Omega)
  R := BaseRing(Omega);
  n := Nrows(Omega);
  assert Transpose(Omega) eq - Omega;
  assert Rank(R) eq n;
  for i in [1..n] do
    assert &+[Deriv(Omega[i,j],j): j in [1..n]] eq 0;
  end for;
  M := Matrix(R,n,n,[&+[Deriv(Omega[i,r],s)*Deriv(Omega[j,s],r)
          :r,s in [1..n]]:i,j in [1..n]]);
  N := [Matrix(R,n,n,[&+[Deriv(M[i,j],r)*Omega[r,k]
          :r in [1..n]]:j,k in [1..n]]): i in [1..n]];
  inv4 := &+[Deriv(Deriv(M[i,j],r),s)*Deriv(Deriv(M[r,s],i),j)
          :i,j,r,s in [1..n]];
  inv6 := &+[Deriv(Deriv(Deriv(N[i,j,k],r),s),t)
            *Deriv(Deriv(Deriv(N[r,s,t],i),j),k)
	  :i,j,k,r,s,t in [1..n]];
  c4 := 3*(n-2)^2/(2^4*n*Binomial(n+3,5))*MC(inv4,1);
  c6 := -(n-2)^3/(2^6*n*Binomial(n+5,7))*MC(inv6,1);
  return [BaseRing(R)|c4,c6];
end function;

// The Omega matrix associated to a Weierstrass equation as defined
// in Theorem 7.1.

function OmegaFromWeierstrass(K,n,ainv:R := PolynomialRing(K,n))
  a1,a2,a3,a4,a6 := Explode(ainv);
  assert Rank(R) eq n;
  xx := func<i|i eq 0 select R.1 else i in [2..n] select R.i else 0>;
  aa := [a1,a2,a3,a4,0,a6];
  xxdot := func<m|(m/2)*(2*xx(m+1) + a1*xx(m) + a3*xx(m-2))
         + (IsOdd(m) select 1 else 0)
          *&+[(-1)^i*(m-(i/2))*aa[i]*xx(m+1-i) : i in [1..6]]>;
  xxbar := func<m|(1/2)*(2*xx(m+1) + a1*xx(m) + a3*xx(m-2))
         + (IsOdd(m) select 1 else 0)
          *&+[(-1)^i*aa[i]*xx(m+1-i) : i in [1..6]]>;
  A := Matrix(R,n,n,[xx(r)*xxdot(s) - xxdot(r)*xx(s)
                             : r,s in [0] cat [2..n]]);
  B := Matrix(R,n,n,[&+[Sign(k+1/2)
          *(xx(r + 2*k)*xxbar(s - 2*k) - xx(s + 2*k)*xxbar(r - 2*k))
        : k in [-100..100] ] : r,s in [0] cat [2..n]]);
  return n*B - 2*A,A,B;
end function;

// The degenerate Omega matrix in Lemma 8.4

function OmegaForNeronPolygon(K,n:R := PolynomialRing(K,n))
  return Matrix(R,n,n,[(Sign(j-i)*n - 2*(j-i))*R.i*R.j: i,j in [1..n]]);
end function;  

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

print "";
print "Section 3 : Curves of small degree";

printf "\nCase n = 3 : Plane cubics\n";
model := GenericModel(3);
F := Equation(model);
F1,F2,F3 := Explode([Derivative(F,i): i in [1..3]]);
Omega := Matrix(3,3,[0,F3,-F2,-F3,0,F1,F2,-F1,0]);
assert OmegaMatrixInvariants(Omega) eq cInvariants(model);

printf "\nCase n = 4 : Quadric intersections\n";
model := RandomModel(4);
// model := GenericModel(4); // takes a few minutes
F1,F2 := Explode(Equations(model));
Omega := ZeroMatrix(PolynomialRing(model),4,4);
for sigma in AlternatingGroup(4) do
  i,j,k,l := Explode([x^sigma : x in [1..4]]);
  Omega[i,j] := Deriv(F1,k)*Deriv(F2,l) - Deriv(F1,l)*Deriv(F2,k);
end for;
assert OmegaMatrixInvariants(Omega) eq cInvariants(model);

printf "\nCase n = 5 : Pfaffian models\n";
model := RandomModel(5);
Phi := Matrix(model);
P := Vector(Equations(model)); // Pfaffians
assert P*Phi eq 0;
assert Transpose(Phi) eq -Phi;
Omega := ZeroMatrix(PolynomialRing(model),5,5);
for sigma in AlternatingGroup(5) do
  i,j,k,l,m := Explode([x^sigma : x in [1..5]]);
  Pk := Matrix(1,5,[Deriv(f,k): f in Eltseq(P)]);
  Phil := Matrix(5,5,[Deriv(Phi[r,s],l): r,s in [1..5]]);
  Pm := Matrix(5,1,[Deriv(f,m): f in Eltseq(P)]);
  Omega[i,j] +:= (1/3)*(Pk*Phil*Pm)[1,1];
end for;
assert OmegaMatrixInvariants(Omega) eq cInvariants(model);

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

print "";
print "Direct verification (for some small values of n) that the";
print "invariants are correctly scaled";

printf "\nStarting with a general Weierstrass equation\nn =";
k<a1,a2,a3,a4,a6> := FunctionField(Rationals(),5);
E := EllipticCurve([a1,a2,a3,a4,a6]);
for n in [3..6] do
  printf " %o,",n;
  R<[x]> := PolynomialRing(k,n);
  Omega := 1/(n-2)*OmegaFromWeierstrass(k,n,aInvariants(E):R:=R);
  assert OmegaMatrixInvariants(Omega) eq cInvariants(E);  
end for;

printf "\n\nStarting with a shorter Weierstrass equation\nn =";
k<a,b> := FunctionField(Rationals(),2);
E := EllipticCurve([a,b]);
for n in [3..8] do
  printf " %o,",n;
  R<[x]> := PolynomialRing(k,n);
  Omega := 1/(n-2)*OmegaFromWeierstrass(k,n,aInvariants(E):R:=R);
  assert OmegaMatrixInvariants(Omega) eq cInvariants(E);  
end for;
print "";

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

print "";
print "Section 8 : Proof of Theorem 1.2";

printf "\nChecking the statement of Lemma 8.4\nn =";
for n in [3..8] do
  printf " %o,",n;
  Lambda := OmegaForNeronPolygon(Rationals(),n);
  assert OmegaMatrixInvariants(Lambda) eq [(n-2)^4,-(n-2)^6];
end for;  

printf "\n\nChecking the proof of Lemma 8.4\nn =";
for n in [3..15] do
  printf " %o,",n;
  la := func<r,s|Sign(s - r)*n - 2*(s - r)>;
  mu := func<i,j|2*(n - 2)*((i - j)^2 - n*Abs(i - j)) + n*(n - 1)*(n - 2)/3>;
  nu := func<i,j,k|mu(i,j)*(la(i,k) + la(j,k))>;
  LHS := 4*&+[mu(i,j)^2 : i,j in [1..n] | i le j];
  RHS := (16/3)*n*(n - 2)^2*Binomial(n + 3,5);
  assert LHS eq RHS;
  LHS := 4*&+[(nu(i,j,k) + nu(j,k,i) + nu(k,i,j))^2
                      : i,j,k in [1..n] | i le j and j le k];
  MID := 64*(n-2)^2*&+[(i - 2*j + k)^2*(n + i + j - 2*k)^2*(n + 2*i - j - k)^2
              : i,j,k in [1..n] | i le j and j le k];
  RHS := 64*n*(n-2)^3*Binomial(n+5,7);
  assert LHS eq MID and MID eq RHS;
end for;
print "";

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

print "";
print "Direct verification (for some small values of n) that the matrix"; 
print "in Theorem 7.1 satisfies the conclusions of Lemma 4.5";
print "";

k<a,b> := FunctionField(Rationals(),2);
E := EllipticCurve([a,b]);
K<x,y> := FunctionField(E);
xx := [1,x,y] cat (&cat[[x^(m+1),x^m*y] : m in [1..10]]);

for n in [4..8] do

  print "n =",n;
  R<[x]> := PolynomialRing(k,n);
  Omega := OmegaFromWeierstrass(k,n,aInvariants(E):R:=R);

  psi := map<E->Proj(R)|[xx[i]: i in [1..n]]>;
  C := Image(psi);
  quads := MinimalBasis(Ideal(C));
  assert #quads eq n*(n-3)/2;

  r := Floor((n-2)/2);
  ff := quads;     
  for rho in [1..r] do
    assert IsHomogeneousBasis(ff);
    assert forall{f : f in ff | TotalDegree(f) eq rho + 1};
    assert #ff eq beta(rho+1,n);  
    printf "I(Sec^%o C) generated by %o %o-ics\n",rho,#ff,rho+1;
      if rho eq r then break; end if;
      if n*beta(rho+1,n) gt Binomial(n+rho+1,rho+2) then 
        ff := NextSecantVariety1(ff);
      else
        ff := NextSecantVariety2(ff);
      end if;
  end for;
  
  if IsOdd(n) then
    assert #ff eq n;
    pf := PfSubmax(Omega);
    F := (1/n)*&+[R.i*pf[i] : i in [1..n]];
    for i in [1..n] do
      assert Derivative(F,i) eq pf[i];
    end for;
    mons := MD(R,n-1);
    forms := [ff[i]*ff[j] : i,j in [1..n] | i le j];
    V := VectorSpace(k,#mons);
    U1 := sub<V|[V![MC(f,mon) : mon in mons] : f in forms]>;
    U2 := sub<V|[V![MC(f,mon) : mon in mons] : f in pf]>;
    assert U2 subset U1; // see Proposition 9.3
    ff := [F];
  end if;

// Checking the conclusions of Lemma 4.5
  assert #ff in [1,2];
  for f in ff do
    for k in [1..n] do
      assert &+[Derivative(f,j)*Omega[j,k]: j in [1..n]] eq 0;
    end for;
  end for;
 
end for;

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

print "";
print "Checking the last part of the proof of Theorem 7.1";

k<a1,a2,a3,a4,a6> := FunctionField(Rationals(),5);
E := EllipticCurve([a1,a2,a3,a4,a6]);
n := 6;
R<[x]> := PolynomialRing(k,n);
_,_,B := OmegaFromWeierstrass(k,n,aInvariants(E):R:=R);

KE1<xP,yP> := FunctionField(E);
xxP := [1,xP,yP] cat (&cat[[xP^(m+1),xP^m*yP] : m in [1..10]]);
xxPn := [xxP[i] : i in [1..n]];
E2 := BaseChange(E,KE1);
KE2<xQ,yQ> := FunctionField(E2);
xxQ := [1,xQ,yQ] cat (&cat[[xQ^(m+1),xQ^m*yQ] : m in [1..10]]);

for r,s in [1..n-1] do
  lambda := (yP + yQ + a1*xQ + a3)/(xP - xQ);
  BB := lambda*(xxQ[r]*xxP[s] - xxP[r]*xxQ[s]);
  assert Evaluate(BB,E2![xP,yP]) eq Evaluate(B[r,s],xxPn);
end for;
print "";

x0,x2,x3,x4,x5 := Explode(x);
assert B[1,3] eq 2*x0*x4 + x2^2 - a1*x0*x3 + 2*a2*x0*x2 + a4*x0^2;
assert B[2,3] eq 2*x2*x4 - x3^2 - a1*x2*x3 + a2*x2^2 - a6*x0^2;

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

/*

Magma V2.23-8     Fri Jun 15 2018 16:48:11 on blackcap [Seed = 1879635127]
Type ? for help.  Type <Ctrl>-D to quit.
> load "omega.m";
Loading "omega.m"
Loading "/homes/taf1000/public_html/papers/kleinmatrix.m"

Section 3 : Curves of small degree

Case n = 3 : Plane cubics

Case n = 4 : Quadric intersections

Case n = 5 : Pfaffian models

Direct verification (for some small values of n) that the
invariants are correctly scaled

Starting with a general Weierstrass equation
n = 3, 4, 5, 6,

Starting with a shorter Weierstrass equation
n = 3, 4, 5, 6, 7, 8,

Section 8 : Proof of Theorem 1.2

Checking the statement of Lemma 8.4
n = 3, 4, 5, 6, 7, 8,

Checking the proof of Lemma 8.4
n = 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,

Direct verification (for some small values of n) that the matrix
in Theorem 7.1 satisfies the conclusions of Lemma 4.5

n = 4
I(Sec^1 C) generated by 2 2-ics
n = 5
I(Sec^1 C) generated by 5 2-ics
n = 6
I(Sec^1 C) generated by 9 2-ics
I(Sec^2 C) generated by 2 3-ics
n = 7
I(Sec^1 C) generated by 14 2-ics
I(Sec^2 C) generated by 7 3-ics
n = 8
I(Sec^1 C) generated by 20 2-ics
I(Sec^2 C) generated by 16 3-ics
I(Sec^3 C) generated by 2 4-ics

Checking the last part of the proof of Theorem 7.1

> 

*/
