/*

VERSION: 18 February 2004

========================================================================
This Magma program finds all pointless genus-4 curves over F_53.
In the paper

  Everett W. Howe, Kristen E. Lauter, and Jaap Top:
  Pointless curves of genus three and four,
  preprint, 2004
  
it is shown that such a curve must be a double cover of an elliptic 
curve E over F_53 having 42 rational points.  There are exactly four E's
with exactly 42 points. For each such E, we must consider rational 
points Q whose images in E(K) modulo 3*E(K) meet all of the orbits
of this group under the action of Aut E.  This means that for each E
we will consider 2 Q's, one lying in 3*E(K) and one not.  We will 
choose our Q's to not be 2-torsion elements.

Suppose that we are given an E and a Q.  Then we list the degree-8
functions f that are regular outside the identity and that have a
double zero at Q.  The covers we must consider are obtained by 
adjoining a root of z^2 = f to the function field of E.

The four E's we consider can all be written in the form
  y^2 = x^3 + a*x + b.
Suppose we write Q = (x0,y0).  Then the tangent line to E at Q is given
by f0 = 0, where f0 is the linear polynomial 
  f0 = 2*y0*(y - y0) - (3*x0^2 + a)*(x - x0).
The degree-8 functions on E that we must consider are all of the form
  f = (x - x0)^2 * (a degree-2 polynomial in x)
      + (x - x0) * (y - y0) * (a polynomial in x of degree < 2)
      + f0 * (constant).
      
Let c be the (nonzero) coefficient of x^2 in the degree-2 polynomial
of x appearing above, and let t be a uniformizing parameter of E at
infinity (so we may take t = y/x, for example).  The Laurent expansion
of f at infinity is c*t^-8 + (higher-order terms), so the function f
is a square locally at infinity if and only if c is a square.  Since
we are searching for double covers of E with no rational points, we
may assume that c is not a square, and by multiplying f by a constant
square (which doesn't change the cover), we may assume that c is 2.

For each (E, Q, f) we then check to see whether there is a rational 
point P on E, different from Q and infinity, for which f(P) is a 
square.  Such a point would give a rational point on the double cover
associated to f.  If there is no such point P, we then check to see
whether the intersection divisor of f with the equation for E has 
exactly one double point (namely Q) and six simple points.  If so,
we print out E, Q, and f, and indicate that this triple will require
further examination. 

We will represent our f's as vectors in K^6.  A vector 
[c1,c2,c3,c4,c5,c6] will correspond to the function
c1 * x^2 * (x-x0)^2 
  + c2 * x * (x-x0)^2
  + c3 * (x-x0)^2
  + c4 * x * (x-x0) * (y-y0)
  + c5 * (x-x0) * (y-y0)
  + c6 * f0.
We will normalize our functions so that the value at the first
point on E is a fixed nonsquare.  To list the functions, we will 
choose 5 nonsquares in K, assign those values to the next 5 affine
points on E, solve for the vector [c1,...,c5,c6], and then check to
see whether the function evaluates to nonsquares on the other points
of E (including the infinite point).

========================================================================
*/

/*
VERSION: 18 February 2004.
First public version.
*/

K := GF(53);
nonsquare := K!2;
assert not IsSquare(nonsquare);

nonsquares := [x : x in K | not IsSquare(x)];

R<x,y>:=PolynomialRing(K,2);

function firstcheck(fvector, pointvalues)
// Here we check whether f(P) is a nonzero square for any point
// P on E (other than Q or infinity).
// Return false if we find such a point, true if not.
// pointvalues is a list of vectors, one for each affine point 
// on E other than Q plus one for the infinite point.  For the
// affine points, the vector is the vector of values of the
// basis functions at P.  For the infinite point, the vector
// is [1,0,0,0,0,0], so that we check to see whether the leading
// coefficient of the funcion is a nonsquare.
  for P in pointvalues do
    if IsSquare((fvector,P)) then
      return false;
    end if;
  end for;
  return true;
end function;

function secondcheck(f, eqn, Q)
// Given a polynomial f in K[x,y] and an equation for E in K[x,y], 
// we check to see that f, viewed as a function on E, has a double
// root at Q and four simple roots elsewhere.  Return false if we
// know that our f does not have this shape.  Return true if we have
// to look more closely at f to decide.
//
// We assume that the degree of f is equal to 8 and that it is
// regular outside infinity.
//
// Instead of using Magma's curve facilities, we will consider
// various projections of the intersection divisor (in P^2) of
// f and the equation for E onto the x-line.  Generically, the
// projection will have the same shape as the original intersection
// divisor --- the only way it can change is if we project at
// a slope that happens to be the slope of the line between two of
// the intersection points.
//
// Since there are at most 7 intersection points, there are at
// most (7 choose 2) = 21 slopes to avoid.
//
// So our strategy will be to compute the projection of the 
// intersection divisor to the x-line via 22 different slopes.  
// The form of the generic intersection divisor is then easy to spot:
//  If we list the divisors as 
// [{# of simple intersections}, {# of double intersections}, etc.]
// then the generic form will be the lexicographically greatest.

  formlist := [];
  for slope in [0..21] do
    newf := Evaluate(f,[x+slope*y,y]);
    neweqn := Evaluate(eqn,[x+slope*y,y]);
    pol := UnivariatePolynomial(Resultant(newf, neweqn, y));
    polfac := Factorization(pol);
    n1 := 0;
    n2 := 0;
    for fac in polfac do
      case fac[2]:
         when 1: n1 +:= Degree(fac[1]);
         when 2: n2 +:= Degree(fac[1]);
      end case;
    end for;
    formlist cat:= [ [n1,n2] ];
  end for;
  Sort(~formlist); // default order is lexicographic.
  realform := formlist[#formlist];
  return realform eq [6,1] and 0 eq Evaluate(f,[Q[1],Q[2]]);
end function;


function bothchecks(E,Q)
// E is the elliptic curve.
// Q is the Q to try.
// List the appropriate f's, and apply tests to (E,Q,f) triples.
// Print the triples that pass both tests.
// Return true if no triples pass, false if not.

  bool := true; 
  elist := Eltseq(E);
  eqn := y^2 + elist[1]*x*y + elist[3]*y 
             - (x^3 + elist[2]*x^2 + elist[4]*x + elist[5]);
  x0 := Q[1];
  y0 := Q[2];

  for m in K do
    b := y0 - m*x0;
    pol := UnivariatePolynomial(Evaluate(eqn,[x,m*x+b]));
    if pol mod UnivariatePolynomial(x-x0)^2 eq 0 then
      m0 := m;
      b0 := b;
      break;
    end if;
  end for;
  f0 := y - m0*x - b0;
  // f0 defines the tangent line to E at Q.
  
  // Make the list of values of the basis functions.
  // First the special vector for the infinite point.
  notQ := [ [K|1,0,0,0,0,0] ];
  // Then the vectors for the affine points other than Q.
  notQ cat:= [ [Evaluate(g,[P[1],P[2]]) :
                 g in [x^2 * (x-x0)^2 ,
                       x * (x-x0)^2,
                       (x-x0)^2,
                       x * (x-x0) * (y-y0),
                       (x-x0) * (y-y0),
                       f0]
                ] : P in E | not P in {Identity(E),Q}];
                
  matrix := Transpose(Matrix(K,6, &cat[notQ[i] : i in [1..6]]));
  testpoints := [Vector(K,6,notQ[i]) : i in [7..#notQ]];
 
  for v1 in nonsquares do
  time 
  for v2, v3, v4, v5 in nonsquares do
    valuevector := Vector(K,6,[nonsquare,v1,v2,v3,v4,v5]);
    fvector := Solution(matrix,valuevector);
     
    if firstcheck(fvector,testpoints) then
      // Convert the vector representation of f into a polynomial.
      f :=  fvector[1] * x^2 * (x-x0)^2 
          + fvector[2] * x * (x-x0)^2
          + fvector[3] * (x-x0)^2
          + fvector[4] * x * (x-x0) * (y-y0)
          + fvector[5] * (x-x0) * (y-y0)
          + fvector[6] * f0;
      if secondcheck(f,eqn,Q) then
        bool := false;
        print "This one needs further analysis: ";
        print E,Q,f;
        print "";
      end if;
    end if;
  end for; end for;
  return bool;
end function;


E1 := EllipticCurve([K|0,0,0,1,20]);
Q1a := E1![5,16];
Q1b := E1![8,13];

E2 := EllipticCurve([K|0,0,0,1,24]);
Q2a := E2![12,11];
Q2b := E2![ 0,17];

E3 := EllipticCurve([K|0,0,0,2,5]);
Q3a := E3![3,12];
Q3b := E3![2,21];

E4 := EllipticCurve([K|0,0,0,8,43]);
Q4a := E4![1,23];
Q4b := E4![0,19];

// Verify that these E's are in the right isogeny class.
assert (#E1 eq 42) and (#E2 eq 42) and (#E3 eq 42) and (#E4 eq 42);

// Verify that we are getting all four elements of the isogeny class.
assert #{jInvariant(E) : E in [E1,E2,E3,E4]} eq 4;

// Verify that the classes of E(K)/3E(K) are represented
// by the Q's and their negatives.
assert #&join[ {P+Q1a,P+Q1b,P-Q1b} : P in E1] eq 42;
assert #&join[ {P+Q2a,P+Q2b,P-Q2b} : P in E2] eq 42;
assert #&join[ {P+Q3a,P+Q3b,P-Q3b} : P in E3] eq 42;
assert #&join[ {P+Q4a,P+Q4b,P-Q4b} : P in E4] eq 42;




print "Checking first curve...";
ck1 := bothchecks(E1,Q1a) and bothchecks(E1,Q1b);
print "Checking second curve...";
ck2 := bothchecks(E2,Q2a) and bothchecks(E2,Q2b);
print "Checking third curve...";
ck3 := bothchecks(E3,Q3a) and bothchecks(E3,Q3b);
print "Checking fourth curve...";
ck4 := bothchecks(E4,Q4a) and bothchecks(E4,Q4b);

print"";

if ck1 and ck2 and ck3 and ck4 then 
  print "There are no pointless genus-4 curves over F_53.";
else
  print "I unexpectedly found a function that needs further study.";
end if;