AN EXPOSITION OF THE SEA ALGORITHM

JÁNOS A. CSIRIK

Abstract. The Schoof–Elkies–Atkin algorithm is an efficient way to count the number of points on an elliptic curve defined over a large prime field. This expository paper describes the algorithm in sufficient detail to allow a reader not familiar with arithmetic geometry to implement the algorithm. The mathematical background for the technique is then given.

1. Introduction

Let p be a large (odd) prime and let

E : y^{2} = x^{3} + a_{4}x + a_{6}

be an elliptic curve, where a_{4} and a_{6} are given fixed integers. In the case where p does not divide 4a_{4}^{3} +27a_{6}^{2}, E can be reduced to an elliptic curve over 𝔽_{p}. The number of points of E over 𝔽_{p}, denoted by #E(𝔽_{p}), is of cryptographic interest, since the properties of this number determine the security of elliptic curve cryptosystems based on E against various known attacks.

The first polynomial time algorithm for determining the number of rational points on an elliptic curve defined over a finite field is due to Schoof. He used calculations with torsion points on the curve to arrive at the number of points. At first Schoof’s algorithm was considered impractical, but Elkies suggested the use of “good” primes (now known as Elkies primes), where isogenies and modular curves can be involved to speed up the calculation. Atkin also made a number of important contributions to the algorithm, which then became known as the Schoof–Elkies–Atkin (SEA) algorithm. Further improvements were later proposed by Dewaghe and Couveignes–Dewaghe–Morain. The SEA algorithm was implemented by Morain, Müller, and Izu et al.

Schoof’s seminal paper [18] describes the original algorithm. He later also published a paper [19] that is a lovely overview of the developments in the subject up to 1995. Elkies’ paper [9] describes the ideas of his original manuscript [8] and contains many other theoretical insights and illuminating examples. The implementations of Morain and Müller are described in [15] and [16]. The implementation of Izu, Kogure, Noro and Yokoyama, which focuses on speeding up the algorithm as much as possible, is described in [13]. Dewaghe’s improvement is published in [7], see also Section 3.13. The improvement by Couveignes–Dewaghe–Morain is published in [5], see also Section 4.3. Atkin never formally published his contributions described in [1], but they are discussed extensively in [9, 19].

This paper, which is not aimed at the experts in the area, describes in detail a reasonably fast implementation of the SEA algorithm that is closely modeled upon Morain’s. The algorithm considered below is probabilistic and, for a 200-bit prime p, succeeds with a probability of about 3/4 (which can be brought arbitrarily close to 1 by enlarging the set A of auxiliary primes below). The algorithm implemented on a typical personal computer takes several minutes to find the number of points on a typical curve over 𝔽_{p}, where p has 200 bits.

It is known that

#E(𝔽_{p}) = p + 1 − t,

where t is an integer which satisfies the Hasse bound

−2\sqrt{p} ≤ t ≤ 2\sqrt{p}.

The algorithm works by calculating t modulo several small auxiliary primes ℓ. When the product of the auxiliary primes exceeds 4\sqrt{p}, the Chinese Remainder Theorem is used to recover the exact value of t, and hence that of #E(𝔽_{p}).

The algorithm works its way though a fixed list of 40 candidates for auxiliary primes given below. For each candidate ℓ, a calculation has to be carried out to generate a certain polynomial Ψ_{ℓ} that is necessary for further calculations with this ℓ. These polynomials Ψ_{ℓ} do not depend on the curve E under consideration and hence might be precomputed and stored if memory allows. Then for any elliptic curve E we can quickly decide if our algorithm applies (the probability that the algorithm applies for a specific E and ℓ is 1/2). For those curves where the algorithm applies, we can determine t modulo ℓ.

When we finished with all our candidates for the auxiliary primes, we can look at the elliptic curve and check whether the product of auxiliary primes that worked exceeds 4\sqrt{p} or not. In the former case, we succeeded in determining t.

A typical application for this point counting would be to take a random prime p and a random elliptic curve E over 𝔽_{p}, with the intention of finding an E with #E(𝔽_{p}) = xr, where r is a prime and x is small. Given such a curve, a point P of order r can be located easily and the pair (E, P) could be used for a number of cryptographic algorithms, such as Diffie-Hellman key exchange, El Gamal encryption, etc. If we use 200-bit primes for p and require x ≤ 32, then the probability that #E = xr is about 2.5%, so we expect to have to run our algorithm on about 55 curves.^{1}

Section 2 describes the algorithm in detail. Section 3 presents the mathematical background of the algorithm. Section 4 presents ideas by which the algorithm could be improved. Section 5 contains certain tables of data that need to be hardwired into a program implementing this algorithm.

2. The Algorithm

2.1. Overview. The set A of potential auxiliary primes is the union of the set A_{s} of small primes and the set A_{l} of larger primes. For each ℓ ∈ A, we need to determine a polynomial Ψ_{ℓ}(F, J) ∈ ℤ[F, J]. For ℓ ∈ A_{s}, this is stored in the program. For ℓ ∈ A_{l}, Ψ_{ℓ} must be calculated by determining a number of coefficients of a certain q-series f(q) ∈ ℤ[[q]] and carrying out certain algebraic operations on it. The polynomials Ψ_{ℓ} do not depend on the elliptic curve under consideration and therefore may be pre-calculated and stored if there is enough space for them (they require just under a half megabyte to store).

We start out with a given prime p and an elliptic curve

E : y^{2} = x^{3} + a_{4}x + a_{6}.

1The chance that a randomly selected integer of similar size is of the form xr as above is about 3%. For more details, see [12, 10]

We need to check a few simple things and calculate the j-invariant j = j(E) = 6912a_{4}^{3}/(4a_{4}^{3} + 27a_{6}^{2}) ∈ 𝔽_{p}.

Working over 𝔽_{p}, we plug in J = j ∈ 𝔽_{p} into Ψ_{ℓ}(F, J) and find all roots f ∈ 𝔽_{p} that satisfy

Ψ_{ℓ}(f, j) = 0.

For each of these f, we need to find all j̃ such that

Ψ_{ℓ}(f, j̃) = 0.

For all possible pairs (f, j̃), we do various operations involving Ψ_{ℓ} and its partial derivatives to obtain the quantities ã_{4}, ã_{6}, p_{1}. For any quintuplet (a_{4}, a_{6}, ã_{4}, ã_{6}, p_{1}), we can determine whether it is “valid”, and if so, we can generate a “kernel polynomial” h(X) ∈ 𝔽_{p}[X] of degree d = (ℓ − 1)/2. Given a kernel polynomial h(X), we can determine the “eigenvalue” e ∈ 𝔽_{p} and conclude that t ≡ e + p/e (mod ℓ).

An application of the Chinese Remainder Theorem completes the calculation.

The remaining subsections of Section 2 explain the details of the steps described above.

2.2. Determination of Ψ_{ℓ} for ℓ ∈ A_{s}. Set A_{s} = {3, 5, 7, 11, 13, 17, 19, 23}. For ℓ ∈ A_{s}, the polynomial Ψ_{ℓ} is to be found in the table of Section 5.3.

2.3. Determination of f for ℓ ∈ A_{l}. Let A_{l} be the set of prime numbers from 29 through 197, except for 37, 43, 67 and 163. We need to calculate using q-series (power series in the variable q, possibly with a finite number of negative powers of q). Let σ_{k}(n) denote the sum of kth powers of the positive divisors of n. The following q-series will be needed.

E_{4}(q) = 1 + 240 ∑_{n=1}^{∞} σ_{3}(n)q^{n} = 1 + 240q + 2160q^{2} + 6720q^{3} + ⋯

E_{6}(q) = 1 − 504 ∑_{n=1}^{∞} σ_{5}(n)q^{n} = 1 − 504q − 16632q^{2} − 122976q^{3} − ⋯
j(q) = 1728E_{4}(q)^{3}/(E_{4}(q)^{3} − E_{6}(q)^{2}) = q^{−1} + 744 + 196884q + 21493760q^{2} + ⋯
η̄(q) = ∏_{n=1}^{∞}(1 − q^{n}) = ∑_{k=−∞}^{∞}(−1)^{k}q^{(3k^{2}+k)/2} = 1 − q − q^{2} + q^{5} + ⋯

For each ℓ, we need to look up the suitable polynomial P_{ℓ} from the table in Section 5.1. Let v denote the degree of P_{ℓ}. Calculate the coefficients a_{k} ∈ 𝔽_{ℓ} such that

P_{ℓ}(j(q))η̄(q)η̄(q^{ℓ}) ≡ ∑_{k=−v}^{2ℓv−v} a_{k}q^{k} + O(q^{2ℓv−v+1}) (mod ℓ).

(This can be done by multiplying out the product on the left hand side modulo ℓ and reading off the coefficients.) For each a_{k}, let b_{k} be the least absolute remainder of a_{k} mod ℓ (the integer with the smallest possible absolute value that reduces to a_{k} modulo ℓ). Now the f we are searching for is the q-series

( ∑_{k=−v}^{2ℓv−v} b_{k}q^{k} + O(q^{2ℓv−v+1})) /(η̄(q)η̄(q^{ℓ})).

This will enable us to determine the coefficients of q^{k} in f(q) for all k ≤ 2vℓ − v.

2.4. Determination of Ψ_{ℓ} for ℓ ∈ A_{l}. This calculation involves manipulation of q-series. To simplify notation, we write all q-series as though all the coefficients were known, but in fact we need to calculate only as many coefficients as our knowledge of f from the previous step will allow. The previous step calculated just enough coefficients in f for this step to work correctly (see Section 3.14 for details).

For the calculation that follows, we shall need the q-expansions of f, f^{2}, …, f^{ℓ}. Let

(f(q))^{k} = ∑_{n} a(n, k)q^{n}.

For each 1 ≤ k ≤ ℓ, let

s_{k}(q) = ∑_{n} ℓa(ℓn, k)q^{n}.

For each 1 ≤ k ≤ ℓ, let

c_{k}(q) = − (s_{k}(q) + ∑_{r=1}^{k−1} c_{k−r}(q)s_{r}(q)) /k.

Finally, let C_{1}(q) = −f + c_{1}(q), C_{ℓ+1}(q) = −fc_{ℓ}(q) and for each 2 ≤ k ≤ ℓ,

C_{k}(q) = −fc_{k−1}(q) + c_{k}(q).

It turns out that for each 1 ≤ k ≤ ℓ + 1, there is a polynomial G_{k} such that

G_{k}(j(q)) ≡ C_{k}(q) (mod p).

Since the q-expansion of j starts with q^{−1} + ⋯, the coefficients of G_{k} are easy to determine. Indeed, let d be the order of the pole of C_{k} at q = 0, so that C_{k}(q) = aq^{−d} +⋯, where a is some non-zero constant. Then C_{k}(q) = aj(q)^{d} +C′_{k}(q), where C′_{k}(q) = C_{k}(q) − aj(q)^{d} is a polynomial in j(q), with a pole of order strictly less than d. Iterating this procedure, we obtain all the coefficients of G_{k}.

Once we determined all the G_{k}, we can write down Ψ_{ℓ} as

Ψ_{ℓ}(F, J) = F^{ℓ+1} + ∑_{i=1}^{ℓ} G_{i}(J)F^{ℓ+1−i}.

2.5. Initialization. The prime p is intended to be large: indeed, we need to check that ℓ < p for all ℓ ∈ A for the algorithm to work at all.

We also need to confirm that p does not divide 4a_{4}^{3} + 27a_{6}^{2}, otherwise E is not an elliptic curve when reduced modulo p.

Finally, calculate and store the j-invariant j = j(E) = 6912a_{4}^{3}/(4a_{4}^{3}+27a_{6}^{2}) ∈ 𝔽_{p}. We discard E if its j-invariant equals 0 or 1728. This condition would only be

satisfied for a very small number of curves. Also later in the calculation, we shall abandon the branch of the calculation where an elliptic curve with j-invariant 0 or 1728 arises.

We also discard E if it is supersingular. Supersingular curves are the ones for which #E(𝔽_{p}) = p+1. We can check E for supersingularity by taking a few random points on it and checking if p+1 annihilates them. If not, then the curve is definitely not supersingular and we can go on. If yes, then E is probably supersingular so we throw it away. Once again, we only throw away a few E this way. Furthermore, this check needs to be carried out only for our original E and not for any other curves that arise during the calculations later on.

2.6. Determination of the Possible Pairs (f, j̃). Since we are working with a single ℓ here, let us denote Ψ_{ℓ} by Ψ. Let Ψ_{1} denote the first partial derivative of Ψ, let Ψ_{12} denote the second partial derivative of Ψ_{1}, and so on.

If we plug in J = j into Ψ(F, J), we obtain a univariate polynomial Ψ(F, j) in the variable F. The roots f ∈ 𝔽_{p} of this polynomial are our candiates for f. If Ψ_{1} or Ψ_{2} vanishes at (f, j), then that f must be discarded from our list. This will happen very rarely but must be caught to avoid a division by zero later on.^{2}

For each one of our remaining candidates f ∈ 𝔽_{p}, we should find all the roots j̃ ∈ 𝔽_{p} of the polynomial Ψ(f, J). These are our candidates for j̃ with the given f. If Ψ_{1} or Ψ_{2} vanishes at (f, j̃), then this pair (f, j̃) must be dropped from our list. This will happen very rarely but must be caught to avoid a division by zero later on.

Also discard all the pairs (f, j̃) where j̃ = 0 or j̃ = 1728.

2.7. Determination of ã_{4}, ã_{6}, p_{1}, Given (f, j̃). We need to calculate the following quantities, which are elements of 𝔽_{p}. Let E_{4} = −48a_{4}, E_{6} = 864a_{6} and

f′ = \frac{E_{6}}{E_{4}} j \frac{Ψ_{2}(f, j)}{Ψ_{1}(f, j)},  Q = \frac{f′}{j̃} \frac{1}{ℓ} \frac{Ψ_{1}(f, j̃)}{Ψ_{2}(f, j̃)},  Ẽ_{4} = \frac{j̃}{j̃− 1728}Q^{2}.

Let Ẽ_{6} = Ẽ_{4}Q and

t_{1} = \frac{1}{Ψ_{1}(f, j)} (−f′Ψ_{11}(f, j) + 2jΨ_{12}(f, j)\frac{E_{6}}{E_{4}} − \frac{E_{6}^{2}}{f′E_{4}^{2}} (jΨ_{2}(f, j) + j^{2}Ψ_{22}(f, j))) ,

t_{2} = \frac{1}{Ψ_{1}(f, j̃)} (−f′Ψ_{11}(f, j̃) + 2ℓj̃Ψ_{12}(f, j̃)\frac{Ẽ_{6}}{Ẽ_{4}} − ℓ^{2} \frac{Ẽ_{6}^{2}}{f′Ẽ_{4}^{2}} (j̃Ψ_{2}(f, j̃) + j̃^{2}Ψ_{22}(f, j̃))) .

Let

t_{3} = \frac{E_{6}}{3E_{4}} − \frac{E_{4}^{2}}{2E_{6}} , t_{4} = ℓ(\frac{Ẽ_{6}}{3Ẽ_{4}} − \frac{Ẽ_{4}^{2}}{2Ẽ_{6}} ).

Finally, let p_{1} = ℓ(t_{2} + t_{4} − t_{1} − t_{3})/4, ã_{4} = −ℓ^{4}Ẽ_{4}/48 and ã_{6} = ℓ^{6}Ẽ_{6}/864.

All of the intermediate results, except for ã_{4}, ã_{6} and p_{1}, can be discarded before moving on.

2.8. Determination of the Kernel Polynomial h(X) given ã_{4}, ã_{6} and p_{1}. Recall our notation d = (ℓ − 1)/2. Here we are aiming to produce a polynomial h(X) of degree d. We need to pick an integer S, which is a small positive integer determining the number of “extra” terms we carry for a certain calculation. These extra terms will have to vanish at the end to assure us that the triple (ã_{4}, ã_{6}, p_{1}) we started with is “valid”. The choice S = 3 works fine.

First we set

p_{0} = d,
p_{2} = ((1 − 10d)a_{4} − ã_{4})/30,
p_{3} = ((1 − 28d)a_{6} − 42p_{1}a_{4} − ã_{6})/70.

2If one, but not both, of Ψ_{1} and Ψ_{2} vanishes at (f, j) then our formulae can be rewritten to make sense and avoid division by zero. However, even writing this footnote seems like excessive effort in the face of an eventuality this unlikely.

Recall that p_{1} was already given as part of our data. Then set

c_{1} = 6p_{2} + 2a_{4}d,
c_{2} = 10p_{3} + 6a_{4}p_{1} + 4a_{6}d,

and for each 2 ≤ r ≤ d − 1 + S, let

c_{r+1} = \frac{3 ∑_{n=1}^{r−1} c_{n}c_{r−n} − (2r − 1)(r − 1)a_{4}c_{r−1} − (2r − 2)(r − 2)a_{6}c_{r−2}}{(r − 1)(2r + 5)} .

Next, for each 3 ≤ n ≤ d − 1 + S, let

p_{n+1} = \frac{1}{4n + 2} (c_{n} − (4n − 2)a_{4}p_{n−1} − (4n − 4)a_{6}p_{n−2}) .

The quantity p_{i} is the sum of the ith powers of the roots of h(X). From here, we can obtain the coefficients in the usual manner (see [4, Prop. 4.3.3.]). Specifically, set s_{0} = 1, and for all 1 ≤ i ≤ d + S, let

s_{i} = \frac{−1}{i} ∑_{k=1}^{i}(−1)^{k}p_{k}s_{i−k}.

We need to check if s_{d+1} = s_{d+2} = ⋯ = s_{d+S} = 0. If not, we conclude that our triplet (ã_{4}, ã_{6}, p_{1}) was not valid, and we go on to the next (f, j̃) pair.

If our triplet was deemed to be valid, then

h(X) = ∑_{i=0}^{d}(−1)^{i}s_{i}X^{d−i}.

If the triplet (ã_{4}, ã_{6}, p_{1}) is not valid, then we can go on to the next pair (f, j̃) to get our next triplet (ã_{4}, ã_{6}, p_{1}), and so on. If no valid triplet (ã_{4}, ã_{6}, p_{1}) can be found, then the algorithm does not apply to the elliptic curve under consideration and this ℓ, so we find out nothing about t modulo ℓ.

2.9. Determination of the Eigenvalue e Given the Kernel Polynomial h(X). First of all, factor h(X) and replace it by any one of its (non-trivial) factors. From now on, h(X) will denote the factor that we chose. Note that in this step, the algorithm can be speeded up appreciably by rearranging the order in which some steps are carried out. For details, see Section 4.1.

Use the table in Section 5.2 to pick an integer s based on the value of ℓ. For any ℓ, we get an s ≤ 11. Given this s, we need four polynomials, a_{s}(X), b_{s}(X), c_{s}(X) and d_{s}(X). The way to obtain these is explained in Section 2.10.

The procedure we need to follow here is slightly different for the cases where the degree of h is even or odd. Let us assume for now that deg(h) is even.

Start by calculating and storing Q_{1}(X) = X^{p} mod h(X) and Q_{2}(X) = (X^{3} + a_{4}X + a_{6})^{(p−1)/2} mod h(X).

Throughout this calculation, we need to maintain a pair of polynomials (P_{1}(X), P_{2}(X)) and an integer e. Initially, (P_{1}(X), P_{2}(X)) = (X mod h(X), 1) and e = 1.

If at some stage (P_{1}(X), P_{2}(X)) = (Q_{1}(X), ±Q_{2}(X)) then we may stop. If P_{2}(X) = Q_{2}(X) then the current value of e is the eigenvalue we are aiming to calculate. If P_{2}(X) = −Q_{2}(X) then the eigenvalue is −e.

Otherwise, do the following replacements simultaneously (in other words, the values of P_{1} and P_{2} to be used on the right hand side are the old ones):

P_{1}(X) ← \frac{a_{s}(P_{1}(X))}{b_{s}(P_{1}(X))} mod h(X),
P_{2}(X) ← P_{2}(X)\frac{c_{s}(P_{1}(X))}{d_{s}(P_{1}(X))} mod h(X),
e ← es mod ℓ

Iterate this until the condition (P_{1}(X), P_{2}(X)) = (Q_{1}(X), ±Q_{2}(X)) is satisfied. This should occur after no more than (ℓ − 1)/2 steps.

In the case where deg(h) is odd, the calculation is the same except that we do not have to calculate Q_{2} or P_{2}. Whenever P_{1}(X) = Q_{1}(X) is achieved, the correct value of the eigenvalue is

e^{s(e)} (\frac{r}{ℓ} ) e,

where r is the resultant of h(X) and w(X) = X^{3} + a_{4}X + a_{6}, and s(x) is the semi-order of x modulo ℓ, i.e., the smallest positive n such that x^{n} ≡ ±1 (mod ℓ).

2.10. Determination of the Polynomials a_{s}(X), b_{s}(X), c_{s}(X) and d_{s}(X). The polynomials we will obtain depend on E, but not on ℓ. Since they are of relatively low degree, it is probably a good idea to calculate and store them once for each E and just look up the right ones for each ℓ.

Choose an integer R such that we need to calculate the polynomials a_{s}(X), b_{s}(X), c_{s}(X) and d_{s}(X) for s ≤ R. For example, R = 11 will do with the set of ℓ that are used in this algorithm. (But R = 5 will suffice with the modification described in Section 4.2.)

We need to calculate a number of polynomials with coefficients in 𝔽_{p}.

w(X) = X^{3} + a_{4}X + a_{6},
f_{1}(X) = 1,
f_{2}(X) = 2,
f_{3}(X) = 3X^{4} + 6a_{4}X^{2} + 12a_{6}X − a_{4}^{2},
f_{4}(X) = 4X^{6} + 20a_{4}X^{4} + 80a_{6}X^{3} − 20a_{4}^{2}X^{2} − 16a_{4}a_{6}X − 4a_{4}^{3} − 32a_{6}^{2}.

We can now determine the desired polynomials for s = 2 as

a_{2}(X) = 4Xw(X) − f_{3}(X),
b_{2}(X) = 4w(X),
c_{2}(X) = f_{4}(X)/4,
d_{2}(X) = 8w(X)^{2}.

For each 5 ≤ n ≤ R + 2, we need to calculate the polynomial f_{n}(X). If n = 2m, let

f_{n} = f_{m}(f_{m+2}f_{m−1}^{2} − f_{m−2}f_{m+1}^{2})/2.

If n = 2m + 1 with m even, let

f_{n} = w^{2}f_{m+2}f_{m}^{3} − f_{m−1}f_{m+1}^{3}.

If n = 2m + 1 with m odd, let

f_{n} = f_{m+2}f_{m}^{3} − w^{2}f_{m−1}f_{m+1}^{3}.

Now for any odd^{3} s in the range 2 < s ≤ R, we can write down

a_{s}(X) = Xf_{s}(X)^{2} − w(X)f_{s−1}(X)f_{s+1}(X),
b_{s}(X) = f_{s}(X)^{2},
c_{s}(X) = f_{s+2}(X)f_{s−1}(X)^{2} − f_{s−2}(X)f_{s+1}(X)^{2},
d_{s}(X) = 4f_{s}(X)^{3}.

2.11. The Last Step. For each elliptic curve E we end up knowing the important integer t modulo various small primes ℓ. Since the absolute value of t is bounded by 2\sqrt{p}, we can use the Chinese Remainder Theorem to deduce t once the product of the moduli exceeds 4\sqrt{p}.

For any specific curve E, once we have narrowed down the possibilities for t to a manageable number, we can use another method to arrive at the correct value of t faster. For details, consult Section 4.5.

3. Mathematical Background

3.1. A Few Unusual Curves. The theory is a lot cleaner if we assume that for any elliptic curve E that occurs later, j(E) is not 0 or 1728, and E is not supersingular (recall that this means #E(𝔽_{p}) = p + 1). The first condition assures that E does not have too many automorphisms, the second one assures that the endomorphism ring of E has rank 2 over ℤ (as opposed to rank 4 for supersingular curves).

Since each point on a supersingular curve is annihilated by p + 1 in the group law, whereas most points on most non-supersingular curves are not annihilated by p + 1, supersingularity can be effectively checked by multiplying a few points on the curve by p + 1. For details, see Section 2.5.

The curves just mentioned are in fact known to have a discrete logarithm problem that is easier than usual to solve, and hence their exclusion only does good to our intended applications in cryptography.

3.2. Reducing Point-Counting to Considering Isogenies. The idea of the algorithm is that since our elliptic curve E is defined over 𝔽_{p}, the Frobenius map F : (x, y) ↦ (x^{p}, y^{p}) acts on it. In fact, it satisfies

F^{2} − tF + p = 0.

The same relation is also valid for the action of F on the set of ℓ-torsion points of E, denoted by E[ℓ]. We could therefore obtain t′ = t mod ℓ by checking which 0 ≤ t′ ≤ ℓ − 1 satisfies (F^{2} + p)Q = t′F Q for all Q ∈ E[l]. This is the original idea of Schoof’s algorithm in [18].

Carrying out this plan quickly leads to arithmetic with polynomials of degree (ℓ^{2} − 1)/2 over 𝔽_{p}. Therefore, in practice it is better to use an idea of Elkies, which works only for half the primes ℓ, but then only requires the use of polynomials of degree (ℓ − 1)/2 over 𝔽_{p}.

The method works if t^{2} − 4p is a square modulo ℓ. This condition is satisfied for (asymptotically) half of all prime numbers ℓ. In this case, the action of F on E[ℓ] splits into (usually) two eigenspaces. If we somehow knew a subgroup B of E[ℓ] where F acts as a scalar, then we could check which 0 ≤ e ≤ ℓ − 1 has the property that

3Similar formulae exist for even s > 2, but those will not be needed in this algorithm.

eQ = F Q

for all Q ∈ B. Then we could use the congruence t′ ≡ e + p/e (mod ℓ) to get t′. Such a subgroup B can be described by a degree (ℓ − 1)/2 polynomial h(X) that vanishes exactly at the x coordinates of points in B. Indeed, given h(X) what we need to check is whether

e(X, Y) ≡ (X^{p}, Y^{p}) (mod h(X), Y^{2} − X^{3} − a_{4}X − a_{6}).

(Here the left hand side means (X, Y) ⊕ ⋯ ⊕ (X, Y) (e times), where ⊕ denotes addition according to the group law of E.)

Any subgroup B of size ℓ in E determines a degree ℓ isogeny (to be referred as an ℓ-isogeny) E → E/B. The group B being (in) an eigenspace of F is equivalent to this isogeny being defined over 𝔽_{p}. So we have reduced our problem to finding an explicitly described ℓ-isogeny E → E_{1} defined over 𝔽_{p}, if such a thing exists. (This happens exactly if t^{2} − 4p is a square modulo ℓ, but we do not know t^{2} − 4p yet!)

3.3. Parametrizing ℓ-isogenies. Let us think about ℓ-isogenies of elliptic curves defined over ℂ first. Section 3.10 shall indicate how these considerations can be made to apply to elliptic curves defined over 𝔽_{p}.

Any elliptic curve E over ℂ can be regarded as a quotient of ℂ by a rank 2 lattice L. Multiplying all elements of the lattice by a complex constant does not change the isomorphism class of the curve, so we may assume that the lattice L is spanned by 1 and τ ∈ ℂ, with τ in the upper half plane. Hence

E ≅ ℂ/(ℤ ⊕ τℤ).

However, for our purposes it is useful to choose another uniformization of E, that obtained via the map z ↦ exp(2πiz). It is clear that this gives an isomorphism

E ≅ ℂ^{∗}/q^{ℤ},

where q = exp(2πiz) and q^{ℤ} = {q^{n} : n ∈ ℤ}. Any compact Riemann surface (such as ℂ^{∗}/q^{ℤ}) is a projective algebraic variety (c.f. [11, B, Theorem 3.1]), and it turns out that ℂ^{∗}/q^{ℤ} is (isomorphic to) the elliptic curve

E : Y^{2} = X^{3} − \frac{E_{4}(q)}{48} + \frac{E_{6}(q)}{864},

where E_{4}(q) and E_{6}(q) are the q-series mentioned in Section 2.3. This E is in fact ℓ-isogenous to E_{1} = ℂ^{∗}/q^{ℓℤ}, via the map z ↦ z^{ℓ}. We can once again use q-series to describe the equation of E_{1}: indeed,

E_{1} : Y^{2} = X^{3} − \frac{E_{4}(q^{ℓ})}{48} + \frac{E_{6}(q^{ℓ})}{864}.

We also need the kernel polynomial h(X), which vanishes exactly at the x coordinates of points in the kernel of our isogeny (these points are just the images of the ℓth roots of 1 in ℂ^{∗}/q^{ℤ}). In fact, as explained in Section 3.9, it suffices to calculate here the first coefficient p_{1} of h(X), because the rest of the coefficients can be obtained from p_{1} and the data we already have. As for p_{1}, it turns out that

p_{1} = \frac{ℓ}{24} (E_{2}(q) − ℓE_{2}(q^{ℓ})) ,

where

E_{2}(q) = 1 − 24 ∑_{n=1}^{∞} σ_{1}(n)q^{n} = 1 − 24q − 72q^{2} − 96q^{3} − ⋯

Hence, if we are given an elliptic curve

E : Y^{2} = X^{3} + a_{4}X + a_{6},

we can find ℓ-isogenies by locating various q such that E ≅ ℂ^{∗}/q^{ℤ}. For each such q, we can find the equation of E_{1} ≅ ℂ^{∗}/q^{ℓℤ} as well as the kernel of the isogeny.

3.4. Locating the Parameters that Belong to a Certain Isogeny. The method suggested in Section 3.3 is impractical in the sense that infinitely many different choices of q give the same elliptic curve ℂ^{∗}/q^{ℤ} (up to isomorphism). Two elliptic curves are isomorphic if and only if they have the same j-invariant. The jinvariant of ℂ^{∗}/q^{ℤ} turns out to be the q-series j(q) = 1728E_{4}(q)^{3}/(E_{4}(q)^{3} −E_{6}(q)^{2}) (the same j(q) that occurs in Section 2.3). Any complex number can occur as the j-invariant of an elliptic curve, and the q-series j(q) can actually take any value as we vary q. Therefore, we could just pick a q for each possible j-invariant and restrict our attention to just those q, since all other ones would just clutter our calculations without yielding new elliptic curves.

In fact, we can do a little better. Instead of restricting attention to those q that give us all the isomorphism classes of elliptic curves E as ℂ^{∗}/q^{ℤ}, we take a bigger set of q so that the corresponding ℓ-isogenies ℂ^{∗}/q^{ℤ} → ℂ^{∗}/q^{ℓℤ} give us all isomorphism classes of ℓ-isogenies E → E_{1}. It turns out this set can be given the structure of an algebraic curve called Y_{0}(ℓ). The addition of two extra points (called the cusps) compactifies the curve, and we obtain a smooth projective curve X_{0}(ℓ). The curve X_{0}(ℓ) is a modular curve, so called since its points are moduli (parameters) for ℓ-isogenies of elliptic curves.

The fact that X_{0}(ℓ) has geometric structure is a powerful tool. For example, the q-series j(q) is just a (rational) function on X_{0}(ℓ), for each point giving us the jinvariant of the “source” curve of the isogeny parametrized by that point. Similarly, j(q^{l}) is a rational function on X_{0}(ℓ) that gives the the j invariant of the target. Also, given an ℓ-isogeny ϕ : E → E_{1}, we can take the dual isogeny ϕ^{∨} : E_{1} → E. This operation is reflected on X_{0}(ℓ) by a morphism w : X_{0}(ℓ) → X_{0}(ℓ). If Q is the point that parametrizes ϕ, then wQ parametrizes ϕ^{∨}. Since (ϕ^{∨})^{∨} = ϕ, we conclude that the morphism w ∘ w is just the identity. The morphism w is called the Atkin-Lehner involution. (It swaps the two cusps.)

Now we can get rid of another unwieldy feature of the process suggested in Section 3.3: that of determining q from E_{4}(q) and E_{6}(q) and using this value to evaluate other q-series. That would be a very tedious task which can be sidestepped as follows. The functions j(q) and j(q^{ℓ}) are known to satisfy a polynomial relation, the Kronecker relation:

Φ(j(q), j(q^{l})) = 0,

for a certain Φ = Φ_{ℓ} ∈ ℤ[X, Y] that can be determined explicitly. The polynomial Φ has degree ℓ + 1 in each variable.

To explain how to use the Kronecker relation, we need some notation. For any q-series g, denote the q-series g(q^{l}) by g̃. Furthermore, for any q-series g(q), let g′(q) = q(dg/dq). In fact g′(q) = (1/2πi)(dg/dτ) (recall that q = exp(2πiτ)), so this g′(q) can be used as a usual derivative of g. When both ′ and ˜ are involved, the ˜ is to be taken last, so g̃′ is intended to mean ˜(g′). This is an important point, since (g̃)′ = ℓ˜(g′).

The q-series j, E_{4} and E_{6} enjoy a number of identities, for example

\frac{j′}{j} = −\frac{E_{6}}{E_{4}} , \frac{j′}{j − 1728} = −\frac{E_{4}^{2}}{E_{6}} .

Therefore, if we know the values of j(q), E_{4}(q) and E_{6}(q) (but not that of q), we can still find the value of j′(q).

Now we can put the Kronecker relation to good use. Assume that we know the values of E_{4}(q) and E_{6}(q) and we want to describe the isogeny ℂ^{∗}/q^{ℤ} → ℂ^{∗}/q^{ℓℤ}. We can immediately determine j(q) = 1728E_{4}(q)^{3}/(E_{4}(q)^{3} − E_{6}(q)^{2}), and by the identities above, also the value of j′(q). Plugging the value of j(q) into the Kronecker relation, we obtain a polynomial whose ℓ + 1 roots are the possible values of j(q^{ℓ}) = j̃. Differentiating the Kronecker relation (using the notation for partial derivatives introduced in Section 2.6), we obtain

j′Φ_{1}(j, j̃) + ℓj̃′Φ_{2}(j, j̃) = 0.

For each candidate for j̃, this allows us to deduce the value of j̃′, and hence that of Ẽ_{4} and Ẽ_{6} by the identities above (which of course remain valid if we plug in q^{ℓ} instead of q). We can also get the value of E_{2} − ℓẼ_{2} for p_{1} by looking at the second derivative of the Kronecker relation and then utilizing some more complicated identities to get at E_{2}, but we won’t do that since there is an even better way. The aspect of this method that needs improving is its use of Φ_{ℓ}: that equation has extremely large coefficients (which would actually be fine with us since we are only trying to compute it modulo p in the end), and determining it requires that we calculate using approximately 2ℓ^{2} coefficients of a certain power series. The improved method cuts that 2ℓ^{2} to about 2ℓv, where v is less than ℓ/24 + 1.

What one should do instead is to choose a suitable q-series f that is a function on X_{0}(ℓ). We also require that f be invariant under the Atkin-Lehner involution w. For technical reasons (that will be explained in Section 3.5), we pick f to have poles only at the cusps, and we want those poles to have relatively low order. Then we find by the process of Section 3.5 a polynomial relation

Ψ_{ℓ}(f, j) = 0.

By the w-invariance of f, we also have

Ψ_{ℓ}(f, j̃) = 0.

Hence given j, we can find all possible values of f, and for each possible value of f, we can get all values of j̃. We can take partial derivatives of the above two relations and gradually deduce all the information we need to describe our isogeny.

It is possible for distinct pairs (j_{1}, j̃_{1}) and (j_{2}, j̃_{2}) to correspond to the same value of f. Therefore, starting from j_{1}, we might get this f and end up with j̃_{2}. However, the p_{1}, p_{2}, … of Section 2.8 obtained from (j_{1}, f, j̃_{2}) will not be power sums of the roots of a degree d polynomial, and hence will be deemed invalid because s_{d+i} ≠ 0 for some small i. Readers not satisfied with this test are referred to Section 3.11.

The calculations needed are explained in Morain’s [15, 3.2.2]. (The second formula on page 273 contains a typo: a denominator should have a 6 instead of a 2. However, it is evident that the correct version was used as input for the subsequent calculations. Also, the numerical example on top of page 274 is of course for ℓ = 31, not 23.) Morain’s notation is similar to ours, the following conversion rules apply:

Morain here  Morain  here  Morain here
F  f  ∂_{F F} Ψ_{11}(f, j) ∂_{F F ∗} Ψ_{11}(f, j̃)
J  j  ∂_{F J} Ψ_{12}(f, j) ∂_{F J∗} Ψ_{12}(f, j̃)
J̃  j̃  ∂_{JJ} Ψ_{22}(f, j) ∂_{JJ∗} Ψ_{22}(f, j̃)
∂_{F}  Ψ_{1}(f, j) ∂_{F ∗}  Ψ_{1}(f, j̃)
∂_{J}  Ψ_{2}(f, j)  ∂_{J∗}  Ψ_{2}(f, j̃)

The polynomials Φ_{ℓ} and Ψ_{ℓ} are not the only ones that can be used here (although Ψ_{ℓ} is the optimal one in a certain sense). For an analysis of these and related polynomials and algorithms for computing them, see [2].

3.5. Obtaining Ψ_{ℓ} from f. Using the terminology of Section 3.4, we introduce another modular curve. It is X(1), and it parametrizes isomorphism classes of elliptic curves. As a Riemann surface, X(1) is just the complex numbers plus a cusp that is considered to lie “at infinity”: each finite number corresponds to the elliptic curve with that j-invariant. (The cusp does not parametrize an elliptic curve, it is just there to make our object compact and geometrically nice.) If we agree that j = ∞ at the cusp, then j gives an isomorphism between X(1) and the complex projective line.

There is a “degeneracy” morphism ϕ : X_{0}(ℓ) → X(1). It sends the point parametrizing E → E_{1} to the point parametrizing E (and sends both cusps of X_{0}(ℓ) to the unique cusp of X(1)). The map ϕ is a ramified Galois covering. Using the coordinate τ on X_{0}(ℓ) (recall that q = exp(2πiτ)), for any value of τ, the set {τ, −1/τ, −1/(τ + 1), …, −1/(τ + ℓ − 1)} is invariant under the deck transformation group of ϕ. So if g is any function on X_{0}(ℓ), then any elementary symmetric polynomial s in {g(τ), g(−1/τ), g(−1/(τ + 1)), …, g(−1/(τ + ℓ − 1))} is actually a pullback of a function on X(1). Therefore s “is” a rational function in j.

Since wτ = −1/ℓτ, if f is a w-invariant function on X_{0}(ℓ) then all the coefficients of the polynomial

Ψ(F) = (F − f(τ)) ∏_{k=0}^{ℓ−1} (F − f (\frac{−1}{τ + k} ))

= (F − f(τ)) ∏_{k=0}^{ℓ−1} (F − f (\frac{τ + k}{ℓ} ))

are rational functions in j(τ). Additionally, if f has poles only at the cusps, then the coefficients of Ψ(F) must be polynomials in j(q). Hence if we calculate the q-expansions of the coefficients of Ψ(F), we can read off what those polynomials in j(q) must be. Thereby we obtain a polynomial relation

Ψ_{ℓ}(F, J) = 0

which is satisfied by F = f(q) and J = j(q).

It remains to calculate the q-expansions of the coefficients of Ψ(F), given the q-expansion of f. (The question of picking a suitable f with a known q-expansion will be addressed in Section 3.6.) Let ζ = exp(2πi/ℓ) and

f^{r}(q) = ∑_{n} a(n, r)q^{n}.

Then

f (\frac{τ + k}{ℓ} )^{r} = f^{r}(q^{1/ℓ}ζ^{k}) = ∑_{n} a(n, r)q^{n/ℓ}ζ^{kn}

Summing the above from k = 0 to k = ℓ − 1, the roots of unity cause the fractional powers of q to cancel each other out, and we are left with

s_{r}(q) = ∑_{k=0}^{ℓ−1} f (\frac{τ + k}{ℓ} )^{r} = ℓ ∑_{n} a(ℓn, r)q^{n}.

These s_{r}(q) are exactly the power sums of the roots of the polynomial ∏_{k=0}^{ℓ−1} (F − f (\frac{τ+k}{ℓ} )). From there we can obtain its coefficients c_{r}(q) in the usual manner, from which the coefficients C_{r}(q) of Ψ(F) come out. The details are given in Section 2.4.

3.6. Finding a good f(q) using hands-on methods. Let X_{0}^{+}(ℓ) denote the smooth curve that is the quotient of X_{0}(ℓ) by the Atkin-Lehner involution w. The requirements of Section 3.4 and 3.5 then translate simply to finding (a lot of initial coefficients of) the q-expansion of a certain non-constant function f on X_{0}^{+}(ℓ). This f is required to have a pole only at the (unique) cusp of X_{0}^{+}(ℓ), and the order v of this pole should be as small as possible. Denote the genus of X_{0}^{+}(ℓ) by g_{0}^{+}(ℓ), and the genus of X_{0}(ℓ) by g_{0}(ℓ). A simple application of the Riemann–Roch theorem tells us that we can always find such a function f with v ≤ g_{0}^{+}(ℓ) + 1. By [15, Thm. 2.1] and the fact that g_{0}(ℓ) is about ℓ/12, we confirm that we can always get a v no bigger than about ℓ/24. (For each ℓ, the geometric properties of X_{0}^{+}(ℓ) determine what the smallest v is. The definition of the cusp being a Weierstrass point of the curve is that v < g_{0}^{+}(ℓ) + 1 can be achieved.

In fact, if ℓ is very small then finding an f is straightforward. In this case Ψ_{ℓ} also has very small coefficients, so we may as well just determine it from the f in advance, and hardwire the polynomial Ψ_{ℓ} into our program.

For instance, if ℓ ∈ {3, 5, 7, 13}, the curve X_{0}(ℓ) is rational, isomorphic to the projective line via

h(q) = (\frac{η(q)}{η(q^{ℓ})})^{24/(ℓ−1)} .

(Here we use η(q) = q^{1/24}η̄(q), c.f. [9, The rational case in Section 4].) This h has a single pole at q = 0. Since (w^{∗}h)h = l^{12/(ℓ−1)}, we conclude that

f(q) = h(q) + ℓ^{12/(ℓ−1)}/h(q)

is a degree-1 function on X_{0}(ℓ) that is invariant under the action of w. Therefore it is a function on the rational curve X_{0}^{+}(ℓ), with a single pole at the cusp. The qexpansion of f is easily calculated to any reasonable precision by using the explicit definition of h(q). For example, for ℓ = 3 we get

f(q) = q^{−1} − 12 + 783q + 8672q^{2} + 65367q^{3} + ⋯

For ℓ ∈ {11, 17, 19, 23}, we use ad hoc methods to produce our function f.

The case ℓ = 11 is treated explicitly in [15, Section 2.3.1]. A suitable f is given as

f(q) = (\frac{θ([2, 1; 1, 6], q)}{η(q)η(q^{11})})^{2} − 1 = q^{−1} + 5 + 17q + 46q^{2} + ⋯

(For information on θ-series and our notation, refer to [14, Section 4.9, especially Corollary 4.9.5.(3)].)

Now consider the case ℓ = 17. First we write down the q-series

f_{1}(q) = (17E_{2}(q^{17}) − E_{2}(q))/24,
f_{2}(q) = 17^{2}E_{4}(q^{17}) + E_{4}(q),
f_{3}(q) = θ([12, 1, 0, 3; 1, 2, 2, 1; 0, 2, 8, 3; 3, 1, 3, 4], q).

All three of these functions are modular forms with respect to the group Γ_{0}(17). The forms f_{1}, f_{2}, f_{3} have weight 2, 4, 2 respectively. The functions f_{1}, f_{3} are antiinvariant under the action of the Atkin-Lehner involution, whereas f_{2} is invariant. Therefore the q-series

a(q) = f_{3}^{2}/f_{2} = 1/290 + (46/4205)q + ⋯,
b(q) = f_{1}/f_{3} = 2/3 − (1/3)q + ⋯

are functions on X_{0}^{+}(17). They satisfy the equation

(b^{3} − \frac{19}{3}b^{2} + \frac{15}{4}b − \frac{1}{4})a + (\frac{1}{360}b − \frac{1}{1080}) = 0.

The above equation can be considered a singular equation of the curve X_{0}^{+}(17). Therefore its normalization (desingularization) is isomorphic to X_{0}^{+}(17). The given equation is only singular at the point at infinity in the direction b = 0, where it has an ordinary triple point. By blowing up the singularity, we see that the normalization is a space curve that lies on a nonsingular the quadric surface in ℙ^{3}. The normalization itself is a (1,3)-curve on the surface. (For details of the blowing up process and the naming convention for curves on the nonsingular quadric surface, refer to [11, Section I.4 and (II, Ex. 5.6)].) In the end it turns out that b gives an isomorphism between X_{0}^{+}(17) and the projective line. The cusp q = 0 corresponds to b = 2/3. Since we want a pole at the cusp, we set

f = \frac{−1/3}{b − 2/3} = q^{−1} + 3 + 7q + 14q^{2} + 29q^{3} + ⋯

The cases ℓ = 19 and ℓ = 23 may be handled in a similar way.

3.7. Finding a good f(q) by lifting it from characteristic ℓ. The difficulty in finding a good f is essentially that we have to know about some complicated geometric properties of the curve X_{0}^{+}(ℓ). An idea of Atkin is very useful in this context. Instead of trying to find the coefficients in the q-expansion of f, we try to find them modulo some prime number and use those values to guess the original coefficients in f. With any luck, the reductions of the coefficients of f will have to do with a simpler curve and will therefore be simpler to calculate.

We shall reduce X_{0}^{+}(ℓ) modulo ℓ. Formally, this amounts to taking the model X_{0}^{+}(ℓ) that Deligne and Rapoport define over ℤ_{ℓ}, and considering its special fiber X_{0}^{+}(ℓ)_{/𝔽_{ℓ}}. The curve thus obtained is not smooth, but it is very simple: it is isomorphic to the projective j-line over 𝔽_{ℓ} with a certain number of simple double points. (For details, see [6, Thm. 6.9].)

The double points result from the identification of the points x and x^{ℓ} for each of x ∈ S, for a certain set S. To define S, we need to talk about supersingular j-invariants. Let K be any field of characteristic ℓ. In this more general setting, there is still a notion of whether an elliptic curve defined over K is ordinary or supersingular. For example, we may say that a curve E is supersingular if E[ℓ] = 0, otherwise it is ordinary. This is consistent with our previous definition of supersingular in the case where K = 𝔽_{ℓ} (for ℓ ≥ 5), since for E defined over such an 𝔽_{ℓ}, E[ℓ] = 0 occurs if and only if^{4} #E(𝔽_{ℓ}) = ℓ + 1. (Recall that the notation E[ℓ] stands for the points of order ℓ on E defined over the algebraic closure of the base field.) It is an interesting fact that for any field K of characteristic ℓ and any elliptic curve E defined over that field, the curve is supersingular if and only if j(E) ∈ T, for a certain finite set T ⊆ 𝔽_{ℓ^{2}}. Then S is just the set of those elements of T that are not in 𝔽_{ℓ}. We shall discuss how to find the set S explicitly in Section 3.8.

Therefore, a function f on X_{0}^{+}(ℓ)_{/𝔽_{ℓ}} that has its only pole at the cusp is nothing but a polynomial function P (j) in j with the property that P (x) = P (x^{ℓ}) for each x ∈ S.

Assume now that we have already found a suitable function f on X_{0}^{+}(ℓ). Then its reduction modulo ℓ is just a function P (j) on X_{0}^{+}(ℓ)_{/𝔽_{ℓ}}. From what was described above (and with our explicit knowledge of the set S), we can determine the nonconstant P of lowest degree v that can occur in this way. Then we could just evaluate

P (j(q)) mod ℓ

and lift the coefficients back to characteristic 0 to obtain our original f. So far, so good, but there are two difficulties that we need to address. Firstly, we need some explicit bounds on the coefficients of f to make sure that

our lifting is unique. However, we expect the coefficients of f to increase exponentially, and therefore there is no chance that we could lift the (at least) hundreds of coefficients that we need. Atkin’s laundering process solves this problem. If (ℓ + 1)/24 > v, then we can look at the function

g(q) = f(q)η(q)η(q^{ℓ}).

It is a cusp form of weight 1, level ℓ (and some character), and therefore its nth coefficient is only O(n^{ε}). Furthermore, if g(q) is also a newform then we have an explicit bound: the nth coefficent is no more than the sum of divisors of n. In fact, if g(q) is a linear combination (with small coefficients) of newforms, then we still have a good explicit bound on the size of the coefficients, enabling us to reconstruct a lot of them from their mod ℓ reductions, which we obtain by calculating

P (j(q))η(q)η(q^{ℓ}) mod ℓ.

A suggestion of Atkin on how to make g(q) a linear combination of newforms with small coefficients is to choose P in such a way that the coefficient of q^{(ℓ+1)/24} in P (j(q))η(q)η(q^{ℓ}) mod ℓ is 0. This amounts to modifying the constant term of P (X).

A second problem has to do with the possibility that the q-expansion of f might start with some terms that are divisible by ℓ and therefore would be visible in the mod ℓ picture. This does not seem to occur for any of the ℓs in use. Atkin, Elkies and yours truly have calculated a number of examples and confirmed that this phenomenon did not occur there.

Of course, once we have a lot of coefficients of g(q), we can find a lot of coefficients of

f(q) = g(q)/(η(q)η(q^{ℓ})).

4For ℓ = 2, 3, we might also have #E(𝔽_{ℓ}) = 1 or #E(𝔽_{ℓ}) = 2ℓ + 1.

The above process works only if the condition (ℓ + 1)/24 > v is satisfied. This is the case if ℓ ≥ 29, with the exceptions ℓ = 37, 43, 67, 163.

Finally, we remark that a suitable g can also be found either as a combination of suitably chosen theta-series (see [1]) or by considering the action of a certain Hecke operator on the modular form η(q)η(q^{ℓ}) (as in [16]).

3.8. Finding functions on X_{0}^{+}(ℓ)_{/𝔽_{ℓ}}. Let H_{ℓ}(X) ∈ 𝔽_{ℓ}[X] denote the polynomial whose roots are exactly the set T of supersingular j-invariants in characteristic ℓ. This polynomial is easily calculated in a number of ways, for example by using [15, Thm. 2.2].

Since all elements of T are contained in 𝔽_{ℓ^{2}}, the polynomial H_{ℓ}(X) splits completely into linear and quadratic factors. Let H_{ℓ}^{∗}(X) denote the product of the quadratic factors (which is easily calculated once H_{ℓ}(X) is known). The roots of H_{ℓ}^{∗}(X) are exactly the elements of the set S mentioned in Section 3.7.

We are aiming to find a non-constant polynomial P ∈ 𝔽_{ℓ}[X] of the lowest possible degree that satisfies

P (x) = P (x^{ℓ})

for all x ∈ S. This is equivalent to saying that

P (X^{ℓ}) − P (X) ≡ 0 (mod H_{ℓ}^{∗}(X)).

Indeed, this is the same as saying that P (X) is in the kernel of the Berlekamp matrix associated to H_{ℓ}^{∗}(X). Therefore, we merely need to find a vanishing linear combination of the quantities X^{ℓ} − X mod H_{ℓ}^{∗}(X), X^{2ℓ} − X^{2} mod H_{ℓ}^{∗}(X), X^{3ℓ} − X^{3} mod H_{ℓ}^{∗}(X), …, and the corresponding linear combination of X, X^{2}, X^{3}, … will give us P (X).

Finally, we also need to adjust the constant term of P (X) as indicated in Section 3.7.

3.9. Obtaining the rest of the coefficients of h(X). In [9], in the subsection titled The kernel of the isogeny in Section 3, Elkies carefully explains how plugging the q-expansions for the x and y coordinates on E_{1} = ℂ^{∗}/q^{ℓℤ} into the derivative of the Weierstrass equation of that curve gives a recursion on the coefficients of h(X) that enables one to find all coefficients once a_{4}, a_{6}, p_{1}, ã_{4}, ã_{6} are known. The actual formulae deduced from there are listed in Section 2.7.

3.10. Isogenies of elliptic curves over 𝔽_{p}. Although in the previous sections we worked over ℂ, exactly the same formulae work over 𝔽_{p}. We can justify this using Deuring’s lifting theorem: for any isogeny E → E_{1} over 𝔽_{p}, we can just lift the coefficients to a number field and reduce them again to when we worked out all the data we needed. For details, consult [19, §7]. We could also appeal to the algebraic definition of modular forms as sections of certain sheaves and arrive at their q-expansions that way.

Secondly, for an ordinary elliptic curve defined over 𝔽_{p} with j-invariant not equal to 0 or 1728, an isogeny is defined over 𝔽_{p} if and only if j, f, j̃ ∈ 𝔽_{p} (see [19, Prop. 6.1a]).

3.11. Checking the correctness of h(X). Note that once a proposed h(X) is determined, it is easy for the suspicious user to check (without using any methods based on modular curves) that it has the claimed properties.

Firstly, calculating ℓ(X, Y) (mod h(X), Y^{2} − a_{4}X − a_{6}), we can confirm that it is neutral element of the elliptic curve, and thereby show that h(X) is indeed the polynomial associated to some order ℓ subgroup B of E. Secondly, the existence of some e such that e(X, Y) ≡ (X^{p}, Y^{p}) (mod h(X), Y^{2} − a_{4}X − a_{6}) shows that this subgroup B is invariant under the action of the Frobenius morphism.

3.12. Determining e given h(X). We have now described a subgroup B of E[ℓ] where Frobenius acts by multiplication by e. Determining the e for which

e(X, Y) ≡ (X^{p}, Y^{p}) (mod h(X), Y^{2} − X^{3} − a_{4}X − a_{6})

amounts to determining e such that eQ = F Q for all Q ∈ B. But since B is cyclic of prime order, it suffices to check eQ = F Q for any one nonzero Q ∈ B. Therefore we may replace h by any nontrivial factor h′ of h in our calculations, since this just amounts to checking eQ = F Q for just those (at least one nonzero) Q ∈ B where h′ vanishes on the x-coordinate.

It is evident that h(X) has no double roots. It is also true that h(X) will split into factors each of which have degree s(e), where s(x), the semi-order of x, is the smallest positive integer n such that x^{n} ≡ ±1 (mod ℓ). Indeed, let Q = (x_{0}, y_{0}) be any particular element of B. Then e(x_{0}, y_{0}) = (x_{0}^{p}, y_{0}^{p}) and therefore (x_{m}, y_{m}) := e^{m}(x_{0}, y_{0}) = (x_{0}^{p^{m}}, y_{0}^{p^{m}}). Therefore x_{0}^{p^{m}} = x_{m} = x_{0} if and only if e^{m} ≡ ±1 (mod ℓ). By definition, the smallest m for which this occurs is s(e), so the minimal polynomial of x_{0} has degree s(e) as claimed.

To calculate the e for which e(X, Y) ≡ (X^{p}, Y^{p}) holds, we just take some number s with the property that the coordinates sP are easy to calculate from the coordinates of P, for any point P on the curve. For example, for s = 2 we have the duplication formulas. This s will also have to have the property that s^{k} hits all possible residue classes modulo ℓ. Actually, since P and −P have the same x coordinates and opposite y coordinates, so it is sufficient for s^{k} to hit either a residue class or its negative. A number with such a property may be called a semiprimitive root for ℓ. Determining a good small s for each ℓ is trivial. For this algorithm, s ≤ 11 will suffice (see Section 5.2).

As for the calculation, observe that all can be done in the ring 𝔽_{p}[X]/(h(X)). We represent each point that occurs as (P_{1}(X), Y P_{2}(X)). The appropriate formulas for taking s times a point are easily derived from the formulas on page 105 in [20].

Note that finding e (given h(X)) can be carried out faster by using a baby step– giant step or similar method, with the arithmetic carried out as in [3, Chapter IV].

3.13. Dewaghe’s improvement. L. Dewaghe proposed a method that makes finding e significantly faster when deg(h) is odd. This condition is bound to occur when ℓ ≡ 3 (mod 4), since deg(h) divides (ℓ − 1)/2 which is odd in this case. Of course deg(h) might happen to be odd also when ℓ ≡ 1 (mod 4).

According to [7, Theorem 1 of Section 3.2], if deg(h) is odd and the x-coordinates of λ(X, Y) and (X^{p}, Y^{p}) match, then

e = λ^{s(λ)} (\frac{r}{p}) λ.

Here r is the resultant of h(X) and w(X) = X^{3} + a_{4}X + a_{6}.

3.14. The accuracy to which f(q) needs to be calculated. We will demonstrate that it is sufficient to evaluate the coefficients of f(q) up to and including that of q^{2ℓv−v}.

Let g(q) be any q-series that is only evaluated to finite precision. Let v(g) be the order of the pole at q = 0 and let t(g) be the smallest integer such that the coefficient of q^{t(g)} is not known. Mnemonically,

g(q) = cq^{−v(g)} + ⋯ + O(q^{t(g)}),

where c is a nonzero constant. We know that v(f) = v. Assume now that t(f) = 2ℓv − v + 1. Using the notation of Section 2.4, we shall prove that t(C_{k}) > 0 for 0 ≤ k ≤ ℓ + 1, and therefore that Ψ_{ℓ} can indeed be determined from the coefficients of f that are known.

For any 1 ≤ k ≤ ℓ, evidently v(f^{k}) = kv(f) = kv, and

t(f^{k}) = 2ℓv − v + 1 − (k − 1)v = v(2ℓ − k) + 1.

Let ⌊x⌋ denote the largest integer not exceeding x. By the definition of s_{k}, for all 1 ≤ k ≤ ℓ, we conclude that

v(s_{k}) ≤ ⌊kv/ℓ⌋
t(s_{k}) = ⌊v(2ℓ − k)/ℓ⌋ + 1.

The inequality sign above reflects the fact that a “randomly picked” coefficient in g^{k} might happen to be zero, in which case s_{k} might have a pole of lower order than might have been supposed at first.

Now looking at the definition of c_{k} and using the elementary fact that for any pair a, b of integers, we always have ⌊a/ℓ⌋ + ⌊b/ℓ⌋ ≤ ⌊(a + b)/ℓ⌋, we conclude that

v(c_{k}) ≤ ⌊kv/ℓ⌋.

(This can be proved by induction on k, starting with k = 1.)

We can also show by induction that

t(c_{k}) ≥ t(s_{k}) = ⌊v(2ℓ − k)/ℓ⌋ + 1.

Indeed, the base case k = 1 is obvious. Then by the same elementary fact,

t(c_{k−r}s_{r}) = min(t(c_{k−r}) − v(s_{r}), t(s_{r}) − v(c_{k−r})) ≥ t(s_{k}),

and the rest is clear.

Looking at the definition of C_{k}, we see that if suffices to show that for all 1 ≤ k ≤ ℓ,

t(fc_{k}) = min(t(f) − v(c_{k}), t(c_{k}) − v(f)) ≥ 1.

This is equivalent to two inequalities, the first one of which is immediate. The second one comes down to

⌊v(2ℓ − k)/ℓ⌋ + 1 ≥ v + 1,

which follows from k ≤ ℓ.

4. Possible Improvements

4.1. Factoring h(X). The following considerations can be used to speed up the determination of e given h(X). As mentioned in Section 2.9, we might replace h(X) by any of its factors h′(X). This has the advantage that the lower the degree of h′(X) we use, the faster the process in the rest of Section 2.9 runs. Dewaghe’s method gives an additional boost of speed if deg(h′) is odd.

However, factoring h(X) to obtain h′(X) might take a long time. Therefore, an alternative strategy might be to not factor h(X), set h′(X) = h(X), and hope that we gain more time on factoring than we lose later.

In a test implementation, these two methods were about equally fast. Here is a method that is faster than either of the above. Recall that h(X) splits into irreducible factors of equal degree. For each ℓ, choose integers o_{ℓ} and e_{ℓ}.

The method is this. For each odd divisor N ≤ o_{ℓ} of deg(h) = (ℓ − 1)/2, check if h(X) has a factor h′(X) of degree N. If not, for each even divisor N ≤ e_{ℓ} of deg(h) = (ℓ − 1)/2, check if h(X) has a factor h′(X) of degree N. If still no divisor h′ is found, set h′ = h.

The constants o_{ℓ} and e_{ℓ} need to be chosen based on how long the particular implementation at hand can find factors of degree N of h(X) for various N, and the speed the rest of Section 2.9 takes as a function of the degree of h.

4.2. Doing with s ≤ 5. A simple calculation shows that for s ≥ 3 odd,

deg(a_{s}) = s^{2},
deg(b_{s}) = s^{2} − 1,
deg(c_{s}) = 3(s^{2} − 1)/2,
deg(d_{s}) = 3(s^{2} − 1)/2.

Hence carrying out the step in Section 2.9 will require plugging polynomials into other polynomials of degree up to 180. (Degree 180 occurs for ℓ = 109, s = 11.) However, we can get along with using s ≤ 5, in which case we have to plug into polynomials of degree no more than 36. This is because whenever s > 5 occurs (ℓ = 41, 109), it turns out that 6 is a semiprimitive root modulo ℓ. We could use s = 6 here, but also we can just break down multiplication by 6 to a multiplication by 2 followed by a multiplication by 3, and therefore we can always restrict ourselves to s ∈ {2, 3, 5}.

4.3. Isogeny cycles. If ℓ is a good auxiliary prime, we can often continue by finding a chain of isogenies from E to E_{1} to E_{2} to … to E_{k} etc. Then we can get a degree ℓ^{k−1}(ℓ − 1)/2 kernel polynomial in basically no time and we can determine t modulo ℓ^{k}. For example, when this works for ℓ = 11, we can find the eigenvalue modulo 121, which takes relatively long (about as long as for ℓ = 127), but at least Ψ_{11} does not have to be calculated. Note that here the arguments in Section 3.12 about factoring h(X) do not apply. Therefore we must either not factor h(X) or be very careful.

4.4. Atkin primes. The algorithm can be made to use slightly fewer auxiliary primes (resulting in a reduction of both running time and memory use) by using an idea of Atkin. Even if Ψ_{ℓ}(F, j) has no roots in 𝔽_{p}, we can get some restrictions on t^{2} mod ℓ by looking at the factorization of Ψ_{ℓ}(F, j). It will usually split into factors which all have the same degree r. The knowledge of r will usually eliminate about half the possibilities for t mod ℓ, but can even determine ±t mod ℓ if r is less than 5 (which is not that common). For details, consult [1] and [19, Proposition 6.2].

4.5. Catching kangaroos. Running our algorithm for many (but not all) auxiliary primes, we arrive at a stage where there are only a relatively small number of possible t (and they form an arithmetic sequence). Since #E(𝔽_{p}) = p + 1 − t, for the correct t the number p + 1 − t annihilates all points on E (under the group law of E). A t for which p + 1 − t annihilates several random points on E is then probably the correct one. (We use the fact that E(𝔽_{p}) “tends” to be cyclic for more certainty, we could also check the quadratic twist E′ of E, which has p + 1 + t points. For details, see [19, Section 3, especially Theorem 3.2].)

For a random point P on E, we can use either Shanks’ baby step–giant step method or Pollard’s lambda (kangaroo trap) method to determine for which of the possible ts will p + 1 − t annihilate P.

A more significant reduction in running time would result from not dropping quite as many auxiliary primes as our kangaroo trap would allow. Then the probability of the algorithm succeeding would go up, and therefore we wouldn’t have to try quite as many curves to get one with the properties described in Section 1.

5. Tables

5.1. Table of the P_{ℓ}(J). Note that the primes above 197 do not appear in the sets A_{s} or A_{l}. They are provided here for the reader’s convenience in implementing the algorithm for use with primes p with more than 200 bits.

ℓ  P_{ℓ}(J)
29  J + 11
31  J + 1
41  J − 5
47  J + 9
53  J^{2} − 3J + 26
59  J + 24
61  J^{2} − 23J − 1
71  J − 33
73  J^{3} + 32J^{2} − 30J + 1
79  J^{2} + 14J − 1
83  J^{2} + 7J − 2
89  J^{2} + 26J − 17
97  J^{4} + 32J^{3} + 42J^{2} − 24J − 2
101  J^{2} + 27J − 13
103  J^{3} + 34J^{2} − 7J − 2
107  J^{3} + 16J^{2} − 32J + 11
109  J^{3} − 51J^{2} + 52J
113  J^{4} − 37J^{3} + 24J^{2} − 3J − 36
127  J^{4} − 54J^{3} − 41J^{2} − 32J − 2
131  J^{2} − 47J − 51
137  J^{5} − 20J^{4} − 23J^{3} + 53J^{2} + 65J + 52
139  J^{4} − 56J^{3} − 18J^{2} + 40J + 1
149  J^{4} + 5J^{3} − 61J^{2} + 48J − 57
151  J^{3} + 34J^{2} − 7J − 1

ℓ  P_{ℓ}(J)^{157  J^{6} − 67J^{5} − 48J^{4} − 70J^{3} + 30J^{2} − 5J + 2}^{167  J^{3} − 60J^{2} + 3J − 14}^{173  J^{4} − 34J^{3} − 60J^{2} − 74J − 22}^{179  J^{3} − 83J^{2} + 18J − 62}^{181  J^{6} + 62J^{5} + 12J^{4} + 82J^{3} − 10J^{2} + 51J − 1}^{191  J^{3} + 60J^{2} − 25J + 56}^{193  J^{6} − 24J^{5} − 15J^{4} − 70J^{3} − 53J^{2} + 58J + 1}^{197  J^{6} + 68J^{5} + 14J^{4} − 66J^{3} + 3J^{2} − 59J + 59}^{199  J^{5} + 62J^{4} + 19J^{3} + 46J^{2} + 87J + 1}^{211  J^{7} + 68J^{6} + 88J^{5} + 77J^{4} + 67J^{3} − 10J^{2} − 82J + 2}^{223  J^{7} − 78J^{6} + 108J^{5} − 41J^{4} + 92J^{3} + 78J^{2} − 17J + 1}^{227  J^{5} − 88J^{4} + 45J^{3} + 62J^{2} + 84J − 73}^{229  J^{7} + 60J^{6} + 39J^{5} − 37J^{4} − 114J^{3} − 42J^{2} + 59J + 1}^{233  J^{7} − 81J^{6} + 91J^{5} + 30J^{4} + 115J^{3} + 105J^{2} + 96J + 30}^{239  J^{4} − 107J^{3} − 39J^{2} − 9J − 28}^{241  J^{6} + 115J^{5} − 54J^{4} − 105J^{3} + 40J^{2} + 85J + 1}^{251  J^{4} + 36J^{3} + 87J^{2} − 8J + 62}^{257  J^{6} − 94J^{5} − 112J^{4} − 51J^{3} − 124J^{2} + 42J − 122}^{263  J^{5} − 37J^{4} − 94J^{3} − 40J^{2} − 78J + 74}^{269  J^{6} + 109J^{5} − 38J^{4} − 7J^{3} + 62J^{2} + 4J − 120}^{271  J^{6} − 127J^{5} − 69J^{4} + 132J^{3} + 46J^{2} − 63J − 2}^{277 J^{9} − 48J^{8} − 70J^{7} + 3J^{6} − 87J^{5} + 25J^{4} + 105J^{3} + 39J^{2} + 24J + 1}^{281  J^{7} + 131J^{6} − 34J^{5} − 36J^{4} + 101J^{3} + 114J^{2} + 115J + 53}

5.2. Table of the s for each ℓ ∈ A.

\begin{matrix}ℓ & s \\ 3, 5, 7, 11, 13, 19, 23, 29, 47, 53, 59, 61, 71, 79, 83, 101 & 2 \\ 103, 107, 131, 139, 149, 167, 173, 179, 181, 191, 197 & 2 \\ 17, 31, 89, 113, 127, 137 & 3 \\ 73, 97, 151, 157, 193 & 5 \\ 41 & 7 \\ 109 & 11 \end{matrix}

5.3. Table of the Ψ_{ℓ} for each ℓ ∈ A_{s}.

Ψ_{3}(F, J) = F^{4} + (−J + 792)F^{3} + (−36J + 221400)F^{2}

+ (1916J + 24690528)F + (J^{2} + 50976J + 803894544)

Ψ_{5}(F, J) = F^{6} + (−J + 780)F^{5} + (−30J + 218940)F^{4}

+ (310J + 25968800)F^{3} + (13700J + 1177897200)F^{2}

+ (38424J + 22576632000)F + (J^{2} − 614000J + 155720872000)

Ψ_{7}(F, J) = F^{8} + (−J + 776)F^{7} + (−28J + 217756)F^{6}

+ (21J + 26195512)F^{5} + (6328J + 1276406726)F^{4} + (39361J + 31050881848)F^{3} + (−240492J + 404938789276)F^{2} + (−2176581J + 2721214073864)F + (J^{2} − 1711008J + 7427483226241)

Ψ_{11}(F, J) = F^{12} + (−J + 684)F^{11} + (55J + 157410)F^{10}

+ (−1188J + 12515580)F^{9} + (12716J + 75763215)F^{8} + (−69630J + 76077144)F^{7} + (177408J − 207606564)F^{6} + (−133056J − 34321320)F^{5} + (−132066J + 418524975)F^{4} + (187407J − 477130500)F^{3} + (−40095J + 270641250)F^{2} + (−24300J − 82012500)F + (J^{2} + 6750J + 11390625)

Ψ_{13}(F, J) = F^{14} + (−J + 772)F^{13} + (−26J + 216424)F^{12}

+ (−156J + 26333528)F^{11} + (1508J + 1359640022)F^{10} + (21658J + 39120460496)F^{9} + (39624J + 716780223796)F^{8} + (−612742J + 8956723925032)F^{7} + (−3355976J + 79070093432161)F^{6} + (454779J + 500196729175884)F^{5} + (43741490J + 2260671730897788)F^{4} + (95939974J + 7142292018579744)F^{3} + (−41335164J + 15009662255513328)F^{2} + (−291162600J + 18874201488396480)F + (J^{2} − 174668400J + 10755802087387200)

Ψ_{17}(F, J) = F^{18} + (−J + 690)F^{17} + (51J + 160191)F^{16}

+ (−1105J + 12849212)F^{15} + (13243J + 77940903)F^{14} + (−95659J − 24306702)F^{13} + (424065J − 489756655)F^{12} + (−1110355J + 856070496)F^{11} + (1454945J + 247945272)F^{10} + (−73746J − 4127455840)F^{9} + (−2450210J + 10326614640)F^{8} + (3131026J − 15993234432)F^{7} + (−1104830J + 18158824448)F^{6} + (−1073992J − 15889021440)F^{5} + (1392232J + 10788499200)F^{4} + (−557600J − 5622784000)F^{3} + (−2720J + 2154240000)F^{2} + (67200J − 537600000)F + (J^{2} − 16000J + 64000000)

Ψ_{19}(F, J) = F^{20} + (−J + 664)F^{19} + (76J + 143260)F^{18}

+ (−2622J + 9204360)F^{17} + (54454J − 176115066)F^{16} + (−761425J + 1108178952)F^{15} + (7598556J − 1742337316)F^{14} + (−55989713J − 13420942600)F^{13} + (310967414J + 79673435585)F^{12} + (−1317638334J − 133492721376)F^{11} + (4284347658J − 271425795648)F^{10} + (−10696404825J + 1738318231104)F^{9} + (20413753140J − 3257912161280)F^{8} + (−29485216120J + 528231178240)F^{7} + (31694225470J + 10718241992704)F^{6} + (−24698209440J − 26958821326848)F^{5} + (13397395520J + 36334713176064)F^{4} + (−4738229120J − 31060143636480)F^{3} + (973578240J + 16944463872000)F^{2} + (−91238400J − 5430382166016)F + (J^{2} + 1769472J + 782757789696)

Ψ_{23}(F, J) = F^{24} + (−J + 720)F^{23} + (23J + 179952)F^{22}

+ (−161J + 17282016)F^{21} + 441081120F^{20} + (3864J + 5678198784)F^{19} + (−5681J + 45492865088)F^{18} + (−46644J + 252605710080)F^{17} + (53084J + 1038071734272)F^{16} + (393024J + 3294356631552)F^{15} + (19136J + 8309302456320)F^{14} + (−1978368J + 16991995871232)F^{13} + (−2689666J + 28563290271744)F^{12} + (2882544J + 39839110889472)F^{11} + (11625488J + 46370418130944)F^{10} + (11002464J + 45154515419136)F^{9} + (−3833824J + 36762400456704)F^{8} + (−19783680J + 24919460020224)F^{7} + (−21906304J + 13946021740544)F^{6} + (−11787776J + 6353857806336)F^{5} + (−1554432J + 2304837156864)F^{4} + (2213888J + 642483486720)F^{3} + (1648640J + 129654325248)F^{2} + (516096J + 16911433728)F + (J^{2} + 65536J + 1073741824)

References

[1] A. O. L. Atkin. Several public email messages. unpublished, 1991–1992.

[2] Ian F. Blake, János A. Csirik, Michael Rubinstein, and Gadiel Seroussi. On the computation of modular polynomials. preprint, 1999.

[3] Ian F. Blake, Gadiel Seroussi, and Nigel Smart. Elliptic Curves in Cryptography, volume 265 of London Mathematical Society Lecture Note Series. Cambridge University Press, 1999.

[4] Henri Cohen. A Course in Computational Algebraic Number Theory, volume 138 of Graduate Texts in Mathematics. Springer, 1993.

[5] J.-M. Couveignes and François Morain. Schoof’s algorithm and isogeny cycles. In Leonard M. Adleman and Ming-Deh Huang, editors, Algorithmic number theory, Proceedings of the First International Symposium (ANTS-I) held at Cornell University, Ithaca, New York, May 6–9, 1994, number 877 in Lecture Notes in Computer Science, pages 43–58. Springer, Berlin, 1994.

[6] P. Deligne and M. Rapoport. Les schémas de modules de courbes elliptiques. In P. Deligne and W. Kuyk, editors, Modular Functions of One Variable II, volume 349 of Lecture Notes in Mathematics, pages 143–316. Springer, 1972.

[7] L. Dewaghe. Remarks on the Schoof-Elkies-Atkin algorithm. Mathematics of Computation, 67(223):1247–1252, 1998.

[8] Noam D. Elkies. Explicit isogenies. manuscript, Boston, MA, 1992.

[9] Noam D. Elkies. Elliptic and modular curves over finite fields and related computational issues. In D. A. Buell and J. T. Teitelbaum, editors, Computational Perspectives in Number Theory: Proceedings of a Conference in Honor of A. O. L. Atkin, (Chicago, IL, 1995), pages 21–76. AMS, 1998.

[10] Steven Galbraith and James McKee. The probability that the number of points on an elliptic curve over a finite field is prime. preprint CORR 99–51, University of Waterloo, Centre for Applied Cryptographic Research, 1999.

[11] Robin Hartshorne. Algebraic Geometry, volume 52 of Graduate Texts in Mathematics. springer, 1977.

[12] E. Howe. On the group orders of elliptic curves over finite fields. Compositio Math., 85:229– 247, 1993.

[13] T. Izu, J. Kogure, M. Noro, and K. Yokoyama. Efficient implementation of Schoof’s algorithm. In Kazuo Ohta and Dingyi Pei, editors, Proceedings of ASIACRYPT ’98, volume 1514 of Lecture Notes in Computer Science, pages 66–79. Springer, 1998.

[14] Toshitsune Miyake. Modular Forms. Springer, 1989.

[15] François Morain. Calcul du nombre de points sur une courbe elliptique dans un cops fini: aspects algorithmiques. Journal de Théorie des Nombres de Bordeaux, 7:255–282, 1995.

[16] Volker Müller. Ein Algorithmus zur Bestimmung der Punktzahl elliptischer Kurven über endlichen Körpern der Charakteristik größer drei. PhD thesis, Universität des Saarlandes, 1995.

[17] J.-M. Pollard. Monte Carlo methods for index computation (modp). Mathematics of Computation, 32(143):918–924, July 1978.

[18] René Schoof. Elliptic curves over finite fields and the computation of square roots mod p. Mathematics of Computation, 44:483–494, 1985.

[19] René Schoof. Counting points on elliptic curves over finite fields. Journal de Théorie des Nombres de Bordeaux, 7:219–254, 1995.

[20] Joseph H. Silverman. The Arithmetic of Elliptic Curves, volume 106 of Graduate Texts in Mathematics. Springer, 1986.

AT&T Shannon Lab, 180 Park Ave, Florham Park NJ 07932-0971

E-mail address: janos@research.att.com