/*

VERSION: 18 February 2004

========================================================================
This Magma program checks whether or not there is a curve of genus three
over F_27 with no rational points.  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_27 having exactly 20 rational points.  Up to Galois
conjugacy, there are two E's to consider.  For each of these E's, we
must consider rational points Q that represent the classes of E(K)
modulo 2*E(K), where K = F_27.  We can choose representative Q's that
are not killed by 2.  Given E and Q, we list the degree-6 functions f 
on E that are regular outside infinity and that have a zero of degree
at least 2 at the point Q.  The covers we must consider are obtained 
by adjoining a root of z^2 = f to the function field of E.

The two E's we consider are
  y^2 = x^3 - x^2 + 1
  y^2 = x^3 - x^2 + a
where a^3 - a + 1 = 0.  Suppose we write Q = (x0,y0) with y0 nonzero,
so that the tangent line to E at Q is given by f0 = 0, where f0 is 
the linear polynomial 
  f0 = y0*(y - y0) + x0*(x - x0).
Then the degree-6 functions on E that we must consider are all of the
form
  f = (x - x0)^2 * (a degree-1 polynomial in x)
      + (x - x0) * (y - y0) * (constant)
      + f0 * (constant).
      
Let c be the (nonzero) coefficient of x in the degree-1 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^-6 + (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 -1.

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 y-resultant of f with the equation for E has exactly one
double root (corresponding to Q) and four simple roots.  If so, we
print out E, Q, and f, and indicate that this triple will require
further examination.  

As we will see, no triples get printed out.

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

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

Z<z>:=PolynomialRing(GF(3));
K<a>:=ext<GF(3) | z^3 - z + 1>;
R<x,y>:= PolynomialRing(K,2);

function firstcheck(f, E, Q)
// 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.
  for P in E do
    if not P in {Q,Identity(E)} then
      if IsSquare(Evaluate(f,[P[1],P[2]])) then
        return false;
      end if;
    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 6 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 6 intersection points, there are at
// most (6 choose 2) = 15 slopes to avoid.
//
// So our strategy will be to compute the projection of the 
// intersection divisor to the x-line via all 27 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 K 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 [4,1] and 0 eq Evaluate(f,[Q[1],Q[2]]);
end function;

function bothchecks(E,Qchoices)
// E is the elliptic curve.
// Qchoices is the list of Q's 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; 
  for Q in Qchoices do
    x0 := Q[1];
    y0 := Q[2];
    f0 := y0*(y-y0) + x0*(x-x0);
    // f0 is the tangent line to E at Q.  This depends on
    // E being of the form y^2 = x^3 - x^2 + c.
    
    // We will take f = (-1)*(f1 + c0*f0)
    // where f1 = (x - x0)^2 * (monic linear polynomial in x)
    //          + (x - x0) * (y - y0) * (constant).
    
    for c1, c2 in K do
      f1 := (x - x0)^2 * (x + c1) + (x - x0) * (y - y0) * c2;
      for c0 in K do
        f := -1*(f1 + c0*f0);
        if firstcheck(f,E,Q) then
          elist := Eltseq(E);
          eqn := y^2 + elist[1]*x*y + elist[3]*y 
                 - (x^3 + elist[2]*x^2 + elist[4]*x + elist[5]);
          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;
  end for;
  return bool;
end function;


// The two curves we must examine:

E1 := EllipticCurve([K | 0, -1, 0, 0, 1]);
E2 := EllipticCurve([K | 0, -1, 0, 0, a]);

twoE1 := {2*P : P in E1};
twoE2 := {2*P : P in E2};

Q1a := E1![1,1];
Q1b := E1![a^15,a];
Q1c := E1![a^19,a^16];
Q1d := E1![a^5,a^9];

// These four Q's represent the classes of E1(K) mod 2*E1(K):
assert 0 eq #{Q : Q in E1 | #({Q-Q1a,Q-Q1b,Q-Q1c,Q-Q1d} meet twoE1) eq 0};
assert #E1 eq 4*#twoE1;

Q2a := E2![a^5,a^14];
Q2b := E2![a^16,a^25];

// These two Q's repesent the classes of E2(K) mod 2*E2(K):
assert 0 eq #{Q : Q in E2 | #({Q-Q2a,Q-Q2b} meet twoE2) eq 0};
assert #E2 eq 2*#twoE2;

// Now do the checks:

check1 := bothchecks(E1,[Q1a,Q1b,Q1c,Q1d]);
check2 := bothchecks(E2,[Q2a,Q2b]);

if check1 and check2 then 
  print "There are no choices that require further analysis.";
else
  print "Unexpectedly, there were choices that require further analysis.";
end if;