
// A Magma file accompanying Section 6 of the paper:
//
// "EVERYWHERE LOCAL SOLUBILITY FOR HYPERSURFACES IN PRODUCTS
// OF PROJECTIVE SPACES"
// TOM FISHER, WEI HO, AND JENNIFER PARK
//
// Version : November 2019
//
// We report on some experiments to numerically investigate
// the proportion of (2,2)-forms over Q that are globally soluble.

// The following Magma code reads in the data for each
// experiment, performs some basic checking, and compiles
// the tables in Section 6.1 of the paper.

SetClassGroupBounds("GRH");
MC := MonomialCoefficient;
P1P1<X0,X1,Y0,Y1> := ProductProjectiveSpace(Rationals(),[1,1]);
monsx := [X0^2,X0*X1,X1^2];
monsy := [Y0^2,Y0*Y1,Y1^2];
mons := [mx*my : mx in monsx,my in monsy];
R<x1,x2> := PolynomialRing(Rationals(),2);
RR<y1,y2> := PolynomialRing(R,2);
discr := func<f|MC(f,y1*y2)^2 - 4*MC(f,y1^2)*MC(f,y2^2)>;

function IsRealSoluble(model:prec:=50)
// Testing binary quartics and quadric intersections
// for solubility over the reals
  n := Degree(model);
  Delta := Discriminant(model);
  error if Delta eq 0, "model is singular";
  if IsOdd(n) or Delta lt 0 then return true; end if;
  if n eq 2 then
    if #Eltseq(model) eq 8 then
      model := CompleteTheSquare(model);
    end if;
    F := Equation(model); 
    for xx in [[1,0],[0,1],[1,1],[1,-1]] do
      if Evaluate(F,xx) gt 0 then 
        return true;
      end if;
    end for;
    P := PolynomialRing(Integers()); X := P.1;
    f := Evaluate(F,[X,1]);
    R := RealField(prec);
    return HasRoot(f,R); 
  end if;
  assert Degree(model) eq 4;
  A,B := Explode(Matrices(model));
  for xx in [[1,0],[0,1],[1,1],[1,-1]] do
    if IsPositiveDefinite(xx[1]*A + xx[2]*B) then
      return false;
    end if;
  end for;
  P := PolynomialRing(Integers()); X := P.1;
  M := MatrixAlgebra(P,4);
  f := Determinant(X*M!A + M!B);
  R := RealField(prec);
  rts := Roots(f,R);
  if #rts eq 0 then return true; end if;    
  xx := [[(rts[i][1]+rts[i+1][1])/2,1] : i in [1..#rts-1]];
  for xx in [[1,0],[0,1],[1,1],[1,-1]] do
    if IsPositiveDefinite(xx[1]*A + xx[2]*B) then
      return false;
    end if;
  end for;
  return true;
end function;

function IsELS(model)
// Testing whether a genus one model over the rationals is
// Everywhere Locally Soluble.
  Delta := Integers()!Discriminant(model);
  pp := PrimeDivisors(Delta);
  badp := [ p : p in pp | not IsLocallySoluble(model,p:Transformation) ];
  if not IsRealSoluble(model) then badp := [-1] cat badp; end if;
  return #badp eq 0,badp;
end function;

function CTP(model1,model2)
// A function calling Donnelly's code for computing the
// Cassels-Tate Pairing
  dd := [Degree(model1),Degree(model2)];
  if dd eq [2,2] then
    C1 := HyperellipticCurve(model1);
    C2 := HyperellipticCurve(model2);
    print "Computing 2 x 2 Cassels-Tate pairing";
    return CasselsTatePairing(C1,C2);
  elif dd eq [2,4] then
    C1 := HyperellipticCurve(model1);
    C2 := Curve(model2);
    print "Computing 2 x 4 Cassels-Tate pairing";
    return CasselsTatePairing(C2,C1);
  end if;   
  error "CTP not implemented";
end function;

function IsIndependent(E,models)
// Checking that a list of binary quartics are independent as
// elements of H^1(Q,E[2]). (Not an efficient implementation,
// but sufficient for present purposes.)
  assert forall{ m : m  in models | Degree(m) eq 2};
  BQlist := [GenusOneModel(2,E)];
  for BQ in models do
    if exists{BQ0 : BQ0 in BQlist | IsEquivalent(BQ,BQ0:E := E)} then
      return false;
    end if;
    BQlist cat:= [BQ*BQ0: BQ0 in BQlist];
  end for;
  return true;
end function;

// Loading in the data for the experiments

data := [];
load "experiment-10.m";
data cat:= [<forms,solution_data,rank_data,qilist,ctp24data>];
load "experiment-30.m";
data cat:= [<forms,solution_data,rank_data,qilist,ctp24data>];
load "experiment-100.m";
data cat:= [<forms,solution_data,rank_data,qilist,ctp24data>];
load "experiment-300.m";
data cat:= [<forms,solution_data,rank_data,qilist,ctp24data>];
load "experiment-1000.m";
data cat:= [<forms,solution_data,rank_data,qilist,ctp24data>];

// The elliptic curves for which we ran EightDescent:

ell8desc := [
    [ 1, 1, 0, 89334, 75786117 ],
    [ 1, 1, 0, -965051, -306769654 ]
];

// The elliptic curves for which we ran analytic_rank_upper_bound:

ellrank1or3 := [
    [ 0, 0, 0, -385216, -118546643 ],
    [ 1, -1, 0, -21940631, -10062163381 ],
    [ 1, -1, 0, -130106786, -418444299752 ],
    [ 0, 1, 0, -674939767, 9768411280745 ],
    [ 1, 0, 1, 1365438724, 1088450102306 ],
    [ 0, -1, 0, -2310473362, 35959254765224 ], // (2,2)-form not ELS
    [ 1, 1, 0, 13646956, 36868880351052 ],
    [ 0, 0, 1, 463718380, -1653282652263 ],
    [ 1, 1, 0, -6811523942, 180704627470189 ]
];

sagetable := [
<2.00000000000000,  [1, 3, 3, 5, 3, 3, 3, 5, 5]>,
<2.10000000000000,  [1, 3, 3, 5, 3, 3, 3, 3, 5]>,
<2.20000000000000,  [1, 3, 3, 3, 3, 3, 3, 3, 5]>,
<2.30000000000000,  [1, 3, 3, 3, 3, 3, 3, 3, 5]>,
<2.40000000000000,  [1, 3, 3, 3, 3, 3, 3, 3, 5]>,
<2.50000000000000,  [1, 3, 1, 3, 3, 3, 3, 3, 5]>,
<2.60000000000000,  [1, 1, 1, 3, 3, 3, 3, 3, 3]>,
<2.70000000000000,  [1, 1, 1, 3, 3, 3, 3, 3, 3]>,
<2.80000000000000,  [1, 1, 1, 3, 3, 3, 3, 3, 3]>,
<2.90000000000000,  [1, 1, 1, 3, 3, 1, 3, 3, 3]>,
<3.00000000000000,  [1, 1, 1, 3, 1, 1, 3, 3, 3]>,
<3.10000000000000,  [1, 1, 1, 3, 1, 1, 3, 3, 3]>,
<3.20000000000000,  [1, 1, 1, 3, 1, 1, 3, 3, 3]>,
<3.30000000000000,  [1, 1, 1, 3, 1, 1, 3, 3, 3]>,
<3.40000000000000,  [1, 1, 1, 3, 1, 1, 3, 3, 3]>,
<3.50000000000000,  [1, 1, 1, 1, 1, 1, 3, 3, 3]>,
<3.60000000000000,  [1, 1, 1, 1, 1, 1, 1, 3, 3]>,
<3.70000000000000,  [1, 1, 1, 1, 1, 1, 1, 3, 3]>,
<3.80000000000000,  [1, 1, 1, 1, 1, 1, 1, 3, 3]>,
<3.90000000000000,  [1, 1, 1, 1, 1, 1, 1, 1, 3]>
];

// sage: E1 = EllipticCurve([ 0, 0, 0, -385216, -118546643 ])
// ....: E2 = EllipticCurve([ 1, -1, 0, -21940631, -10062163381 ])
// ....: E3 = EllipticCurve([ 1, -1, 0, -130106786, -418444299752 ])

// sage: E1.analytic_rank_upper_bound(max_Delta=2.0,adaptive=True)
// 1
// sage: E2.analytic_rank_upper_bound(max_Delta=2.5,adaptive=True)
// 3
// sage: E2.analytic_rank_upper_bound(max_Delta=2.6,adaptive=True)
// 1
// sage: E3.analytic_rank_upper_bound(max_Delta=2.5,adaptive=True)
// 1

for i in [1..#ellrank1or3] do
  EE := EllipticCurve(ellrank1or3[i]);
  ss := Sprint(EE);
  ss := &cat[ss[i] : i in [26..#ss-19]];  
  printf "%o & ",ss;
  dd := [d[1] : d in sagetable | d[2,i] eq 1];
  if #dd ne 0 then 
    printf "%o ",RealField(2)!Minimum(dd);
  end if;
  print "\\\\";
end for;

/*
 y^2 = x^3 - 385216*x - 118546643  & 2.0 \\
 y^2 + x*y = x^3 - x^2 - 21940631*x - 10062163381  & 2.6 \\
 y^2 + x*y = x^3 - x^2 - 130106786*x - 418444299752  & 2.5 \\
 y^2 = x^3 + x^2 - 674939767*x + 9768411280745  & 3.5 \\
 y^2 + x*y + y = x^3 + 1365438724*x + 1088450102306  & 3.0 \\
 y^2 = x^3 - x^2 - 2310473362*x + 35959254765224  & 2.9 \\
 y^2 + x*y = x^3 + x^2 + 13646956*x + 36868880351052  & 3.6 \\
 y^2 + y = x^3 + 463718380*x - 1653282652263  & 3.9 \\
 y^2 + x*y = x^3 + x^2 - 6811523942*x + 180704627470189  & \\
*/

// The tables in Section 6.1 :

answers := [
<10,
[ 753, 755, 753, 753, 885 ],
[ 5, 0, 86, 344, 313, 116, 21, 0, 885 ],
[ 5, 0, 16, 285, 310, 116, 21, 0, 753 ]>,
<30,
[ 640, 652, 642, 642, 885 ],
[ 0, 0, 122, 310, 291, 129, 29, 4, 885 ],
[ 0, 0, 1, 208, 274, 126, 29, 4, 642 ]>,
<100,
[ 536, 582, 549, 549, 875 ],
[ 0, 0, 171, 321, 257, 96, 25, 5, 875 ],
[ 0, 0, 0, 205, 221, 93, 25, 5, 549 ]>,
<300,
[ 378, 502, 432, 433, 867 ],
[ 0, 0, 230, 373, 187, 58, 19, 0, 867 ],
[ 0, 0, 0, 210, 151, 52, 19, 0, 432 ]>,
<1000,
[ 275, 464, 357, 464, 879 ],
[ 0, 0, 292, 366, 169, 44, 8, 0, 879 ],
[ 0, 0, 0, 175, 132, 43, 7, 0, 357 ]>
];

for exno in [1..5] do

  H := [10,30,100,300,1000][exno];
  forms,solution_data,rank_data,qilist,ctp24data := Explode(data[exno]);
  assert #forms eq 1000;
  assert #solution_data eq 1000;
  assert #rank_data eq 1000;

  print "Checking QIlist";
  QIlist := [GenusOneModel(Rationals(),4,coeffs): coeffs in qilist];
  time assert forall{QI : QI in QIlist | IsELS(QI)};
  BQfrom4Sel := [DoubleGenusOneModel(QI): QI in QIlist];
  function IsInImageOf4Selmer(BQ)
    if exists{ BQ0 : BQ0 in BQfrom4Sel | IsEquivalent(BQ,BQ0) } then
      return true;
    end if;
    print "Running 4-descent";
    time fd := FourDescent(BQ);
    return #fd gt 0;
  end function;

  bounds1 := [0,1000];                 // initial estimate
  bounds2 := [0,1000];                 // improved estimate
  nels := 0;                           // number which are ELS
  table := ZeroMatrix(Integers(),2,9); // breakdown by rank

  for ctr in [1..1000] do
    F := &+[forms[ctr,i]*mons[i] : i in [1..9]];
    sd := solution_data[ctr];
    print ctr,F;
    F1 := discr(Evaluate(F,[x1,x2,y1,y2]));
    BQ := GenusOneModel(F1); 
    if Discriminant(BQ) eq 0 then
      assert Evaluate(F,sd) eq 0;
      print "Solution =",P1P1!sd;
      soln := "found";
      els := true;
      rank := -1; // to indicate curve is singular
    else
      els,badp := IsELS(BQ);

// Information about the Mordell-Weil group E(Q)

      E := MinimalModel(Jacobian(BQ));
      r2,r4,r8,xx := Explode(rank_data[ctr]); 
      gens := [Points(E,x)[1] : x in xx];
      assert Regulator(gens) gt 0.1; // (Reg <> 0) => indep
      assert RootNumber(E) eq (-1)^r2;
      assert r2 - r4 in [0,2,4];
      assert r4 - #gens in [0,1,2,3];
      if H lt 1000 then
        if r4 - #gens in [0,1] then
	  assert r4 - r8 eq 0;
	else  
          assert exists(d){d : d in ctp24data | d[1] eq ctr};
          assert #TorsionSubgroup(E) eq 1;
          m2 := GenusOneModel(Rationals(),2,d[2]);
          m4 := GenusOneModel(Rationals(),4,d[3]);
          assert IsELS(m2) and IsELS(m4);
          BQlist := [GenusOneModel(2,P): P in gens];
          BQlist cat:= [m2,DoubleGenusOneModel(m4)];
          assert IsIndependent(E,BQlist);
          assert r4 - r8 eq (d[4] eq 0 select 0 else 2);  
	end if;
        if Eltseq(E) in ellrank1or3 then
  	  assert [r2,r4,r8,#gens] eq [3,3,3,1];
        else
	  assert r8 - #gens in [0,1];
	end if;
	rankbounds := [#gens,r8]; // lower and upper bounds
      else // experiment with H = 1000
        assert r8 eq -1; // not computed in this case
	rankbounds := [#gens,r4]; // lower and upper bounds
      end if;
      // using parity to improve the lower bound
      if IsOdd(r2 - rankbounds[1]) then rankbounds[1] +:= 1; end if;
      rank := rankbounds[1]; // our tables use the lower bound

// Information about whether C = { F = 0 } < P^1 x P^1 is soluble
    
      if #sd eq 4 then
        assert Evaluate(F,sd) eq 0;
        print "Solution =",P1P1!sd;
        soln := "found";
	assert els;
        if Eltseq(E) in ell8desc then soln := "eightdescent"; end if;
      elif not els then
        assert sd eq badp;
	print "Not locally soluble at",badp;
        soln := "notELS";
	assert Eltseq(E) notin ellrank1or3
	    or Eltseq(E) eq [0, -1, 0, -2310473362, 35959254765224];
      elif #sd eq 5 then
        BQ1 := GenusOneModel(Rationals(),2,sd);
        assert IsELS(BQ1);
        assert CTP(BQ,BQ1) eq 1;
	print "No solution by 2 x 2 Cassels-Tate pairing";
        soln := "CTP22";
      else
        assert IsInImageOf4Selmer(BQ);
        if #sd eq 20 then
          QI1 := GenusOneModel(Rationals(),4,sd);
          assert IsELS(QI1);
          // assert CTP(BQ,QI1) eq 1;  // can take a few minutes
   	  print "No solution by 2 x 4 Cassels-Tate pairing";
          soln := "CTP24";
        elif Position(ellrank1or3,Eltseq(E)) in [1..8] then
          print "No solution via analytic rank bound";
	  rankbounds[2] := 1;
          soln := "analytic";
        elif r4 - #gens eq 1 then
          print "Soluble by parity conjecture";
          soln := "parity";
          assert r2 - #gens eq 1 
              or (H eq 1000 and [r2,r4,#gens] eq [4,2,1]);
	else
          soln := "unresolved";
        end if; 
      end if;

      if #TorsionSubgroup(E) ne 1 then
        assert soln in {"found","notELS"};
      elif els then 
        V := VectorSpace(GF(2),#gens);
        vecs := [Eltseq(v): v in V | v ne V!0];
        pts := [&+[Integers()!(v[i])*gens[i] : i in [1..#gens]]: v in vecs];
        BQsol := [GenusOneModel(2,E)] cat [GenusOneModel(2,P): P in pts];
        flag := exists{BQ0 : BQ0 in BQsol | IsEquivalent(BQ,BQ0)};
	assert flag eq (soln in {"found","eightdescent"});
        if H eq 300 and r8 eq #gens and not flag and soln eq "unresolved" then
	  soln := "CTP24";
        end if;
      end if;
      
    end if; // Discriminant = 0

// Bookkeeping to produce the tables in the paper

    assert soln in {"found","CTP22","CTP24","notELS","analytic",
                    "eightdescent","parity","unresolved"};
    if els then nels +:= 1; table[1,rank + 2] +:= 1; end if;
    if soln in {"found"} then bounds1[1] +:= 1; end if;
    if soln in {"CTP22","notELS"} then bounds1[2] -:= 1; end if;
    if soln in {"found","eightdescent","parity"} then
       bounds2[1] +:= 1; table[2,rank + 2] +:= 1; end if;
    if soln in {"CTP22","CTP24","analytic","notELS"} then
       bounds2[2] -:= 1;
    end if;

  end for;

// Putting totals in the right hand column

  table[1,9] := &+[table[1,i] : i in [1..8]];
  table[2,9] := &+[table[2,i] : i in [1..8]];

  print "";
  print "H =",H;
  printf "Initial estimate %o - %o    Improved estimate %o - %o    #ELS  %o",
       bounds1[1],bounds1[2],bounds2[1],bounds2[2],nels;
  print "";
  print table;
  print "";

  assert answers[exno] eq <H,bounds1 cat bounds2 cat [nels],
       Eltseq(table[1]),Eltseq(table[2])>;

  print "press [RETURN] to continue";
  read xxx;

end for;


/*

// The unresolved cases described in the last paragraph of the paper.

{* [x[1],x[2],#x[4]] where x is rank_data[i] : i in [1..1000] |
             solution_data[i] eq [] *};
	     
{*
    [ 2, 2, 1 ]^^71,
    [ 3, 3, 1 ]^^90,
    [ 3, 3, 2 ]^^8,
    [ 4, 2, 1 ]^^3,
    [ 4, 4, 1 ],
    [ 4, 4, 2 ]^^14,
    [ 5, 3, 1 ],
    [ 5, 5, 3 ]
*}

*/