SetColumns(0);

// the example from the paper
function elkies_quartic()
  P:=PolynomialRing(Rationals()); x:=P.1;
  return 19435071440*x^4-5351620404*x^3+130338882000*x^2-194951575764*x^1-357457601448;
end function;

function frob(u)
  D := LCM([Integers()!Denominator(uu) : uu in u]);
  u := [Integers()!(Abs(uu)*D) : uu in u];
  D := GCD(u);
  u := [uu div D: uu in u];
  return u;
end function;

// multiply by squares to kill denominator and really-obvious common factors
function quartic_simplify(C)
  dx := LCM([Integers()!Denominator(a) : a in Coefficients(C)]);
  dx1 := dx * Squarefree(dx); _,dd := Squarefree(dx1);
  C := C*dx1;
  dx := GCD([Integers()!t : t in Coefficients(C)]);
  a,b := Squarefree(dx);
  return C/b^2, dd/b;
end function;

function quartic_points(q)
  tbd := 10^2;
  while (2 ne 3) do
    print "before cputime call";
    k1 := Cputime();
    pts := Points(q : Bound:=tbd);
    k2 := Cputime();
    print "cputime was ",k2-k1;
    if (#pts gt 0 or (k2-k1) gt 10.0 or tbd gt 10^6) then
      return pts;
    end if;
    tbd := (3*tbd) div 2; // quadratic
    print "Increasing bound to ",tbd;
  end while;
end function;

P<x,y,z,t> := ProjectiveSpace(Rationals(),3);
cone1 := 17*x^2 + 7*y^2 - 26*y*z + 7*z^2;
cone2 := 13*y^2 - 7*y*z + 13*z^2 - 51*t^2;
PX := [6617,18454,61247,29797];
// the whole point of the Elkies construction is that we have
// a curve, defined as an intersection of two P3 quadrics,
// which happens to lie within the K3 surface x^4+y^4+z^4=1
// and that, being an intersection of two P3 quadrics,
// it's naturally elliptic

// so let's work more boringly

PS2<xx,yy,zz>:=ProjectivePlane(Rationals());
cone1 := 17*xx^2 + 7*yy^2 - 26*yy*zz + 7*zz^2;
cone2 := 13*xx^2 - 7*xx*yy + 13*yy^2 - 51*zz^2;

Evaluate(cone1,[PX[1],PX[2],PX[3]]);
Evaluate(cone2,[PX[2],PX[3],PX[4]]);

 CC := Conic(PS2,cone1);
 happy,pt := HasPoint(CC);
 if happy then
  print "Point on cone1 = ",pt;
  pm := Parametrization(CC,CC!pt);
  k1:=Inverse(pm)([PX[1],PX[2],PX[3]]);
  k1 := k1[1];
  print "Known point arises from k=",k1;
  P<k>:=PolynomialRing(Rationals());
  gp := pm([1,k]);
  
  xx := Evaluate(gp[1],1/k1);
  yy := Evaluate(gp[2],1/k1);
  zz := Evaluate(gp[3],1/k1);
  kkk := 13*yy^2 - 7*yy*zz + 13*zz^2;
  print IsSquare(kkk*51);
//  assert 5 eq 6;

  // y^2 * gp[3]^2 = c5*gp[1]^2 + c6*gp[1]*gp[3] + c7*gp[3]^2

  thing := 13*gp[2]^2 - 7*gp[2]*gp[3] + 13*gp[3]^2;
  crv := 51*thing;
  crv,den := quartic_simplify(crv);
  print "Hunting points on ",crv;
  // we now have to find a point on y^2 = -crv
  H := HyperellipticCurve(crv);
  time Ks := quartic_points(H);
  if (#Ks gt 0) then
   for kp in Ks do
    if (kp[3] ne 0) then
     kv := kp[1]/kp[3];
     x := Evaluate(gp[1],kv);
     y := Evaluate(gp[2],kv);
     z := Evaluate(gp[3],kv);
     t := Roots(k^4-((x^4+y^4+z^4)/18))[1][1];
     print frob([x,y,z,t]);
    end if;
   end for;
   // now use the EC to recover other roots
   E,m1 := EllipticCurve(H, Ks[1]);
   E2,m2 := MinimalModel(E);
   print "Elliptic curve underlying it all is ",E2," of rank ",Rank(E2);
   if (Rank(E2) ne 0) then
    ts,tsg := TorsionSubgroup(E2);
    // make torsion explicit
    ts := [tsg(r) : r in ts]; 
    gg := Generators(E2); 
    zpart := [r : r in gg | Order(r) eq 0];
    for t in [1..#ts] do
     for s in CartesianPower([-3..3],#zpart) do
      q := ts[t] + &+[s[i]*zpart[i] : i in [1..#s]];
      if (q ne E!0) then
       if (Inverse(m2)(q) ne E![0,0,1]) then
        kp := Inverse(m1)(Inverse(m2)(q));
        if (kp[3] ne 0) then
         kv := kp[1]/kp[3];
         xx := Evaluate(gp[1],kv);
         yy := Evaluate(gp[2],kv);
         zz := Evaluate(gp[3],kv);
         tt := Roots(k^4-((xx^4+yy^4+zz^4)/18))[1][1];
         J := frob([xx,yy,zz,tt]);
         print "Lift of ",<t,s>," is ",J;
        end if;
       end if;
      end if;
     end for;
    end for;
   end if;
  else
   print "Couldn't find points on ",H;
   print "Expected rank was ",Rank(AssociatedEC(-crv));
  end if;
 end if;