
////////////////////////////////////////////////////////////////////
//                                                                //
// This Magma file (version 2.21-1) accompanies the paper         //
//                                                                //
// Title : Visualizing elements of order 7 in the                 //
//         Tate-Shafarevich group of an elliptic curve            //
// Author : Tom Fisher                                            //
// Date : August 2016 (updated September 2017)                    //
//                                                                //
////////////////////////////////////////////////////////////////////

CD := CremonaDatabase();
CR := func<E|Conductor(E) lt LargestConductor(CD) select 
    CremonaReference(E) else IntegerToString(Conductor(E)) cat "*">;
MyRowReduce := func<mat|LLL(Matrix(Basis(PureLattice(Lattice(mat)))))>;
PZ := PolynomialRing(Integers());
MC := MonomialCoefficient;
MD := MonomialsOfDegree;
Deriv := Derivative;

load "vis7data.m";

assert #shalist eq 222;
assert #table1 eq 79;
assert #FF eq 5;
assert #table4 eq 35;
assert #conductors eq 35;
assert #MWgens eq 35;
assert #cgdata eq 35;
assert #kumsrf_pts eq 35;

print "---------------------------------------------";
print "Section 3. Visibility in abelian surfaces";

// Theorem 3.1

function X7Twists(a,b)
  K := Universe([a,b]);
  R<x,y,z> := PolynomialRing(K,3);
// Halberstadt, Kraus : Equation for X_E(7)
    F := a*x^4 + 7*b*x^3*z + 3*x^2*y^2 - 3*a^2*x^2*z^2 - 6*b*x*y*z^2 
       - 5*a*b*x*z^3 + 2*y^3*z + 3*a*y^2*z^2 + 2*a^2*y*z^3 - 4*b^2*z^4;
// Poonen, Schaefer, Stoll : Equation for X_E^-(7)
    G := -a^2*x^4 + 2*a*b*x^3*y - 12*b*x^3*z + (-6*a^3 - 36*b^2)*x^2*y^2 
      + 6*a*x^2*z^2 + 2*a^2*b*x*y^3 - 12*a*b*x*y^2*z + 18*b*x*y*z^2 
      + (3*a^4 + 19*a*b^2)*y^4 + (-8*a^3 - 42*b^2)*y^3*z 
      + 6*a^2*y^2*z^2 - 8*a*y*z^3 + 3*z^4;
  return F,G;
end function;

// Example 3.2

a := -4062871;
b := -3152083138;
E := EllipticCurve("3364c");
assert IsIsomorphic(E,EllipticCurve([a,b]));
F,G := X7Twists(a,b);

R<x,y,z> := PolynomialRing(Rationals(),3);
M1 := Matrix(Rationals(),3,3,
  [ 2320, 0, -3509, -2716430, 1682, 4042687, -2, 0, 3 ]);
M2 := Matrix(Rationals(),3,3,
  [ 116, 1160, -145, 0, -1, 0, 134560, 2703815, -169041 ]);
subst1 := [&+[M1[i,j]*R.j : j in [1..3]]: i in [1..3]];
subst2 := [&+[M2[i,j]*R.j : j in [1..3]]: i in [1..3]];
F := (-1/9517173136)*Evaluate(F,subst1);
G := (-1/32015770429504)*Evaluate(G,subst2);

assert F eq x^4 + 2*x^3*y + x^3*z + 3*x^2*y^2 + 3*x^2*y*z + 3*x^2*z^2
  + 2*x*y^3 - 3*x*y^2*z + 6*x*y*z^2 + 12*x*z^3 - 3*y^3*z - 3*y^2*z^2
  + 6*y*z^3 - 2*z^4;

C := Curve(Proj(R),F);
pts := PointSearch(C,10^3);
assert Set(pts) eq {C|[0,1,0],[1,-1,0]};
C := Curve(Proj(R),G);
pts := PointSearch(C,10^3);
assert #pts eq 0;

// Proof of Theorem 3.3

aa := func<E|Vector(GF(7),TracesOfFrobenius(E,10^3))>;
nfail := func<E,F|#[x : x in Eltseq(aa(E)-aa(F)) | x ne 0]>;

goodtally := 0;
multtally := 0;
potgoodlist := [];
newthmlist := [];

P<X> := PolynomialRing(Rationals());
K := NumberField(X^4-7);
OK := RingOfIntegers(K);
pp7 := Factorization(7*OK)[1,1];

for EF in table1 do
  E := EllipticCurve(EF[1]);
  if exists(i){i : i in [1..#FF] | EF[2] eq Sprint(FF[i][1]) cat "*"} then
    F := EllipticCurve(FF[i][2]);
    assert Conductor(F) eq FF[i][1];
  else
    F := EllipticCurve(EF[2]);
  end if;
  assert CR(E) in {EF[1],EF[1] cat "1"};
  assert CR(F) in {EF[2],EF[2] cat "1"};
  assert nfail(E,F) le 1;
  assert EF[1] in shalist;
  if CR(E) ne "88450f1" then 
    bd1,bd2 := RankBounds(F); assert [bd1,bd2] eq [2,2];
  else 
    P1 := F![ 404428, -1584624508 ];
    P2 := F![ -57522281013974704137006860452533791294179689445988/290\
      85417710961952827793222434722107359927969, -2421774220485513671\
      35768271470188936406387471024934845450306478147020662172/156860\
      269909577671322999120267652141626078645010309044603708349297 ];
    assert Regulator([P1,P2]) gt 10^(-5);
    assert #TwoDescent(F) eq 3;
  end if;
  assert forall{ c : c in TamagawaNumbers(E) | c mod 7 ne 0};
  assert forall{ c : c in TamagawaNumbers(F) | c mod 7 ne 0};
  if {ReductionType(E,7),ReductionType(F,7)} eq {"Good"} then
    goodtally +:= 1;
  else
    multflag := false;
    for m in [1,-1,7,-7] do
      E1 := MinimalModel(QuadraticTwist(E,m));
      F1 := MinimalModel(QuadraticTwist(F,m));
      if {ReductionType(E1,7),ReductionType(F1,7)}
                             eq {"Split multiplicative"} then
        multflag := true;
      end if;
    end for;
    if multflag then 
      multtally +:= 1;
    elif ReductionType(E,7) eq "Nonsplit multiplicative" 
          and ReductionType(F,7) eq "Good" then
        Append(~newthmlist,EF[1]);
    else
      EK := MinimalModel(BaseChange(E,K),pp7);
      FK := MinimalModel(BaseChange(F,K),pp7);
      if Valuation(Discriminant(EK),pp7) eq 0 and 
           Valuation(Discriminant(FK),pp7) eq 0 then
        Append(~potgoodlist,EF[1]);
      else 
        assert EF[1] eq "84672kg";
      end if;
    end if;  
  end if;    
end for;

assert goodtally eq 50;
assert multtally eq 20;
assert potgoodlist eq [ "27930bj", "31311h", "71148bh", "84966ea" ];
assert newthmlist eq [ "46704k", "73416g", "75712bz", "92862j" ];

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

print "---------------------------------------------";
print "Section 5. Visibility in abelian threefolds";

/*
print "Computing modular forms";

primelimit := 100;
pp := [p : p in [2..primelimit] | IsPrime(p)];
K<a> := Qrt2;
OK := RingOfIntegers(K);
p7 := ideal<OK|3 + a>;

for i in [1..#table2] do
  dd := table2[i];
  print dd[1];
  E := EllipticCurve(dd[1]);
  apE := TracesOfFrobenius(E,primelimit);
  apf := dd[2];
  assert forall{i : i in [1..#pp] | (apE[i]-apf[i]) in p7};
  assert exists(j0){j : j in [1..#pp] | apf[j] notin Rationals()};
  assert exists(j1){j : j in [j0+1..#pp] | apf[j] notin Rationals()};
  fdata := [<pp[j],MinimalPolynomial(apf[j])>: j in [j0,j1]];
  N := Conductor(E);
  time M := ModularSymbols(N,2,1);     
  time M := CuspidalSubspace(M);       
  time M1 := Kernel(fdata,M); 
  assert Dimension(M1) eq 2;
  r,s := Explode(Eltseq(apf[j0]));
  time Tp := HeckeOperator(M1,pp[j0]);
  a0 := (Tp - r)/s;
  assert a0^2 eq 2;
  mat := Matrix(2,4,[Eltseq(a0^i): i in [0,1]]);
  print "Checking Hecke eigenvalues";
  for j in [1..#pp] do
    printf "%o ",pp[j];
    Tp := HeckeOperator(M1,pp[j]);
    soln := Solution(mat,Vector(4,Eltseq(Tp)));
    r,s := Explode(Eltseq(soln));
    assert Tp eq r + s*a0;
    assert apf[j] eq r + s*a;
  end for;
  print "";
end for;
*/

// Theorem 5.1 (Bending)

PQ<x> := PolynomialRing(Rationals());
function Genus2Curve(A,P,Q,delta)
  R := 4*P;
  B := (Q*(P*A-Q)+4*P^2+1)/P^2;
  C := 4*(P*A-Q)/P;
  f := x^6 + A*x^5 + (A^2*P^3 - 3*A*P^2*Q + A*P*Q + 4*P^3 
       + 4*P^2 + 2*P*Q^2 - 2*P - Q^2 + 1)/P^2*x^4 
       + (A^2*P^2*Q + 4*A*P^2 - 3*A*P*Q^2 + A*P + 4*P^2*Q - 4*P*Q 
       + 2*Q^3 - 2*Q)/P^2*x^3 + (4*A^2*P^3 - 12*A*P^2*Q + A*P*Q 
       + 16*P^3 - 8*P^2 + 8*P*Q^2 + 4*P - Q^2 + 1)/P^2*x^2 
       + (8*A*P - 12*Q)/P*x + 4/P;
  D := HyperellipticCurve(f);
  D := QuadraticTwist(D,delta);
  D := ReducedMinimalWeierstrassModel(D);
  D := SimplifiedModel(D);
  return D;
end function;

// Criterion for a genus 2 Jacobian to be absolutely simple,
// taken from page 158 of Cassels and Flynn.

function MyCount(g,h,p)
  P1 := PolynomialRing(GF(p));
  P2 := PolynomialRing(GF(p^2));
  Cp := HyperellipticCurve([P1!g,P1!h]);
  Cpp := HyperellipticCurve([P2!g,P2!h]);
  return #Cp,#Cpp;
end function;

function ProvedSimple(C,bd)
  P<x> := PolynomialRing(Rationals());
  PP<y> := PolynomialRing(P);
  CC := ReducedMinimalWeierstrassModel(C);
  eqn := Evaluate(Equation(CC),[x,y,1]);
  g := -Coefficient(eqn,0);
  h := Coefficient(eqn,1);
  assert CC eq HyperellipticCurve([g,h]);
  badp := BadPrimes(CC);
  pp := [p : p in [3..bd] | IsPrime(p) and p notin badp];
  ans := false;
  for p in pp do
    n1,n2 := MyCount(g,h,p);
    t := p + 1 - n1;
    n := (1/2)*(n1^2 + n2) - (p+1)*n1 - p;
    flag1 := IsSquare(t^2 - 4*n);
    poly := x^4 - (n/p)*x^3 + ((t^2 - 2*n - 2*p)/p)*x^2 - (n/p)*x + 1;
    flag2 := exists{k : k in [1..12] | GCD(x^k-1,poly) ne 1};
    if (not flag1) and (not flag2) then ans := true; end if;
  end for;
  return ans;
end function;

// Component group orders computed using Donnelly's programs 
// RegularModel and ComponentGroup

function ComponentGroupOrders(C)
  CC := ReducedMinimalWeierstrassModel(C);
  cgos := [];
  for p in BadPrimes(C) do
    try
      if p ne 2 then 
        M := RegularModel(C,p);  
      else
        M := RegularModel(CC,2);  
      end if;
      rr := #ComponentGroup(M);
    catch e;
      rr := 0;
    end try;
//    print p,rr;
    Append(~cgos,rr);
  end for;
  return cgos;
end function;

// Component group orders (and Tamagawa numbers) computed in 
// some simple cases where C has multiplicative reduction.

function Thickness(g,h,p)
  a0,a1,a2 := Explode([Coefficient(g,i): i in [0..2]]);
  b0,b1 := Explode([Coefficient(h,i): i in [0..1]]);
  assert forall{c : c in [a1,b0] | Valuation(c,p) gt 0};
  t := Valuation(a0,p);
  if t in {0,1} then return t,[0,0]; end if;
  delta := GF(p)!(b1^2 + 4*a2);
  assert delta ne 0;
  h := (1/p)*Evaluate(h,p*x);
  g := (1/p)^2*Evaluate(g,p*x);
  a0,a1,a2 := Explode([Coefficient(g,i): i in [0..2]]);
  b0,b1 := Explode([Coefficient(h,i): i in [0..1]]);
  x0 := Integers()!(-(1/delta)*(2*a1 + b0*b1));
  y0 := Integers()!((1/delta)*(a1*b1 - 2*a2*b0));
  g := g - h*y0 - y0^2;
  h := h + 2*y0;
  g := Evaluate(g,x + x0);
  h := Evaluate(h,x + x0);
  return Thickness(g,h,p) + 2;
end function; 

function MyComponentGroupOrder(C,p)
  C := ReducedMinimalWeierstrassModel(C);
  P<x> := PolynomialRing(Rationals());
  PP<y> := PolynomialRing(P);
  eqn := Evaluate(Equation(C),[x,y,1]);
  g := -Coefficient(eqn,0);
  h := Coefficient(eqn,1);
  thicknesses := [];
  R<X,Y> := PolynomialRing(GF(p),2);
  epslist := [];
  all_sing_pts_over_Fp := true;
  for ctr in [1..2] do
    gX := Evaluate(PZ!g,X);
    hX := Evaluate(PZ!h,X);
    Cp := Scheme(Spec(R),Y^2 + hX*Y - gX);
    assert Dimension(Cp) eq 1;
    Sp := SingularSubscheme(Cp);
    if Dimension(Sp) gt 0 then return 0,_; end if; // case not covered 
    singpts := [Eltseq(pt): pt in Points(Sp)];
    if Dimension(Sp) eq 0 and #singpts lt Degree(Sp) then
      all_sing_pts_over_Fp := false;
    end if;
    if ctr eq 2 then 
      singpts := [pt : pt in singpts | pt[1] eq 0]; 
    end if;
    for pt in singpts do
      aa,bb := Explode(ChangeUniverse(pt,Integers()));
      eqn1 := (y+bb)^2 + Evaluate(h,x+aa)*(y+bb) - Evaluate(g,x+aa);
      g1 := -Coefficient(eqn1,0);
      h1 := Coefficient(eqn1,1);
      a0,a1,a2 := Explode([Coefficient(g1,i): i in [0..2]]);
      b0,b1 := Explode([Coefficient(h1,i): i in [0..1]]);
      assert forall{c : c in [a0,a1,b0] | Valuation(c,p) gt 0};
      delta := GF(p)!(b1^2 + 4*a2);
      if delta ne 0 then 
	if p eq 2 then  
          eps := GF(2)!a2 eq 0 select 1 else -1;
        else
  	  eps := LegendreSymbol(Integers()!delta,p);
        end if;
        Append(~epslist,eps);
        tau := Thickness(g1,h1,p);
        Append(~thicknesses,tau);
      else
        return 0,_; // singular point is not a node
      end if;
    end for;
    g := P!(x^6*Evaluate(g,1/x));
    h := P!(x^3*Evaluate(h,1/x));
  end for;
  if #thicknesses eq 0 then
    if all_sing_pts_over_Fp then 
      return 1,1; // good reduction
    else
      return 0,_; // case not covered
    end if;
  end if;
  assert #thicknesses in [2,3];
  assert all_sing_pts_over_Fp;
  assert #Set(epslist) eq 1;
  if #thicknesses eq 2 then Append(~thicknesses,[0]); end if;
  t1,t2,t3 := Explode(thicknesses);
  GG<x1,x2,x3> := FreeAbelianGroup(3);
  RR := sub<GG|t1*x1 - t2*x2,t2*x2 - t3*x3,x1 + x2 + x3>;
  Phi := quo<GG|RR>;
  assert #Phi eq t1*t2 + t2*t3 + t3*t1;
  cp := epslist[1] eq 1 select #Phi else #{x : x in Phi | x eq -x};
  return #Phi,cp;
end function;

print "Some checking of Table 4";
//printf "(The messages \"WARNING: Extending the base field to ...\"";
//print " are generated by the call to \"RegularModel\")";
primelimit := 100; // Extending to 1000 takes a few minutes
pp := [p : p in [2..primelimit] | IsPrime(p)];
print "";

for ctr in [1..#table4] do

  gg := table4[ctr];
  crlab := gg[1];
  assert crlab in shalist;
  kappa := gg[2];
  A,P,Q,delta := Explode(gg[3]);
  ntors,rank := Explode(gg[4]);
  C := HyperellipticCurve(gg[7]);
  E := EllipticCurve(crlab);
  crlab1 := crlab;
  if exists{xx : xx in table1 | crlab eq xx[1] } then
    crlab1 := "(" cat crlab cat ")";
  end if;   
  printf "%o %o %o %o %o %o\n",crlab1,kappa,
    [A,P,Q,delta],[ntors,rank],gg[5],gg[6];

// Checking the ratio of levels
  N := Conductor(E);
  assert conductors[ctr] eq (kappa*N)^2;

// Checking the equation for C
  assert IsIsomorphic(C,Genus2Curve(A,P,Q,delta));

// Comparing the traces of Frobenius (up to primelimit)
  coeffs_f := [Coefficient(gg[7],i): i in [0..6]];
  aplist := TracesOfFrobenius(E,primelimit);
  badp := BadPrimes(E) cat BadPrimes(C);
  for i in [1..#pp] do
    p := pp[i];
    if p notin badp then
      Cp := HyperellipticCurve(PolynomialRing(GF(p))!coeffs_f);
      n1 := #Cp;
      n2 := #ChangeRing(Cp,GF(p^2));
      ap := aplist[i];
      trace := (p + 1) - n1;
      norm := Integers()!((1/2)*(n1^2 + n2)) - (p + 1)*n1 - p;
      assert (ap^2 - trace*ap + norm) mod 7 eq 0;
    end if;
  end for;

// Checking the Mordell-Weil group of F = Jac C.
  F := Jacobian(C);  
  pts := [elt<F|[xx[1],xx[2]],xx[3]>: xx in MWgens[ctr]];
  T,Tmap := TorsionSubgroup(F);
  assert ntors eq #T;
  assert rank eq #pts;
  assert rank in {0,4} 
     or (crlab in {"23025c","75712bz"} and rank eq 2)
     or (crlab in {"25296a"} and rank eq 6);

// Checking the Selmer rank (takes a few minutes in total)
  T := TorsionSubgroup(F);
  t := Valuation(#[x : x in T | 2*x eq T!0],2);
//  print "Computing 2-Selmer group";
//  S2,map := TwoSelmerGroup(F);
//  selrk := #Generators(S2) - t;
//  print "dim S^(2)(F/Q) =",#Generators(S2);
//  print "2-Selmer rank =",selrk;
//  assert selrk eq rank;

// Checking the generators are independent (takes less than a minute in total)
  if crlab notin ["14938n","20502ba","70950bc"] then 
    H := HeightPairingMatrix(pts);
    det := Determinant(H);
    assert det gt 10^(-5);  
  end if;

// Checking the component group orders (takes less than a minute in total)
//  cgos := ComponentGroupOrders(C);
//  assert cgos eq cgdata[ctr][4];

// Checking the primes with Tamagawa number divisible by 7
  assert forall{ c : c in TamagawaNumbers(E) | c mod 7 ne 0};
  badp := BadPrimes(C);
  assert cgdata[ctr][2] eq badp;
  assert forall{i : i in [2,3,4] | #cgdata[ctr][i] eq #badp};
  badcp := [];
  for i in [1..#badp] do
    p := badp[i];
    cgordersage := cgdata[ctr][3][i];
    cgordermagma := cgdata[ctr][4][i];
    cgordermult,cp := MyComponentGroupOrder(C,p);
    ccc := {cgordersage,cgordermagma,cgordermult};
    assert #ccc eq 1 or (#ccc eq 2 and 0 in ccc); // answers agree
    if p eq 2 and crlab in ["9510e","26166h","76608bl"] then
      continue; 
    // "9510e","26166h" : Genus 2 curve has "type C" multiplicative 
    // reduction. A separate calculation (not reproduced here) shows 
    // that the component group order is a power of 2.
    // "76608bl" : Genus 2 Jacobian has potential good reduction.
    end if;
    assert exists(cgorder){c : c in ccc | c ne 0};
    if cgorder mod 7 eq 0 then 
      assert cgordermult ne 0; 
      if cp mod 7 eq 0 then Append(~badcp,p); end if;
    end if;
  end for;
  assert badcp eq gg[5];

// Checking the reduction behaviour at p = 7
  redtype := gg[6];
  if redtype eq "g" then 
    assert ReductionType(E,7) eq "Good"; 
    assert 7 notin BadPrimes(C);
  elif redtype in ["m","m^*"] then 
    E1 := E;
    C1 := C;
    if redtype eq "m^*" then 
      E1 := MinimalModel(QuadraticTwist(E,7));
      C1 := ReducedMinimalWeierstrassModel(QuadraticTwist(C,7));
    end if;
    assert ReductionType(E1,7) in 
      {"Split multiplicative", "Nonsplit multiplicative"};
    assert MyComponentGroupOrder(C1,7) ne 0;
  elif redtype in ["g^{(3)}","g^{(4)}"] then 
    assert ReductionType(E,7) eq "Additive";
    r := redtype eq "g^{(3)}" select 3 else 4;
    PQ<X> := PolynomialRing(Rationals());
    K<u> := NumberField(X^r-7);
    assert Norm(Conductor(ChangeRing(E,K))) mod 7 ne 0;
    OK := RingOfIntegers(K);
    pr7 := Decomposition(OK,7)[1,1];
    _,resmap := ResidueClassField(pr7);
    PK<X> := PolynomialRing(K);
    P7<x> := PolynomialRing(GF(7));
    f := Evaluate(-Equation(C),[X,0,1]);
    assert P7!f ne 0;
    ff := Factorization(P7!f);
    a := Integers()!Coefficient(ff[1][1],0);
    if r eq 3 then 
      g := (1/7^2)*Evaluate(f,u*X-a);
    else
      g := (1/u^10)*Evaluate(f,u^2*X-a);
    end if; 
    assert forall{c : c in Eltseq(g) | c in OK};
    g7 := P7![resmap(c): c in Eltseq(g)];
    assert Discriminant(g7) ne 0;
  else
    assert redtype eq "";
    assert crlab in ["75712bz","84474w"];
    if crlab eq "75712bz" then 
      assert ReductionType(E,7) eq "Nonsplit multiplicative";
      assert 7 notin BadPrimes(C); // good reduction
    else
      assert ReductionType(E,7) eq "Good";
      assert 7 in BadPrimes(C); // potential good reduction
    end if;
  end if;

// Checking the Jacobian is absolutely simple
  assert ProvedSimple(C,50) or crlab eq "78400le";

end for;

// Remark 5.4

E := EllipticCurve("78400le");
K<u> := QuadraticField(-7);
a := -13462750;
b := -19012918750;
assert E eq EllipticCurve([a,b]);
R<x,y,z> := PolynomialRing(Rationals(),3);
FF := R!X7Twists(a,b);
M1 := Matrix(Rationals(),3,3,
	     [ -14700, 0, -17150, 31727500, -85750, 35500500, 7, 0, 8 ]);
F := (1/(2^3*5^9*7^9))*FF^M1;
assert F eq -14*x^3*y - 9*x^3*z + 21*x^2*y^2 - 6*x^2*y*z 
   + 18*x^2*z^2 - 14*x*y^3 + 36*x*y^2*z - 48*x*y*z^2 + 22*x*z^3 
   - 16*y^3*z + 18*y^2*z^2 + 8*y*z^3 + 4*z^4;

X := Curve(Proj(R),F);
pts := <[0,1,0],[1,0,0],[2,0,-3],[1+u,-1+u,0]>;
pts1 := [X(K)!P : P in pts];
pts := <[&+[M1[i,j]*P[j]: j in [1..3]]: i in [1..3]] : P in pts>;
assert forall{P : P in pts | Evaluate(FF,P) eq 0};

function Covariants(F)
  Deriv := Derivative;
  Hmat := Matrix(3,3,[Deriv(Deriv(F,i),j): i,j in [1..3]]);
  Psi_6 := (-1/54)*Determinant(Hmat);
  vv := Matrix(3,1,[Deriv(Psi_6,i): i in [1..3]]);
  HHmat := HorizontalJoin(Hmat,vv);
  HHmat := VerticalJoin(HHmat,Matrix(1,4,Eltseq(vv) cat [0]));
  Psi_14 := (1/9)*Determinant(HHmat);
  Jmat := Matrix(3,3,[Deriv(f,i): i in [1..3],f in [F,Psi_6,Psi_14]]);
  Psi_21 := (1/14)*Determinant(Jmat);
  return Psi_6,Psi_14,Psi_21;
end function;

cubicforms := [
      2*(a*x^2 + 3*b*x*z + 3*y^2 + 2*a*y*z)*z,
      2*(-b*x^3 - 2*a*x^2*y - 6*b*x*y*z - a*b*x*z^2 - 8*y^3 
         - 8*a*y^2*z - 2*a^2*y*z^2 + (4*a^3 + 28*b^2)*z^3),
      2*(-b*x^3 - a^2*x^2*z - 12*b*x*y*z - 8*a*b*x*z^2 - 2*y^3 
         + a*y^2*z + 4*a^2*y*z^2 + (4*a^3 + 16*b^2)*z^3),
      2*(3*a*b*x^3 - 2*a^2*x^2*y + (a^3 + 21*b^2)*x^2*z - 2*a*b*x*y*z 
	  - 2*a^2*b*x*z^2 - a^2*y^2*z + 6*b^2*y*z^2 + a*b^2*z^3) ];

E2 := EllipticCurve([ 0, -1, 0, -9581304833, 376825066377537 ]);
E3 := EllipticCurve([ 0, -1, 0, -30961949539208, -70450654083418969338 ]);
E4 := EllipticCurve([0,5*(7+u),0,70*(1+u),0]);
EE := <E,E2,E3,E4>;

for i in [1..4] do
  pt := pts[i];
  assert Evaluate(FF,pt) eq 0;
  Psi_6,Psi_14,Psi_21 := Covariants(FF);
  cc4 := Evaluate(Psi_14,pt);
  cc6 := Evaluate(Psi_21,pt);
  dd := [Evaluate(d,pt): d in cubicforms];
  assert exists(d){d : d in dd | d ne 0};
  assert forall{d1 : d1 in dd | d1 eq 0 or IsSquare(d*d1)};
  assert Evaluate(Psi_6,pt) ne 0;
  E1 := EllipticCurve([-27*cc4*d^2,54*cc6*d^3]);
  assert IsIsomorphic(E1,EE[i]);
end for;

print "---------------------------------------------";
print "Section 6. Checking the congruences";
print "";

for ctr in [1..35] do

  gg := table4[ctr];
  crlab := gg[1];
  kappa := gg[2];
  A,P,Q,delta := Explode(gg[3]);
  ntors,rank := Explode(gg[4]);
  f := gg[7];
  C := HyperellipticCurve(gg[7]);
  E := EllipticCurve(crlab);
  crlab1 := crlab;
  if exists{xx : xx in table1 | crlab eq xx[1] } then
    crlab1 := "(" cat crlab cat ")";
  end if;   
  printf "%o %o %o %o %o %o\n",crlab1,kappa,
    [A,P,Q,delta],[ntors,rank],gg[5],gg[6];

// Checking surjectivity of rhobar_{E,7}.  
  assert not HasComplexMultiplication(E);
  aaa := TracesOfFrobenius(E,300);
  ppp := [p : p in [2..300] | IsPrime(p)];
  ss := {};
  surj := false;
  NE := Conductor(E);
  for i in [1..#ppp] do
    p := ppp[i];
    ap := aaa[i];
    if (7*NE mod p) ne 0 and (ap mod 7) ne 0 then
      Include(~ss,LegendreSymbol(ap^2 - 4*p,7));
    end if;
  end for;
  assert {-1,1} subset ss;

  J := Jacobian(C);
  kk := kumsrf_pts[ctr];
  assert kk[1] eq CR(E);
  assert kk[2] eq f;
  assert #kk[4] eq 5;

// Our data file specifies an integral basis in a compact format, 
// which we now unpack.

  mat1 := Matrix(Rationals(),8,8,kk[3]);
  RR<aa,bb,cc> := PolynomialRing(Rationals(),3);
  LHS := [1,aa,bb,cc,aa^2,bb^2,cc^2,aa*bb,bb*cc,cc*aa,aa^3,bb^3];
  IdMat := IdentityMatrix(Rationals(),8);
  mat0 := Matrix(4,8,[IdMat[i]: i in [1..4]]);
  mat2 := VerticalJoin(mat0,mat1);
  mat3 := MyRowReduce(KernelMatrix(mat2));
  pols := [&+[mat3[i,j]*LHS[j]: j in [1..12]]: i in [1..4]];
  II := ideal<RR|pols>;
  X := Scheme(Spec(RR),II);
  assert Degree(X) eq 8;
  assert IsIrreducible(X);
  pol := Basis(EliminationIdeal(II,{aa}))[1];
  P<T> := PolynomialRing(Rationals());
  pol := Evaluate(pol,[T,0,0]);
  F<a> := NumberField(pol);
  aa,bb,cc := Explode(Eltseq(Points(X,F)[1]));
  LHS := [aa^2,bb^2,cc^2,aa*bb,bb*cc,cc*aa,aa^3,bb^3];
  matinv := mat1^(-1);
  basOF := [&+[matinv[i,j]*LHS[j]: j in [1..8]]: i in [1..8]];
  OF := LLL(MaximalOrder(Order(basOF)));
  M := Matrix(Integers(),8,8,[Eltseq(OF!x): x in basOF]);
  assert Determinant(M) in {-1,+1};
  assert Matrix(8,1,LHS) eq ChangeRing(mat1,F)*Matrix(8,1,basOF);

// Our point P on the Kummer surface is defined over L 
// where [L:K] = 3 and [K:Q] = 8. Our data file specifies each
// coordinate of P by giving its minimal polynomial over K.

// Looking up a point P on the Kummer surface
  PF<X> := PolynomialRing(F);
  cubicpolys := [&+[d[i][j]*basOF[j]*X^(4-i)
                       : i in [1..4],j in [1..8]]: d in kk[4]];
  K := KummerSurface(J);
  L<x0> := NumberField(cubicpolys[1]);
  KL := BaseChange(K,L);
  x1_list := [r[1]: r in Roots(cubicpolys[2],L)];
  x2_list := [r[1]: r in Roots(cubicpolys[3],L)];
  x3_list := [r[1]: r in Roots(cubicpolys[4],L)];
  xx_list := [r[1]: r in Roots(cubicpolys[5],L)];
  pts := [[x0,x1,x2,x3]: x1 in x1_list,x2 in x2_list,x3 in x3_list];
  assert #pts eq 27;
  ff := DefiningPolynomial(K);
  pts := [pt : pt in pts | Evaluate(ff,pt) eq 0];
  assert #pts eq 1;
  P := KL!pts[1];

// Checking 7*P = 0
  assert P ne KL![0,0,0,1];
  assert 7*P eq KL![0,0,0,1];
  assert AbsoluteDegree(L) eq 24;
  assert Degree(MinimalPolynomial(P[2]/P[1],Rationals())) eq 24;

// Checking P is defined over the field we want
  poly24 := DivisionPolynomial(E,7);
  assert IsIrreducible(poly24);
  assert Evaluate(poly24,xx_list[1]) eq 0;

end for;

print "---------------------------------------------";
print "Section 7. Proof of Theorem 1.2 ";
print "";

PQ<x> := PolynomialRing(Rationals());
A,P,Q := Explode([7,9,37]);
R := 4*P;
B := (Q*(P*A-Q)+4*P^2+1)/P^2;
C := 4*(P*A-Q)/P;
g := x^6 + A*x^5 + (A^2*P^3 - 3*A*P^2*Q + A*P*Q + 4*P^3 
       + 4*P^2 + 2*P*Q^2 - 2*P - Q^2 + 1)/P^2*x^4 
       + (A^2*P^2*Q + 4*A*P^2 - 3*A*P*Q^2 + A*P + 4*P^2*Q - 4*P*Q 
       + 2*Q^3 - 2*Q)/P^2*x^3 + (4*A^2*P^3 - 12*A*P^2*Q + A*P*Q 
       + 16*P^3 - 8*P^2 + 8*P*Q^2 + 4*P - Q^2 + 1)/P^2*x^2 
       + (8*A*P - 12*Q)/P*x + 4/P;

assert 9*(x+1)^6*Evaluate(g,(-x-2)/(x+1)) eq x*(x+4)*(x^4+2*x^3-x-3);

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

E := EllipticCurve("67080r1");

D := -191;
ED := MinimalModel(QuadraticTwist(E,D));
r,l := AnalyticRank(E:Precision:=40);
rD,lD := AnalyticRank(ED:Precision:=40);
assert [r,rD] eq [0,1];
w1,w2 := Explode(Periods(E));
HeegnerHeight := l*lD*Sqrt(Abs(D))/(2*w1*Im(w2));
print "h_K(y_K) =",HeegnerHeight;
// 217958.077779060986262802052278

R<x1,x2,x3,x4> := PolynomialRing(Rationals(),4);
Q1 := x1^2 + 6*x1*x2 + 7*x1*x3 + 23*x1*x4 - 3*x2^2 
           + 12*x2*x3 + 8*x2*x4 - 30*x3^2 + 54*x3*x4 - 60*x4^2;
Q2 := 6*x1^2 - x1*x2 + 89*x1*x3 - 37*x1*x4 + 37*x2^2 
           - 86*x2*x3 + 22*x2*x4 - 166*x3^2 + 114*x3*x4 + 199*x4^2;
QI := GenusOneModel([Q1,Q2]);
C4,E2,pi := nCovering(QI);
P := C4![ 847793227 , 2227947281 , 2665508726 , 1875455642 ];
_,psi := IsIsogenous(E2,ED);
P := psi(pi(P));

print "h_K(P) = 2*h(P) =",2*Height(P);
// 1112.03100907684176664694924631
print "h_K(y_K)/h_K(P) =",HeegnerHeight/(2*Height(P));
// 196.000000000000000000000000000

print "Checking P is a generator for ED(Q) mod torsion";

pts := Saturation([P],10^3);
pts := [P : P in pts | Order(P) eq 0];
assert pts eq [P];
bd := SiksekBound(ED:Torsion);
B := Ceiling(Exp(bd + 1));
pts1 := Points(ED: Bound := B);
assert pts1 eq {ED!0};

print "Equation for C (a twist of the Klein quartic)";

c4,c6 := Explode(cInvariants(E));
a := -27*c4;
b := -54*c6;
assert [a,b] eq [1013948336592,471906827379024768];
F,G := X7Twists(a,b);
R<x,y,z> := PolynomialRing(Rationals(),3);
M1 := Matrix(Rationals(),3,3,[102188568,-22341576,-103840044,
 73939908228192,94676367452256,135756936508464,-28,46,-151]);
F := (R!F)^M1;
F := 6*F/MonomialCoefficient(F,x^4);
assert F eq 6*x^4 + 25*x^3*y - 39*x^3*z - 3*x^2*y^2 
  - 15*x^2*y*z + 27*x^2*z^2 - 42*x*y^2*z - 141*x*y*z^2 
  - 87*x*z^3 + 8*y^4 - 6*y^3*z - 21*y^2*z^2 + 127*y*z^3 - 15*z^4;
C := Curve(Proj(R),F);
P1 := C![4,9,2];
P2 := C![5,-10,3];
print F;
M2 := Matrix(Rationals(),3,3,[ -2970312, -10081176, -3614004, 
 2, -4, 9, -5559004047744, 3336437174688, 10576760019552 ]);
G := (R!G)^M2;
G := 5*G/MonomialCoefficient(G,x^4);
assert G eq 5*x^4 - 26*x^3*y - 24*x^3*z + 36*x^2*y^2 - 54*x^2*y*z 
  + 36*x^2*z^2 - 78*x*y^3 - 12*x*y^2*z - 132*x*y*z^2 + 14*x*z^3 
  + 11*y^4 + 174*y^3*z + 120*y^2*z^2 - 10*y*z^3 - 90*z^4;
CC := Curve(Proj(R),G);
assert not IsLocallySoluble(CC,2);

jE := jInvariant(E);
R<X> := PolynomialRing(Rationals());
f := (X^2 + 5*X + 1)^3*(X^2 + 13*X + 49) - jE*X;
K1 := OptimizedRepresentation(NumberField(f));

mat1 := Matrix(Rationals(),8,8,[ 
  17, 9, 4, 2, -1, -2, -1, -3, 
  494, 14, -3, 18, -4, -2, -5, -3, 
  -67, -25, 51, -40, 17, -4, -10, 8, 
  69, 19, 1, 0, 1, -3, 5, -4, 
  541, -6, 5, -10, -8, 1, -4, -5, 
  50, 13, -1, -1, -1, -1, 3, -1, 
  1451, 211, 76, 32, -50, -69, 32, -46, 
  6960, 36, 759, -258, 83, -64, -60, -187 
]);
RR<aa,bb,cc> := PolynomialRing(Rationals(),3);
LHS := [1,aa,bb,cc,aa^2,bb^2,cc^2,aa*bb,bb*cc,cc*aa,aa^3,bb^3];
IdMat := IdentityMatrix(Rationals(),8);
mat0 := Matrix(4,8,[IdMat[i]: i in [1..4]]);
mat2 := VerticalJoin(mat0,mat1);
mat3 := MyRowReduce(KernelMatrix(mat2));
pols := [&+[mat3[i,j]*LHS[j]: j in [1..12]]: i in [1..4]];
II := ideal<RR|pols>;
X := Scheme(Spec(RR),II);
assert Degree(X) eq 8;
assert IsIrreducible(X);
pol := Basis(EliminationIdeal(II,{aa}))[1];
P<T> := PolynomialRing(Rationals());
pol := Evaluate(pol,[T,0,0]);
assert pol eq T^8 - 4*T^7 - 84*T^6 - 3920*T^5 - 74263*T^4
   - 408660*T^3 - 5409068*T^2 - 10750762*T - 12779595;

K := NumberField(pol);
aa,bb,cc := Explode(Eltseq(Points(X,K)[1]));
LHS := [aa^2,bb^2,cc^2,aa*bb,bb*cc,cc*aa,aa^3,bb^3];
matinv := mat1^(-1);
basOK := [&+[matinv[i,j]*LHS[j]: j in [1..8]]: i in [1..8]];
OK := LLL(MaximalOrder(Order(basOK)));
M := Matrix(Integers(),8,8,[Eltseq(OK!x): x in basOK]);
assert Determinant(M) in {-1,+1};
assert IsIsomorphic(K,K1);

print "Checking the class number is odd (under GRH)";

SetClassGroupBounds("GRH");
// SetVerbose("ClassGroup",4);
time AA,AAmap := ClassGroup(K);
assert #AA eq 1;

print "Checking some details of the 2-descent";

// A cubic form f with coefficients in K

fdata := [
    [ 102, -90, 48, -18, 66, -24, -24, 24 ],
    [ 235, -141, 176, -35, 191, -62, -34, -10 ],
    [ 216, -381, -159, 120, -192, -48, 60, 141 ],
    [ 83, -224, -5, 2, -34, -35, 71, -82 ],
    [ -56, 682, 162, 272, -34, 14, 86, -48 ],
    [ 211, -126, 279, 104, 63, -110, 264, -433 ],
    [ 45, -136, 15, -90, 83, -8, -58, 29 ],
    [ 404, -592, 54, 100, -59, -131, 181, -129 ],
    [ -313, -989, -524, -469, -207, 72, -72, -4 ],
    [ 13, -159, 54, -58, 33, -26, 72, -178 ]    
];

gdata := [
    [ -25, 21, 44, -22, -43, -20, -7, -15 ],
    [ 84, -117, 11, -58, 45, -44, -10, 2 ],
    [ 89, -162, 1, -4, -6, -26, -35, 131 ],
    [ 20, -8, -108, -14, -30, -24, 28, -24 ],
    [ -274, 522, -10, 24, -42, 48, -4, 72 ],
    [ 163, -201, -46, 16, 75, 116, -145, -89 ],
    [ -4, 38, -28, 10, 22, 4, -20, 4 ],
    [ -68, 88, 76, 84, 54, -92, -10, -58 ],
    [ 62, 15, 99, -218, 105, -64, -54, -86 ],
    [ 157, -216, 183, -302, 236, -10, 159, -31 ]
];

R3K := PolynomialRing(K,3);
mons := MD(R3K,3);
f := &+[fdata[i,j]*mons[i]*basOK[j]: i in [1..10],j in [1..8]];
g := &+[gdata[i,j]*mons[i]*basOK[j]: i in [1..10],j in [1..8]];
h := 12*x^2*z + 6*x*y^2 + 2*x*y*z + 15*x*z^2 
                 - 3*y^3 + 19*y^2*z - 4*y*z^2 + 3*z^3;
II := ideal<R3K|[Derivative(g,i): i in [1..3]]>;
YY := Scheme(Proj(R3K),II);
assert Dimension(YY) eq 0;
assert Degree(YY) eq 3;
Hmat := Matrix(3,3,[Deriv(Deriv(F,i),j): i,j in [1..3]]);
H := (-1/54)*Determinant(Hmat);
assert [R3K|F,H] subset II;
// The ternary cubic g defines the union of the three tangent
// lines at the points of a triangle (say T1).
// Therefore div(g) = 4 T1 - 3 H (where H = hyperplane section)

function MyDivide(D,n)
  SS,ee := Support(D);
  return &+[Integers()!(ee[i]/n)*SS[i]: i in [1..#SS]];
end function;  

function EffectiveDivisor(D)
  V,Vmap := RiemannRochSpace(D);
  D1 := Divisor(Vmap(V.1)) + D;
  SuppD1 := Support(D1); // computing this speeds things up later
  return D1; 
end function;

KC := FunctionField(C);
P := [Divisor(P1),Divisor(P2)];
H := Numerator(Divisor(KC!(x/y)));
assert Degree(H) eq 4;
D := MyDivide(Divisor(KC!(h/x^3))+3*H,2); 
assert Degree(D) eq 6;
assert IsEffective(D);
// Therefore div(h) = 2 D - 3 H

xx := [ 7798, 7367, 620, 7654, -6519, -602, 6523, -4513 ];
cc := (1/129)*&+[xx[i]*basOK[i]: i in [1..8]];
assert f^2 - cc*g*(R3K!h) in ideal<R3K|F>;
// Therefore div(f) = 2 T1 + D - 3 H

// Determining the set of bad primes S for the 2-descent.

Delta := Discriminant(E);
assert {p : p in BadPrimes(E) | Valuation(Delta,p) mod 7 eq 0} eq {5};
Rp := PolynomialRing(GF(5),3);
Cp := Curve(Proj(Rp),F);
assert not IsSingular(Cp); // C has good reduction at 5

// We take S = {p2} where p2 is the unique prime above 2.

p2 := Decomposition(OK,2)[1,1]; 
LTSM := LocalTwoSelmerMap(p2);
V := VectorSpace(GF(2),10);
LocalTest := func<x|V!Eltseq(x @ LTSM)>; 

// A basis for K(S,2).

xx := [
    [ -1, 0, 0, 0, 0, 0, 0, 0 ],
    [ 2, 0, 0, 0, 0, 0, 0, 0 ],
    [ 111, -2, 1, -9, -2, 0, 0, -1 ],
    [ 0, 8, -7, 8, -3, 1, -1, -4 ],
    [ 5, -8, 0, -2, 3, -1, 2, -2 ],
    [ -26, 12, -9, 9, -5, -1, 3, 1 ]
];

r1,r2 := Signature(K);
assert [r1,r2] eq [2,3];
KS2_basis := [&+[x[i]*basOK[i]: i in [1..8]]: x in xx];
assert [Valuation(x,p2): x in KS2_basis] eq [0,4,3,0,0,1];
assert forall{x : x in KS2_basis | IsSquare(x*OK) or IsSquare(x*p2)};
U := sub<V|[LocalTest(x): x in KS2_basis]>;
assert Dimension(U) eq 6;
// We have a 6-dimensional subspace of K(S,2), and assuming 
// the class number is odd, this is all of K(S,2).
assert [SquareFree(Integers()!Norm(x)): x in KS2_basis] eq [1,1,1,1,1,-1];
Uglob := sub<V|[LocalTest(KS2_basis[i]): i in [1..5]]>;
assert Dimension(Uglob) eq 5;

// The local 2-Selmer conditions at p = 2

p := 2;
Qp := pAdicField(p:Precision:=50);
PP<T> := PolynomialRing(Qp);
vv := [];
for pt in [[4,9,2],[5,-10,3],[1,0,2],[1,2,3]] do
  xx := [pt[1],pt[2],pt[3]+T];
  poly := Evaluate(F,xx);
  rts := [r[1] : r in Roots(poly) | Valuation(r[1]) gt 1];
  assert #rts eq 1;
  localpt := [pt[1],pt[2],pt[3]+rts[1]];
  localpt1 := C(Qp)!localpt;
  xx := Evaluate(f,ChangeUniverse(localpt,Integers()));
  Append(~vv,LocalTest(xx));
end for;
U0 := sub<V|[LocalTest(x): x in [2,-1,5]]>;
Uloc := U0 + sub<V|[vv[i]-vv[1]: i in [1..#vv]]>;
assert Dimension(Uloc) eq 5;
// This is the correct size for the local image (details omitted).

// Checking that the images of G1 = [Lambda - 2 P1] and 
// G2 = [P1 - P2] generate the 2-Selmer group

xx := Evaluate(f,[4,9,2])*Evaluate(f,[5,-10,3]);
assert xx ne 0;
assert IsSquare(86*xx*KS2_basis[3]);
Uknown := sub<V|[LocalTest(x): x in [-1,2,KS2_basis[3]]]>;
assert Dimension(Uknown) eq 3;
assert Uknown eq Uloc meet Uglob;

print "Chabauty's method with p = 7";

R := EffectiveDivisor(D - 3*P[1]);
G1 := (R + P[1] - H);
G2 := P[1] - P[2];
// To see that G1 is the same as above, note that the divisors
// 2 T1 + D - 3 H, H - 2 Lambda and 2 T1 - 3 \Lambda are principal.

function MyRowReducePolys(polys,d)  
  R := Universe(polys);
  mons := MonomialsOfDegree(R,d);
  mat := Matrix(#polys,#mons,[MC(q,mon):mon in mons,q in polys]);
  mat := MyRowReduce(mat);
  polys := [&+[mat[i,j]*mons[j]: j in [1..#mons]]:i in [1..Nrows(mat)]];
  return polys;
end function;

qq := MyRowReducePolys(Basis(Ideal(R)),2);

p := 7;

Rp<u,v,w> := PolynomialRing(GF(p),3);
PP<s> := PolynomialRing(GF(p));
Cp := Curve(Proj(Rp),Rp!F);
Sp := SingularSubscheme(Cp);
assert Dimension(Sp) eq 0;
assert Degree(ReducedSubscheme(Sp)) eq 1;
singpt := Points(Sp)[1];
assert singpt eq Cp![4,3,1];
assert Valuation(Evaluate(F,[4,3,1]),7) eq 1;
param := [
    4*s^4 + 2*s^2 + 2*s + 4,
    3*s^4 + 3*s^2 + 5*s + 5,
    s^4 + 4*s + 2
];
assert Evaluate(F,param) eq 0;

function GetParameter(pt)
  k := Universe(Eltseq(pt));
  kX<X> := PolynomialRing(k);
  ff := ChangeUniverse(param,kX);
  eqns := [ff[i]*pt[j] - ff[j]*pt[i]: i,j in [1..3] | i lt j];
  t := Roots(GCD(eqns))[1,1];
  assert Parent(pt)![Evaluate(f,t): f in ff] eq pt;
  return t;
end function;

function MyIntegral(f)
  P<X> := Parent(f);
  d := Degree(f);
  return &+[(1/i)*Coefficient(f,i-1)*X^i : i in [1..d+1]];
end function;
  
gg := [MyIntegral(f): f in param];
V := VectorSpace(GF(p),3);

function phi(pts)
  tt := [GetParameter(pt) : pt in pts];
  return V![&+[Evaluate(g,t): t in tt]: g in gg];
end function;

P1red := Cp![GF(p)!x : x in Eltseq(P1)];
P2red := Cp![GF(p)!x : x in Eltseq(P2)];

ptsH := Points(Scheme(Cp,Rp.3),GF(p^3));
assert #ptsH eq 4;
ptsR := Points(Scheme(Cp,[Rp!q: q in qq]),GF(p^3));
assert #ptsR eq 3;

phiG1 := phi(ptsR) +  phi([P1red]) - phi(ptsH);
phiG2 := phi([P1red]) - phi([P2red]);

assert phiG1 eq V![3,6,1];
assert phiG2 eq V![6,0,6];

W := sub<V|phiG1,phiG2>;
assert Dimension(W) eq 2;

pts := {P : P in Points(Cp) | not IsSingular(P)};
pts := {P : P in pts | phi([P]) - phi([P1red]) in W};
assert pts eq {P1red,P2red,Cp![1,-1,0]};
assert phi([Cp![1,-1,0]]) - phi([P1red]) eq phiG1 + 4*phiG2;

print "Power series expansion for points near P1 = (4 : 9 : 2)";

PR<aa> := PolynomialRing(Rationals());
RR<t> := PowerSeriesRing(PR);
zz := 0;
for r in [1..20] do
  gg := Evaluate(F,[2*(t + 2),9 + zz + aa*t^r,2]);
  cc := Roots(Coefficient(gg,r))[1][1];
  zz +:= cc*t^r;
end for;
RR<t> := PowerSeriesRing(Rationals());
zz := RR!zz + O(t^8);
pt := [t + 2,(9 + zz)/2,1];
assert Evaluate(F,pt) eq O(t^8);
print pt;

R2<u,v> := PolynomialRing(Rationals(),2);
g := Evaluate(F,[u,v,1]);
assert pt[3] eq 1;
gu := Derivative(g,u);
w3 := Derivative(pt[2])/Evaluate(gu,[pt[1],pt[2]]);
w1 := pt[1]*w3;
w2 := pt[2]*w3;
la := [Integral(w): w in [w1,w2,w3]];

print "Computing divisors D_i such that 7 G_i = [D_i - 3*P1]";

function MySupport(C,D)
  YY := Scheme(Ambient(C),Ideal(D));
  cc := PointsOverSplittingField(YY)[1];
  polys := [MinimalPolynomial(cc[i]/cc[3]): i in [1..3]];
  F := NumberField(polys[1]);
  assert Degree(F) eq 3;
  Q1 := C(F)![Roots(polys[i],F)[1][1]: i in [1..#polys]];
  return Eltseq(Q1);
end function;

RR<X> := PolynomialRing(GF(7));
data := [<G1,[6,3,4]>,<G2,[5,2,0]>,<G1+G2,[4,5,4]>];
for dd in data do
  G,ans := Explode(dd);
  D1 := EffectiveDivisor(G + 3*P[1]);      // G + 3*P1
  D2 := EffectiveDivisor(2*D1 - 3*P[1]);   // 2*G + 3*P1
  D3 := EffectiveDivisor(D1 + D2 - 3*P[1]);// 3*G + 3*P1
  D4 := EffectiveDivisor(2*D2 - 3*P[1]);   // 4*G + 3*P1
  D7 := EffectiveDivisor(D3 + D4 - 3*P[1]);// 7*G + 3*P1
  Q1 := MySupport(C,D7);
  t := (Q1[1]/Q1[3]) - 2;
  assert RR!MinimalPolynomial(t) eq X^3;
  for i in [1..3] do
    laD := Trace(Evaluate(la[i],t));
    assert Integers(7^2)!laD eq 7*ans[i];
  end for;
  assert (ans[1] + ans[2] + 3*ans[3]) mod 7 eq 0;
end for;

print "The Mordell-Weil sieve with p = 41";

p := 41;

Rp<u,v,w> := PolynomialRing(GF(p),3);
Cp := Curve(Proj(Rp),Rp!F);
assert not IsSingular(Cp);

J,_,Jmap := ClassGroup(Cp);
assert InvariantFactors(J) eq [2*3^2*7*547,0];

P1red := Cp![GF(p)!x : x in Eltseq(P1)];
P2red := Cp![GF(p)!x : x in Eltseq(P2)];
Kp := FunctionField(Cp);
H := Numerator(Divisor(Kp!(u/w)));
assert Degree(H) eq 4;  
mm := [Numerator(Divisor(Kp!((Rp!q)/Rp.1^2))): q in qq];
R := GCD(GCD(mm[1],mm[2]),mm[3]);
assert Degree(R) eq 3;

G1 := R + Divisor(P1red) - H;
G2 := Divisor(P1red) - Divisor(P2red);

r := GF(7)!Eltseq(Jmap(G1))[1];
psi := func<x|(GF(7)!Eltseq(Jmap(x))[1])/r>;
assert psi(G1) eq 1;
assert psi(G2) eq 1;

J1 := sub<J|Jmap(G1),Jmap(G2)>;
assert J1 eq sub<J|2*J.1>;

pts := Points(Cp);
assert #pts eq 42;
pts1 := { Q : Q in pts | Jmap(Divisor(Q)-Divisor(P1red)) in J1};
assert #pts1 eq 15;

ans := {psi(Divisor(Q)-Divisor(P1red)) : Q in pts1};
assert ans eq {0,1,2,6};
assert psi(G1 + 4*G2) notin ans;

