Ellpack

The file ellpack gives Maple code for computations with elliptic curves

  E: y^2 = f(x)

where

  f(x) = x^3 + ax^2 + bx + c.

The function Discr(f) computes the discriminant of f. When it is nonzero, E is nonsingular. The function ESum(f,P,Q), for P and Q distinct of E, computes the sum P+Q. The function EDup(f,P) computes 2P. For these computations we assume a, b, and c to be rational numbers, and we assume the same for the coordinates of P and Q.

For the above computations, we assume that the coefficients a, b, c are rational numbers, as are the coordinates of P and Q. It follows from the definition of the group law that R = P+Q has rational coordinates.

The function N(f,p) computes the number of points of E(Z/p), that is, of E considered over the field of integers modulo p. Thus we assume a, b, c to be integers.

Reference: Silverman and Tate, Rational Points on Elliptic Curves (Springer-Verlag 1992).

Notation: E(Z), is the set of points on E with integer coordinates. Likewise, E(Q), E(R), and E(C) are the sets of points on E with rational, real, and complex coordinates.

Examples and Discussion

1. The Discriminant

Let

  f := x -> x^3 - x;

Then Discr(f) = 4, and so E is nonsingular. To visualize the set of real points E(R), you can use pencil and paper or the following:

  with(plots):a := 2;
  b := 2;
  implicitplot( y^2 = f(x), x = -a..a, y = -b..b );

Now let

  f := x -> x^3 - x^2;

You will find that D = 0, so E is singular. In fact, a simple pencil and paper (or pure thought) analysis shows that E has a singular point at the origin. This point is a node: two branches of E(R) with distinct tangents pass through the origin.

2. The Group Law: Using ESum and EDup

The function ESum(f,P,Q) computes the sum R = P+Q for the elliptic curve E(Q). For example, if we set

  f := x -> x^3 + 17;
  P := [-1,4];
  Q := [2,5];

Then

  R := ESum(f,P,Q);

computes R = P + Q. In this case, R = (-8/9, -109/27).

As partial check on computations, set

  g := P -> P[2]^2 - f(P[1]);

and compute

  g(R);

The result should be zero.

Let us now compute S = 2P:

  S:= EDup(f,P);

It should be the case that g(S) = 0. In fact, S = (137/64, 2651/512).

3. Counting the number of points on E(Z/p)

To find the number of points on the elliptic curve

  E: y^2 = x^3 - x

modulo various primes, load the package and try this:

  f := x -> x^3 - x;

  N(f,5);
  N(f,7);
  N(f,11);
  N(f,13);

You should get 8, 8, 12, and 8 for N(f,p).

Discussion

The number of points of E(Z/p) is computed as follows. First, there is the point at infinity. Second, if f(x) = 0, then y = 0, so such an x contributes one point. Third, if f(x) is nonzero, then either it has two square roots mod p, or it has none. In the first case, f(x) is said to be a quadratic reside (it is the square of something). In the second case it is said not to be a quadratic residue. The Legendre symbol, computed using quadratic reciprocity, gives the necessary information:

  legendre(u,p) =  0 if u = 0 mod p
  legendre(u,p) = -1 if u is not a quadratic residue
  legendre(u,p) =  1 if u is a quadratic residue

Ellpack: the code

with(numtheory): # load number theory package

Discr := proc(f)
     local a,b,c;
     a := coeff(f(x),x,2);
     b := coeff(f(x),x,1);
     c := coeff(f(x),x,0);
     RETURN( -4*a^3*c + a^2*b^2 + 18*a*b*c - 4*b^3 - 27*c^2 );
 end;

ESum := proc(f,P,Q)
     local a,b,c,lambda,nu,xx,yy;
     a := coeff(f(x),x,2);
     b := coeff(f(x),x,1);
     c := coeff(f(x),x,0);
     lambda := (Q[2]-P[2])/(Q[1]-P[1]);
     nu := P[2] - lambda*P[1];
     xx := lambda^2 - a - P[1] - Q[1];
     yy := lambda*xx + nu;
     RETURN( [xx,-yy] );
 end:
 
 EDup := proc(f,P)
    local a,b,c,lambda,x,y,numx,denx,fprimex,xx,yy;
     a := coeff(f(x),x,2);
     b := coeff(f(x),x,1);
     c := coeff(f(x),x,0);
     x := P[1]; y := P[2];
     fprimex := subs(u=x,diff(f(u),u));
     lambda := fprimex/(2*y);
     
     numx := x^4 - 2*b*x^2 - 8*c*x + b^2 - 4*a*c;
     denx := 4*(x^3 + a*x^2 + b*x + c);
     xx := numx/denx;
     yy := y + lambda*(xx-x);
 
     RETURN( [xx,-yy] );
 end:

N := proc(f,p)
      local x, y, S;
      S := 0;
      for x from 0 to p-1 do
         y := f(x);
         S := S + 1 + legendre(y,p);
      od;
      RETURN (S + 1);
end: