Elliptic curves with a given number of points

Reinier Bröker and Peter Stevenhagen

Mathematisch Instituut, Universiteit Leiden, Postbus 9512, 2300 RA Leiden, The Netherlands

            reinier@math.leidenuniv.nl psh@math.leidenuniv.nl

Abstract. We present a non-archimedean method to construct, given an integer N ≥ 1, a finite field F_{q} and an elliptic curve E/F_{q} such that E(F_{q}) has order N.

1 Introduction

A classical theorem of Hasse from 1934 states that for an elliptic curve E defined over the finite field F_{q} of q elements, the order of the group E(F_{q}) of F_{q}-rational points is an integer in the Hasse interval

H_{q} = [q + 1 − 2\sqrt{q}, q + 1 + 2\sqrt{q}]

around q. If E is given in some standard way, say by a Weierstrass equation over F_{q}, there are several algorithms that compute the order of E(F_{q}). The 1985 algorithm by Schoof [8, 9] runs in time polynomial in log q, and in small characteristic p there are even faster p-adic algorithms due to Satoh [7] and Kedlaya [5].

The situation is rather different in the case of the following problem, which can be seen as an ‘inverse problem’ to the point counting problem.

Problem 1. Given an integer N ≥ 1, find a a finite field F_{q} and an elliptic curve E/F_{q} for which the number of F_{q}-rational points equals N.

As with other inverse problems, such as in Galois theory, this is mathematically a natural question to ask. In this particular case, an efficient solution to the problem would also be desirable in view of the need in current applications to construct elliptic curves having point groups satisfying various smoothness requirements with respect to their order. It is one of the reasons why we focus on the order N, and do not specify the finite field F_{q} as being part of the input. In addition, we will use the freedom with respect to the choice of a base field F_{q} to our advantage.

A necessary condition for our problem to be solvable for given N is clearly that N is contained in some Hasse interval H_{q}, so we would like the union ⋃_{q} H_{q} over all prime powers q to contain all positive integers. It is easy to see that the contribution to the union coming from the ‘true’ prime powers q that are not primes is negligible: it is contained in a zero density subset of Z_{≥1}. For this reason, we may and will restrict in the sequel to the case where the base field F_{q} is the prime field coming from a prime number q. In this particular case, all integers in H_{q} actually do occur as the group order of E(F_{q}) for some elliptic curve E, so N ∈ H_{q} is sufficient to guarantee the existence of a solution. (For arbitrary prime powers q there are often not enough supersingular curves to realize all orders congruent to 1 modulo the characteristic.)

For the equality Z_{≥1} = ⋃_{q prime} H_{q} we need to show that the primes are not too far apart, i.e., that the gap between consecutive primes q and q′ is roughly bounded by 4\sqrt{q} for large q. This is more than what is currently known to be true: even under assumption of the Riemann hypothesis the gap between consecutive primes can only shown to be of order O(\sqrt{q}(log q)^{2}). However, from a practical, algorithmic point of view there are always lots of primes q for which a large integer N is contained in H_{q}. Indeed, by the prime number theorem, we expect 1 out of every log N integers around N to be prime, so for large N the set of primes q having N ∈ H_{q} is on average of size 4\sqrt{N}/ log N, and finding such q is never a problem in practice.

Once we have found a prime q > 3 for which we have N ∈ H_{q} (we now require N > 1), there is the following naive algorithm to find an elliptic curve having exactly N rational points over F_{q}. Suppose that we are not in the easy cases where we have N = q + 1 (then any supersingular curve over F_{q} will do) or where one of the few curves with j-invariant 0 or 1728 has the right number of points. Then we try

E_{a} : y^{2} = x^{3} + ax − a (a ∈ F_{q}^{∗} \ {−27/4})

for random a as the Weierstrass equation of the desired curve until we find a curve having N points. More precisely, we write N = q +1−t and check whether for our a the point (1, 1) ∈ E_{a}(F_{q}) is annihilated by N = q + 1 − t or q + 1 + t. If it is, we check whether the number of F_{q}-rational points is indeed q + 1 ± t. For order N = q +1−t we are done, for order q +1+t not E_{a} itself but its quadratic twist has N points. Even though the distribution of the group orders of elliptic curves over F_{q} is not quite uniform, we expect to examine O(\sqrt{q}) = O(\sqrt{N}) curves E_{a} before we hit a curve having exactly N points. As the amount of time spent per a is usually very small, and certainly polynomial in log N, this yields a probabilistic algorithm with expected running time O(N^{\frac{1}{2} +o(1)}). It is quite practical for small values of N, but becomes unwieldy for N ≫ 10^{15}.

In the next section we briefly describe a classical deterministic algorithm based on complex multiplication methods which, although not asymptotically faster than the naive algorithm, can be improved in various ways. Our first improvement is a p-adic approach to complex multiplication based on the recent work of Couveignes and Henocq [3]. It is described in section 3, and illustrated by the explicit computation in section 4 of an ‘ANTS 6 curve’ having 2004061320040618 rational points. Our second improvement, in section 5, consists of using ‘small’ modular functions in this p-adic context to push the limits of what is feasible by p-adic methods. Although the resulting algorithm is still far from polynomial, its power is illustrated in section 6 by the computation of an elliptic curve having 10^{30} rational points.

2 Complex multiplication

A deterministic algorithm to produce an elliptic curve E over the prime field F_{q} having N ∈ H_{q} points is provided by the theory of complex multiplication. One writes N = q +1−t and observes that the Frobenius endomorphism F_{q} : E → E on the desired curve E/F_{q} satisfies the quadratic relation F_{q}^{2} − tF_{q} + q = 0 of discriminant Δ = t^{2} − 4q < 0 in End(E) = End_{\overline{F}_{q}}(E). Assume we are not in the supersingular case t = 0. Then our observation gives rise to an embedding

\begin{matrix}End(E) & ⟶ K = Q(\sqrt{Δ}) \\ F_{q} & ⟼ π_{q} = (t + \sqrt{Δ})/2 \end{matrix}

that maps F_{q} to a prime element π_{q} of trace t and norm q in the quadratic order O_{Δ} = Z[(Δ + \sqrt{Δ})/2] ⊂ K of discriminant Δ. By the Deuring lifting theorem, there exist a number field H ⊃ K and an elliptic curve Ẽ/H such that

1. there exists ϕ_{q} ∈ End(Ẽ) satisfying ϕ_{q}^{2} − tϕ_{q} + q = 0;

2. the prime q splits completely in H/Q, and for every prime 𝔮|q in H the reduced curve Ẽ mod 𝔮 is an elliptic curve having q + 1 − t points.

In fact, the reduction of the endomorphism ϕ_{q} ∈ End(Ẽ) above modulo a prime 𝔮|q in H yields the Frobenius endomorphism of the curve Ẽ mod 𝔮 over F_{q}.

The smallest field H ⊃ K over which a curve Ẽ satisfying 1 and 2 can be defined is the Hilbert class field of K. More explicitly, from a list of reduced binary quadratic forms [a, b, c] of discriminant D = disc(K), which can be viewed as a list of elements of the class group Cl(D) of discriminant D, we can form the class polynomial

F_{D} = ∏_{[a,b,c]∈Cl(D)} (X − j(\frac{−b + \sqrt{D}}{2a} )) ∈ Z[X]

of discriminant D. Here j : H → C denotes the well-known elliptic modular function on the complex upper half plane. Using sufficiently accurate complex approximations of the zeroes j(\frac{−b+\sqrt{D}}{2a}) of F_{D}, one may exactly determine F_{D} as it has integral coefficients. Once we have F_{D}, we are essentially done. Indeed, any zero of the irreducible polynomial F_{D} ∈ Z[X] generates the Hilbert class field H of K over K, and modulo q the polynomial F̄_{D} ∈ F_{q}[X] splits into linear factors. (In fact, this property of F_{D} is an excellent check for the correctness of any algorithm to compute F_{D}.) The zeroes of F̄_{D} in F_{q} are the j-invariants of the elliptic curves over F_{q} having endomorphism ring isomorphic to the ring of integers O_{D} of K. If \overline{j} ∈ F_{q} is one of these zeroes, we write down a curve with this j-invariant, and check whether it has q + 1 − t points. If it hasn’t, we have found a curve with q + 1 + t points, and (for j ≠ 0, 1728) its quadratic twist has N = q + 1 − t points.

This deterministic algorithm, although relatively simple, is not much faster than the naive algorithm, as for large D the class polynomial F_{D} is of degree h(D) = O(D^{\frac{1}{2} +o(1)}) and has coefficients of size O(D^{\frac{1}{2} +o(1)}), making the total running time of order O(D^{1+o(1)}). However, if we have the freedom to pick q on input N, it is often possible to find primes q for which the associated discriminant Δ = Δ(q) = (q −1−N)^{2} −4N is not only of size N^{\frac{1}{2} +o(1)} (this is simply done by picking q ∈ [N + 1 − 2\sqrt{N}, N + 1 + 2\sqrt{N}] close to the end points of the interval), but moreover has a large square factor. This leads to a field discriminant D of K which is quite a bit smaller than Δ, and makes the method feasible in cases where the naive method would take too long. As we currently cannot even prove the existence of a single prime in q ∈ [N + 1 − 2\sqrt{N}, N + 1 + 2\sqrt{N}], we certainly cannot prove the existence of q for which Δ(q) has large square factors and D is of order substantially smaller than O(N^{\frac{1}{2} +o(1)}).

3 A non-archimedean approach

The key feature of the complex multiplication method in the previous section is the computation of the class polynomial F_{D} ∈ Z[X] of the order O_{D} for suitable D. As the zeroes modulo q of F_{D} are j-invariants of curves E/F_{q} for which either E or its quadratic twist has exactly N points, this immediately solves our problem. In this section we achieve the computation of F_{D} in an other way, using p-adic instead of complex approximations of the zeroes of F_{D}. Working in a non-archimedean setting has the advantage that we no longer have to cope with the problem of rounding errors that arises in the complex case. It does require a p-adic substitute for the complex analytic method to evaluate j in CM-points of H using Fourier expansions, and this is provided by the recent work of Couveignes and Henocq [3] explained in this section.

Let N = q + 1 − t and Δ = t^{2} − 4q be as in the previous section, and D < −4 the discriminant of Q(\sqrt{Δ}). We first construct an elliptic curve E over a finite field F_{p} which has CM with O_{D}. As we want p to be as small as possible, we let s be the smallest positive integer of the same parity as D for which p = (s^{2} −D)/4 is prime. For D ≡ 1 mod 8 such p does not exist for parity reasons unless D = −7, and we pick the smallest positive s for which p = (s^{2} − 4D)/4 is prime instead. In practice we expect s to be small, at most a power of log |D|, such that p is of the same order of magnitude as D. Unfortunately, even under GRH proven upper bounds [3] for s are much weaker.

As there is a prime element π_{p} = (s + \sqrt{D})/2 (or π_{p} = (s + 2\sqrt{D})/2) of norm p and trace s > 0 in the order O_{D} = Z[π_{p}] (or O_{4D} = Z[π_{p}]), there exists an ordinary elliptic curve over F_{p} having CM by Z[π_{p}] and p + 1 ± s points over F_{p}. We can find such a curve E/F_{p} by applying the naive algorithm, as we saw in section 2 that |D| is much smaller than N. We have End(E) = O_{D} for D ≢ 1 mod 8. For D ≡ 1 mod 8 the ring End(E) is either equal to Z[π_{p}] = O_{4D} or to O_{D} ⊋ Z[π_{p}]. We are in the second case if all 2-torsion of E is F_{p}-rational, and in the first if it isn’t. As for D ≡ 1 mod 8 the class polynomials F_{D} and F_{4D} have the same degree and both generate the Hilbert class field H of K, they are both fine for our purposes. Alternatively, in case we have End(E) = O_{4D} there is a unique point P of order 2 in E(F_{p}), and we can replace E by the 2-isogenous curve E/⟨P ⟩ to achieve End(E) = O_{D}. We assume that we have End(E) = O_{D} in the sequel. Note that D is by definition fundamental.

By the Deuring lifting theorem, there exists a curve Ẽ/H and a prime 𝔭|p such that Ẽ reduces modulo 𝔭 to E and such that we have End(Ẽ) ≅ End(E). As p splits completely in H/Q the curve Ẽ is actually defined over Q_{p}. In fact, Ẽ is the unique elliptic curve over Q_{p} with reduction E and endomorphism ring End(Ẽ) ≅ End(E) ≅ O_{D}. It is this canonical lift Ẽ of E that we want to compute, as its j-invariant is a zero of our class polynomial.

Let Ell_{D}(F_{p}) be the set of \overline{F}_{p}-isomorphism classes of elliptic curves over F_{p} with endomorphism ring O_{D}. The j-invariants of the elements in Ell_{D}(F_{p}) are the zeroes of (F_{D} mod p), so Ell_{D}(F_{p}) is finite of order h(D) = #Cl(D). It can be identified with the similarly defined set Ell_{D}(Q_{p}) of \overline{Q}_{p}-isomorphism classes of elliptic curves over Q_{p} with endomorphism ring O_{D}. The j-invariants of the elements in Ell_{D}(Q_{p}) are the zeroes of F_{D}. Let C_{p} be the completion of an algebraic closure of Q_{p}, and write X_{D}(C_{p}) for the set of isomorphism classes of elliptic curves over C_{p} with the property that their reduction is an element of Ell_{D}(F_{p}). Then X_{D}(C_{p}) is a p-adic analytic space, as the j-invariant identifies it with a subset of C_{p}. It consists of h(D) discs of p-adic radius 1. Each disc contains exactly one element of Ell_{D}(Q_{p}), and this is the subset of C_{p} we want to compute. It consists of the (j-invariants of the) isomorphism classes of elliptic curves in X_{D}(C_{p}) having CM with O_{D}.

Let I ⊂ O_{D} be an O_{D}-ideal prime to p and E/F_{p} an elliptic curve in Ell_{D}(F_{p}). Then there is a separable isogeny E → E_{I} which has the subgroup E[I] of I-torsion points of E as its kernel. In this way, we obtain a bijection ρ̄_{I} : Ell_{D}(F_{p}) → Ell_{D}(F_{p}) that sends the isomorphism class of E to that of E_{I}. We obtain an action of the group I(p) of O_{K}-ideals prime to p on Ell_{D}(F_{p}), and since principal O_{K}-ideals act trivially, this action factors via the quotient map I(p) ↠ Cl(D). This makes Ell_{D}(F_{p}) into a principal homogeneous Cl(D)-space.

The fundamental idea in [3] is that the action of I(p) on Ell_{D}(F_{p}) admits a natural lift to an action of I(p) on X_{D}(C_{p}). More precisely, for I ∈ I(p) the map ρ_{I} : X_{D}(C_{p}) → X_{D}(C_{p}) is a p-adic analytic map that lifts ρ̄_{I} in the sense that on Ell_{D}(Q_{p}) ⊂ X_{D}(C_{p}), the restriction of ρ_{I} is the standard Galois action that factors via Cl(D). If I = (α) is principal with α ∈ O_{D} \ Z, then ρ_{α} = ρ_{I} stabilizes each disk around a CM-point in Ell_{D}(Q_{p}), and has the CM-point in the disk as its unique fixed point.

It is shown in [3] that the derivative of the map ρ_{α} in a point of Ell_{D}(Q_{p}) equals α/ᾱ, and this can be used to compute the j-invariant of a CM-point in Ell_{D}(Q_{p}) starting from an arbitrary point in the disk using a Newton iteration process. If E_{1} is any lift of E to C_{p}, we put

(1)  j(E_{k+1}) = j(E_{k}) − \frac{j(ρ_{α}(E_{k})) − j(E_{k})}{(α/ᾱ) − 1}  for k ∈ Z_{≥1}.

If α is small with respect to p, then (α/ᾱ) − 1 is a unit in Z_{p} and the sequence (1) converges to the j-invariant of the canonical lift Ẽ. In each step, the p-adic precision of the approximation is doubled.

Ell_{D}(F_{p}): X_{D}(C_{p}):

For the definition of ρ_{I} on the isomorphism class of E′ in X_{D}(C_{p}), we note that the subgroup E[I] of I-torsion points of the reduction E/F_{p} of E′ lifts canonically to a subgroup E′[I] of E′, and we put ρ_{I}(E′) = E′_{I}, with E → E′_{I} the isogeny with kernel E′[I]. This provides a lift of E[I] to a group scheme over the p-adic disk in X_{D}(C_{p}) lying over [E] ∈ Ell_{D}(F_{p}). More algorithmically, if E is given by a Weierstrass model, the subgroup E[I] can be described by a separable polynomial f̄_{I} ∈ F_{p}[X] having the x-coordinates of the affine points in E[I] as its zeroes. If I has norm n, then p ∤ n and f̄_{I} divides the n-th division polynomial of E in F_{p}[X]. Choosing a Weierstrass model for E′ reducing to that for E, we can lift f̄_{I} uniquely by Hensel’s lemma to a factor f_{I} of the n-th division polynomial of E′, and this factor describes the subgroup E′[I]. Note that ρ_{I} is the identity on X_{D}(C_{p}) if I is generated by an integer, and that, consequently, ρĪ is the inverse of ρ_{I}.

For the explicit computation of ρ_{α}(E′), we take an element α = a + bπ_{p} in O_{D} \ Z that is sufficiently smooth, i.e., a product of O_{D}-ideals L = (ℓ, a + bπ_{p}) of small prime norm ℓ ≠ p. Such an element is found by sieving in the set

{a + bπ_{p} : a, b ∈ Z, b ≠ 0, (a, b) = 1}.

Again, the smoothness properties are in practice much better than what can be rigorously proved [3]. For the computation of ρ_{L}(E′), we first compute the action of ρ̄_{L} on the reduction E = Ē′ ∈ Ell_{D}(F_{p}), i.e., the j-invariant of E_{L} = Ē′_{L}. The kernel E[L] of the isogeny E → E_{L} is a cyclic subgroup of order ℓ of E[ℓ] that is an eigenspace of the Frobenius morphism with eigenvalue −b/a ∈ F_{ℓ}. The corresponding polynomial f̄_{L} ∈ F_{p}[X] can be computed by the techniques used by Atkin and Elkies to improve Schoof’s original point counting algorithm, see [9]. These techniques also yield a Weierstrass model for E_{L}; we only need the j-invariant j(E_{L}) ∈ F_{p} of E_{L} = ρ̄_{L}(E). From the decomposition of ρ_{α} into ‘prime degree’ maps ρ_{L}, we obtain a cycle of isogenies

(2)  E ⟶^{ρ̄_{L_{1}}} E_{L_{1}} ⟶^{ρ̄_{L_{2}}} E_{L_{1}L_{2}} ⟶ … ⟶^{ρ̄_{L_{t}}} E_{(α)} = E.

To compute the action of L on the lift E′ of E, we do not lift f̄_{L} to some precision to a divisor of the ℓ-th division polynomial of E′, and then find a Weierstrass model and the j-invariant for E_{L} using Vélu’s formulas. As the division polynomial has degree (ℓ^{2} − 1)/2 for ℓ > 2, this would be rather timeconsuming. Instead, we exploit the ℓ-th modular polynomial Φ_{ℓ}(X, Y) ∈ Z[X, Y], which is of degree ℓ + 1 in each of the variables, and which we can precompute for a number of small primes ℓ. As there is an isogeny E′ → E′_{L} of degree ℓ, we know that j(E′_{L}) is a zero of Φ_{ℓ}(j(E′), Y) ∈ Z_{p}[Y]. Since we know the j-invariant of the reduction of E′_{L}, we also know which root to approximate in Z_{p}, and this reduces the lifting process to a simple Hensel lift of a zero of a polynomial of degree ℓ + 1.

For our Newton process in (1), we start from E ∈ Ell_{D}(F_{p}), and compute the cycle of F_{p}-isogenies in (2). We then lift E arbitrarily to a curve E_{1} over Q_{p}, the ‘1 digit precision’ approximation of the canonical lift. Now we compute lifts over Q_{p} of our F_{p}-isogenies in 2 digit precision, using the modular polynomials, and use the value of ρ_{α}(E_{1}) obtained to update E_{1} as in (1) to a 2 digit precision approximation E_{2} of the canonical lift. We continue this process of making (not really closed) cycles over Q_{p}, doubling the precision of the computation at each step, until we have the canonical lift with high enough accuracy (see [2] for an estimate of the required accuracy).

If we know the canonical lift Ẽ, we compute its Galois conjugates again via the modular polynomials. For this, we need small primes L that generate the class group. Under GRH, this can be done using primes L not exceeding the Bach bound 6 log^{2}(|D|), see [2]. In practice, this is never a problem. If we have all the conjugates to the required precision, we find by simple expansion of the product below the class polynomial

F_{D} = ∏_{[I]∈Cl(D)}(X − j(Ẽ_{I})) ∈ Z[X].

4 An elliptic curve for ANTS 6

We illustrate the working of our p-adic method by computing a tailor made elliptic curve for ANTS 6 having exactly N = 2004061320040618 points.

First we look for a small discriminant, so we write N = q + 1 − t for various primes q and search for a large square dividing Δ = t^{2} − 4q. In this example, the choice

N = 2004061230508291 + 1 + 89532326

yields Δ = −2^{3} · 3 · 619^{2} · 22567, so we take D = −2^{3} · 3 · 22567 = −541608 in this section. The corresponding class group Cl(D) has order 132.

Our goal is to compute the class polynomial F_{D}. We will do this p-adically, so we first find an elliptic curve over some F_{p} which has CM with O_{D}. The smallest integer s > 0 for which (s^{2} − D)/4 is prime is s = 2, so we have D = 2^{2} − 4p with p = 135403. We fix this value of p for the rest of this section.

We now apply the naive method to find a curve E_{a} : y^{2} = x^{3} + ax − a over F_{p} with trace of Frobenius 2. We find that E_{1737} has trace 2 and take this as our base curve E/F_{p}.

Next we determine which element α ∈ O_{D} \ Z we will use for the map ρ_{α}. The ideal (α) = (22539 + 4π_{p}) of norm 510353281 = 19^{2} · 29^{2} · 41^{2} factors as

(α) = L_{19}^{2} · L_{29}^{2} · L_{41}^{2} = (19, π_{p} + 6)^{2} · (29, π_{p} − 13)^{2} · (41, π_{p} − 13)^{2}.

We compute the action for the prime ideal L_{19}. The eigenvalue of the action of Frobenius on the 19-torsion is −6 ∈ F_{19}. If we evaluate the modular polynomial Φ_{19}(X, Y) in X = j(E) = 41556 ∈ F_{p}, we get a polynomial which has two roots over F_{p}, namely 19533 and 54827. From this we deduce that L_{19} sends j(E) to one of these two roots; we don’t know which one yet.

We just guess that the correct j-invariant is 54827 ∈ F_{p}. Following Elkies, we now compute the eigenspace S of the 19-torsion corresponding to this isogeny. We get the x-coordinates of the points on E in S as zeroes of

\begin{aligned}f̄_{L_{19}} = X^{9} & +29873X^{8} + 49874X^{7} + 131130X^{6} + 49222X^{5} + 46538X^{4}+ \\ & 111513X^{3} + 68602X^{2} + 126444X + 20947 ∈ F_{p}[X] \end{aligned}

Since we know that the eigenvalue for L_{19} is −6, we can now just check whether

(X^{p}, Y^{p}) = −6 · (X, Y)

holds for points in S, i.e., we compute both (X^{p}, Y^{p}) and −6 · (X, Y) in the ring

F_{p}[X, Y]/(f̄_{L_{19}}(X), Y^{2} − X^{3} + 1737X − 1737).

Note that the · means adding on the curve! In this example, it turns out that (X^{p}, Y^{p}) and −6·(X, Y) are not the same. It follows that the correct j-invariant of the L_{19}-isogenous curve is the other value 19533 ∈ F_{p}.

The action of L_{19} on the curve with j-invariant 19533 is now easier to compute: the modular polynomial has again two roots, but one of the roots has to be j(E). This root corresponds to the action of L̄_{19}, so we pick the other root. If we compute the entire cycle corresponding to L, we get:

41556 ⟶^{L_{19}} 19533 ⟶^{L_{19}} 100121 ⟶^{L_{29}} 86491 ⟶^{L_{29}} 40349 ⟶^{L_{41}} 32517 ⟶^{L_{41}} 41556.

We now lift E/F_{p} to E_{1}/Q_{p} by lifting the coefficients of the Weierstrass equation arbitrarily. The polynomial Φ_{19}(j(E_{1}), Y) ∈ Z_{p}[Y] has exactly two roots, one of which reduces to 19533 modulo p. We compute this root, which is the value of the L_{19}-action on j(E_{1}), to two p-adic digits of precision. We continue to lift the whole cycle to two p-adic digits of precision, and update j(E_{1}) according to formula (1) to obtain j(E_{2}). Starting from j(E_{2}), we now lift the cycle to four p-adic digits of precision, compute j(E_{3}) from this, and so on. We obtain

j(Ẽ) = 41556 + O(p)
= 41556 − 17953p + O(p^{2})
= 41556 − 17953p − 51143p^{2} − 17793p^{3} + O(p^{4})
= 41556 − 17953p − 51143p^{2} − 17793p^{3} + 45123p^{4} + 52596p^{5} + 18237p^{6}
+42211p^{7} + O(p^{8})
= 41556 − 17953p − 51143p^{2} − 17793p^{3} + 45123p^{4} + 52596p^{5} + 18237p^{6}
+42211p^{7} + 45716p^{8} + 58788p^{9} + 18836p^{10} − 4101p^{11} − 60004p^{12}
−24668p^{13} + 27527p^{14} − 58942p^{15} + O(p^{16}).

From the estimate in [2], we know that we need approximately 3000 decimals of accuracy. As we have 10^{3000} ≈ p^{585}, we compute j(Ẽ) up to 585 p-adic digits. The class group Cl(D) ≅ Z/2Z × Z/66Z is generated by a prime of norm 2 and a prime of norm 29, so we find all 132 conjugates of j(Ẽ) under Gal(H/K) up to 585 p-adic digits using the modular polynomials Φ_{2} and Φ_{29}. In the end, we expand the polynomial of degree 132 to find the class polynomial

F_{D} = ∏_{[I]∈Cl(D)}(X − j(Ẽ_{I})) ∈ Z[X]

We now compute a root of F_{D} ∈ F_{q}, and note that we can check our computations so far by testing whether F_{D} splits completely in F_{q}[X]. One of the 132 roots is j = 5215470850369 ∈ F_{q}, and an elliptic curve with this j-invariant is E_{a} : y^{2} = x^{3} + ax − a with a = 27j/(4(1728 − j)) = 1460967812073632 ∈ F_{q}. We know that E_{a} has CM with O_{D} and that its trace of Frobenius equals ±t. To test whether it actually equals t, we look at the order of P = (1, 1) ∈ E_{a}(F_{q}). We see that (q + 1 + t)P ≠ O = (q + 1 − t)P, so E_{a} must have trace equal to t. We conclude that the curve defined by

y^{2} = x^{3} + 1460967812073632x + 543093418434659

over F_{2004061230508291} has exactly 2004061320040618 rational points.

5 Using class invariants

A serious drawback of the complex multiplication method is that it requires the computation of the class polynomial F_{D}, which grows rapidly in size with D, and is already sizable for moderately small discriminants. Around 1900, Weber [11] computed generating polynomials for Hilbert class fields using the values at CM-points of modular functions of higher level instead of the j-function. These techniques have become important again in an algorithmic context, and many complications from Weber’s days are now well understood [4, 10]. The classical theory of class invariants is firmly rooted in complex analytic arguments, but much of it can be made to work in our non-archimedean setting. This section gives an indication of what can be done, leaving a fuller treatment to [1].

The complex multiplication method to compute the class polynomial F_{D} is based on the fact that its zeroes are the j-invariants of the elliptic curves in characteristic zero having complex multiplication with O_{D}. In the complex analytic setting, we simply list the complex lattices (up to scaling) giving rise to such curves and compute their j-invariants to sufficient accuracy using the q-expansion of j. In the p-adic setting, we first compute one such curve in the finite set Ell_{D}(F_{p}). Then we lift the action on Ell_{D}(F_{p}) of the group I(p) of O_{D}-ideals prime to p to an action on the set X(C_{p}) of their C_{p}-lifts in such a way that the induced action on the subset Ell_{D}(Q_{p}) of canonical lifts is the standard Galois action coming from Cl(D). This enables us to compute the finite subset Ell_{D}(Q_{p}) ⊂ X(C_{p}) ⟶^{j} C_{p} in C_{p} consisting of the zeroes of F_{D}. We use a Newton type lift of our curve in Ell_{D}(F_{p}) to Ell_{D}(Q_{p}) together with the explicit Galois action of Cl(D) on Ell_{D}(Q_{p}), which we handle by exploiting modular polynomials.

We now want to replace j by modular functions of higher level. These are elements of the modular function field F = ⋃_{n≥1} F_{n} as defined in [6, Chapter 6, Section 3]. Here F_{n} denotes the field of modular functions of level n over Q. It can be viewed as the function field of the modular curve X(n) over the cyclotomic field Q(ζ_{n}). Over F_{1} = Q(j), it is generated by the Fricke functions of level n, which are normalized x-coordinates of n-torsion points on an elliptic curve with j-invariant j.

The modular functions f we use are integral over Z[j], so they are given as the zero of some irreducible polynomial Ψ_{f}(X, j) ∈ Z[j, X]. If we specialize j to be the j-invariant of a curve Ẽ ∈ Ell_{D}(Q_{p}), then the roots of the polynomial Ψ_{f}(X, j(Ẽ)) ∈ H[X], which has integral coefficients, lie in the ray class field H_{n} of conductor n of K = Q(\sqrt{D}), with n the level of f. It is known that for many choices of ‘small’ f, one or more of these roots are class invariants that actually lie in the Hilbert class field H = H_{1} ⊂ H_{n} of K. If we can determine which roots end up in H, and compute the explicit Galois action of Cl(D) on such a root, we can compute its irreducible polynomial over K or Q just like we did this for j. In the complex analytic setting the tool to perform these tasks is Shimura’s reciprocity law [4, 10]. It tells us in which points τ ∈ K ⊂ H a function f should be evaluated to obtain a class invariant, and describes the conjugates of f(τ) over K as the values of conjugates of f over Q(j) in certain other points τ′ ∈ K. These values can be approximated in C using the q-expansions of f and its conjugates. Once we have computed the irreducible polynomial of f(τ) over K or Q, one can use the relation Ψ_{f}(f(τ), j(τ)) = 0 to obtain information on j(τ) itself.

In a p-adic setting, we cannot deal with j and f as functions on the complex upper half plane, and the expansion of modular functions as Fourier series in q = e^{2πiτ} has no non-archimedean analogue when dealing as we do with CMcurves, which have integral j-invariants. What we do have is an action from class field theory of the O_{D}-ideals coprime to n on the roots of Ψ_{f}(X, j(Ẽ)) ∈ H[X] in H_{n} associated to j(Ẽ). This action factors via the class group Cl(n^{2}D) of the order of discriminant n^{2}D, which is a non-maximal order for n > 1.

Example 1. The modular function j : H → C has a holomorphic cube root γ_{2} : H → C that is modular of level 3 and has Ψ_{γ_{2}}(X, j) = X^{3} − j. It is the unique root of X^{3} − j having a rational q-expansion. If D is not divisible by 3 and we write O_{D} = Z[τ] with τ + τ̄ = 0 mod 3, then γ_{2}(τ) is a class invariant, and its ‘size’ is only one third of that of j(τ).

If Ẽ : y^{2} = x^{3} + ax + b is in Ell_{D}(Q_{p}) and c_{1}, …, c_{4} are the 4 roots of its 3-division polynomial, then

(3)  \frac{−48a}{2a − 3(c_{1}c_{2} + c_{3}c_{4})}

is a cube root of j(Ẽ). As there are 3 ways to divide the roots c_{1}, …, c_{4} in two sets of two roots each, formula (3) yields 3 distinct cube roots of j. Note that there is no obvious way to single out an expression in (3) as being ‘the function’ corresponding to γ_{2}.

If I is an O_{D}-ideal prime to 3p, the isogeny Ẽ → Ẽ_{I} induces an isomorphism Ẽ[3] ⟶^{∼} Ẽ_{I}[3] on the 3-torsion subgroups, and maps each possible cube root of j(Ẽ) in (3) to some well-defined cube root of j(Ẽ_{I}). This is the Galois action of the Artin symbol of I, which maps j(Ẽ) to j(Ẽ_{I}), on these cube roots. It provides an extension of the map ρ_{I} : Ell_{D}(Q_{p}) → Ell_{D}(Q_{p}) to the set of the cube roots of (the j-invariants of) the elements in Ell_{D}(Q_{p}). As all these cube roots are p-adically integral, reduction of this map provides an extension of ρ̄_{I} : Ell_{D}(F_{p}) → Ell_{D}(F_{p}) to the set of cube roots in F̄_{p}, provided we have p ≠ 3. In a similar way, we have an extension of the map ρ_{I} : X_{D}(C_{p}) → X_{D}(C_{p}) to the set of cube roots.

Example 1 nicely illustrates that the cube roots of j are functions on the modular curve X(3), the points of which can be viewed as isomorphism classes of elliptic curves with complete 3-level structure: in order to have a well defined value not only the j-invariant of the curve but also some ‘ordering’ of 3-torsion points is required.

As the field F_{n} is generated over Q(j) by Fricke functions, a modular function f of level n is always a rational expression in j and the roots of the n-th division polynomial of a curve with j-invariant j. As in Example 1, the action ρ_{I} of an O_{D}-ideal I coprime to pn on X_{D}(C_{p}) naturally maps the roots of Ψ_{f}(X, j) for j ∈ X_{D}(C_{p}) to the roots of Ψ_{f}(X, ρ_{I}(j)). This observation suffices to treat modular functions providing class invariants by p-adic methods similar to that in section 3.

Suppose f is an integral modular function of level n that is known to provide class invariants for O_{D} at certain CM-points for O_{D} in H. Then we know that for every curve Ẽ ∈ Ell_{D}(Q_{p}), certain roots of Ψ_{f}(X, j(Ẽ)) ∈ Q_{p}[X] are class invariants, so they lie in Q_{p}. If E ∈ Ell_{D}(F_{p}) is the reduction of Ẽ, then Ψ_{f}(X, j(E)) ∈ F_{p}[X] will have roots in F_{p}, and we want to know which roots in F_{p} arise as the reduction of class invariants. Let β be a root, and assume for simplicity that Ψ_{f}(X, j(E)) is separable. (This is usually the case in practice as p is not too small, of size O(D^{1+ε}).) The roots of Ψ_{f}(X, j(Ẽ)) ∈ Q_{p}[X] all lie in H_{n}, and they are class invariants exactly when they are fixed by the Galois group

Gal(H_{n}/H) ≅ (O_{D}/nO_{D})^{∗}/O_{D}^{∗} = ker[Cl(n^{2}D) → Cl(D)].

It follows that β arises as the reduction of a class invariant if it is fixed by the maps ρ̄_{x} : Ell_{D}(F_{p}) → Ell_{D}(F_{p}) for a set of generators x of (O_{D}/nO_{D})^{∗}/O_{D}^{∗}. We can compute ρ̄_{x}(j(E)) as before when x is sufficiently smooth, and we need the extension of the action to the roots of Ψ_{f}(X, j(E)). In theory this can be done by working with explicit Weierstrass models and an explicit description of f in terms of n-torsion points as in Example 1. In practice we work with modular polynomials for prime degree ℓ isogenies relating the roots of Ψ_{f}(X, j_{1}) to that of Ψ_{f}(X, j_{2}) when j_{1} and j_{2} are j-invariants of ℓ-isogenous curves. Such modular polynomials are in many ways similar to the modular polynomials Φ_{ℓ}(X, Y) arising for the j-function, and they may be found using complex analytic methods involving q-expansions.

Suppose we find that a certain root β ∈ F_{p} of Ψ_{f}(X, j(E)) is the reduction of a class invariant β̃. Then we lift the pair (E, β) to the pair (Ẽ, β̃) consisting of the canonical lift Ẽ/Q_{p} and the class invariant β̃. In theory this can be done by lifting j(E) as in section 3, using cycles of smooth isogenies, and then compute the Hensel lift of β to the root β̃. This method makes use of the modular polynomials Φ_{ℓ} for j, which are big. In many cases the corresponding polynomials Φ_{ℓ}^{f} for f are much more pleasant to work with, so it is better to lift isogeny cycles not in terms of j-invariants but in terms of a root of Ψ_{f}(X, j), and find the resulting j-invariant from its corresponding root.

All that remains is computing the conjugates of the pairs (Ẽ, β̃) under Cl(D). This is done in exactly the same manner as before, using small primes that generate Cl(D). We are only interested in the conjugates of β̃, so we use the modular polynomials Φ_{ℓ}^{f} to compute these conjugates. In the end we expand the polynomial

F_{D}^{f} = ∏_{[I]∈Cl(D)}(X − β̃^{I}) ∈ Z[X] (or O_{D}[X].)

This polynomial splits again completely in F_{q}[X], and from a root in F_{q} we compute the corresponding j-value in F_{q}.

Example 2. We take for f the Weber function 𝔣, which is classically defined on H in terms of the Dedekind η-function as 𝔣(τ) = ζ_{48}^{−1}η(\frac{τ+1}{2})/η(τ). It is a modular function of level 48 of degree 72 over Q(j) with

Ψ_{𝔣}(X, j) = (X^{24} − 16)^{3} − jX^{24} ∈ Z[j, X].

For discriminants D ≡ 1 mod 8 with 3 ∤ D, it yields class invariants when evaluated at appropriate points τ ∈ H. In other cases small powers of 𝔣 often have the same property [4].

In principle we can express 𝔣 in terms of x-coordinates of 48-torsion points, but there is no need to do this. It suffices to find, for E/F_{p} a curve having CM with O_{D}, first a root β ∈ F_{p} of Ψ_{f}(X, j(E)) that is the reduction of a class invariant β̃, then a good approximation of the root β̃ of Ψ_{𝔣}(j(Ẽ), X) in Q_{p}, and finally good approximations in Q_{p} of the conjugates of β̃ under Cl(D). These questions ultimately reduce to computing the action of an ideal L of prime norm ℓ ∤ pn on pairs (E′, β′), with E′ in X_{D}(C_{p}) or Ell_{D}(F_{p}) and β′ a root of Ψ_{f}(X, j(E′) in Q_{p} or F_{p}. For E′ we know how to compute j(Ẽ′_{L}), for (β′)^{L} we use the fact that it is a zero of the modular polynomial Φ_{ℓ}^{𝔣}(β′, X) and of Ψ_{𝔣}(X, j(E′_{L})). Usually there are only two roots Φ_{ℓ}^{𝔣}(β′, X) that we need to consider, the correct one and the image under the action of L̄. It is of great help that the modular polynomials Φ_{ℓ}^{𝔣} are quite a bit smaller than the classical modular polynomials for j. For small ℓ they are really small, like

\begin{aligned}Φ_{5}^{𝔣}(X, Y) & = (X^{5} − Y)(X − Y^{5}) + 5XY \\ Φ_{7}^{𝔣}(X, Y) & = (X^{7} − Y)(X − Y^{7}) + 7(XY − X^{4}Y^{4}). \end{aligned}

For ℓ = 13 it takes at least two of these pages to write down Φ_{ℓ}, but we have

\begin{aligned}Φ_{13}^{𝔣}(X, Y) = & (X^{13} − Y)(X − Y^{13}) + 5 · 13XY \\ & +13(X^{2}Y^{12} + X^{12}Y^{2} + 4X^{10}Y^{4} + 4X^{10}Y^{4} + 6X^{6}Y^{8} + 6X^{8}Y^{6}). \end{aligned}

6 An elliptic curve having 10^{30} points.

We illustrate the power of the method presented in the previous section by constructing an elliptic curve having exactly N = 10^{30} rational points.

Just as in section 4, we first look for a suitable discriminant. We write N = q + 1 − t, and by looking at |t| slightly less than 2\sqrt{N} = 2 · 10^{15}, we find that for trace t = 1999999999167682 the number q = 10^{30} + t − 1 is prime and

Δ = t^{2} − 4q = −2^{12} · 3^{2} · 5^{2} · 17^{2} · 367^{2} · 3943 · 23537

has a large square factor leading to D = −92806391. The corresponding class group Cl(D) has order 15610. Computing the Hilbert class polynomial for our D would require an accuracy of 313618 decimals, which is clearly not practical.

Instead, we notice that we have D ≡ 1 mod 8 and D ≢ 0 mod 3, so we can use the classical Weber function 𝔣 to compute a class polynomial for H(K). We first compute an elliptic curve in Ell_{D}(F_{p}) for a ‘small’ prime p. We have D ≡ 1 mod 8, and the smallest s ∈ Z_{>0} with p = (s^{2} − 4D)/4 prime is s = 132, leading to p = 92810747. The first curve E_{a}/F_{p} of trace 132 we encounter is E_{1086} : y^{2} = x^{3} + 1086x − 1086 of j-invariant 37202456. As E has all three of its two-torsion points defined over F_{p}, its endomorphism ring is O_{D}, not O_{4D}.

We now have to determine which root β of the polynomial Ψ_{𝔣}(X, j(E)) = (X^{24} − 16)^{3} − j(E)X^{24} ∈ F_{p}[X] is the reduction of a class invariant β̃ ∈ Z_{p}. We are lucky since ±21677132 are the only two roots in F_{p}. Since −β̃ is also a class invariant, it does not matter which root we pick. We take β = 21677132.

For the smooth ideal inducing ρ_{α} we pick (α) = (−420+π_{p}), which factors as

(11, 8 + 2π_{p}) · (17, 4 + 2π_{p})^{2} · (23, 16 + 2π_{p})^{2} · (31, 13 + 2π_{p}) · (41, 30 + 2π_{p}).

Just as in section 4, we compute the cycle in F_{p} for the j-invariants:

37202456 ⟶^{L_{11}} 4967239 ⟶^{L_{17}} ⋯ ⟶^{L_{31}} 21402782 ⟶^{L_{41}} 37202456.

Using this cycle, we can also compute the cycle for β. For instance, the modular polynomial Φ_{11}^{𝔣}(β, Y) ∈ F_{p}[Y] has two roots: 32604444 and 60476019. In order to determine which root to take, we note that β is a root of Ψ_{𝔣}(X, j(E_{L_{11}})) = (X^{24} − 16)^{3} − j(E_{L_{11}})X^{24} ∈ F_{p}[X]. We find that 60476019 is the root we need. Continuing like this, we get the following cycle for β in F_{p}:

21677132 ⟶^{L_{11}} 60476019 ⟶^{L_{17}} ⋯ ⟶^{L_{31}} 53004472 ⟶^{L_{41}} 21677132.

Just as in section 4, we lift E/F_{p} to E_{1}/Q_{p} by lifting the coefficients of its Weierstrass equation. We could now compute the canonical lift Ẽ in two p-adic digits just like we did in section 4, and use that information to compute β̃ ∈ Z_{p} in two p-adic digits accuracy. Indeed, β̃ is a root of Ψ_{𝔣}(X, j(Ẽ)) ∈ Z_{p}[X], and since we have (β̃ mod p) = β, this root is just a Hensel lift of β. This approach has the disadvantage that we have to use modular polynomials for the j-function.

Instead, we lift β to a root β_{1} of Ψ_{𝔣}(X, j(E_{1})) ∈ Z_{p}[X]. Now we lift the cycle that we had for β ∈ F_{p} to a cycle for β_{1} ∈ Z_{p} by applying the small modular polynomials for 𝔣 once more. Since we know that β_{1}^{(α)} is a root of Ψ_{𝔣}(X, j(E_{1}^{(α)})), we can compute

j(E_{1}^{(α)}) = \frac{((β_{1}^{(α)})^{24} − 16)^{3}}{(β_{1}^{(α)})^{24}}

and use this value to update j(E_{1}) as in formula (1) to a value of the j-invariant of the canonical lift Ẽ that is accurate to two p-adic digits. Knowing j(Ẽ) mod p^{2}, we can lift β ∈ F_{p} to a root of Ψ_{𝔣}(X, j(Ẽ)) in Z_{p} that is accurate to two p-adic digits. We continue this process of doubling the precision until we have β̃ with sufficient accuracy. The first four cycles yield:

β̃ = 21677132 + O(p)
= 21677132 + 28966941p + O(p^{2})
= 21677132 + 28966941p + 7010373p^{2} + 31182954p^{3} + O(p^{4})
= 21677132 + 28966941p + 7010373p^{2} + 31182954p^{3} − 33808617p^{4}
+27519307p^{5} − 31601027p^{6} − 36195013p^{7} + O(p^{8})
= 21677132 + 28966941p + 7010373p^{2} + 31182954p^{3} − 33808617p^{4}
+27519307p^{5} − 31601027p^{6} − 36195013p^{7} − 8331811p^{8}
−33957007p^{9} − 18191700p^{10} + 5895954p^{11} − 42670221p^{12}
+23637278p^{13} − 40784695p^{14} + 7754196p^{15} + O(p^{16}).

We expect to need ⌈313618/72⌉ = 4356 decimals of accuracy, so we compute β̃ upto 550 p-adic digits. The class group Cl(D), which is cyclic of order 15610, is generated by a prime of norm 11. We can thus compute all the conjugates of β̃ under Gal(H/K) to 550 p-adic digits using the modular polynomial Φ_{11}^{𝔣}. In the end, we expand the polynomial of degree 15610 to find the class polynomial

F_{D}^{𝔣} = ∏_{[I]∈Cl(D)}(X − β̃^{I}) ∈ Z[X],

which we reduce modulo q to compute a root γ ∈ F_{q}. From γ, we compute the corresponding j-value and write down a curve Ē with that j-invariant. We then know that Ē has CM with O_{D}, but we still need to check whether the trace really equals t. This turns out not to be the case, so we conclude that the quadratic twist of Ē, given by

y^{2} = x^{3} + 669397215131271955483581235905x + 363369366443977510319399421188

over F_{q}, with q = 10^{30} + 1999999999167681, has exactly 10^{30} rational points. Checking that this curve indeed has the required number of points is an easy matter for the current point counting algorithms.

References

[1] Bröker, R.: Thesis, Universiteit Leiden, in preparation.

[2] Cohen, H.: A course in computational algebraic number theory, Springer Lecture Notes in Mathematics 138, (1993)

[3] Couveignes, J.-M. and Henocq, T.: Action of modular correspondences around CM points, Algorithmic Number Theory Symposium V. (2002)

[4] Gee, A.: Class invariants by Shimura’s reciprocity law, Journal de Théorie des Nombres de Bordeaux 11, (1999), 45–72

[5] Kedlaya, K.: Counting Points on Hyperelliptic Curves using Monsky-Washnitzer Cohomology, Journal Ramanujan Mathematical Society 16, (2002) 323–338

[6] Lang, S.: Elliptic functions, 2nd edition Springer Lecture Notes in Mathematics 112, (1987)

[7] Satoh, T.: The canonical lift of an ordinary elliptic curve over a finite field and its point counting, Journal Ramanujan Mathematical Society 15, (2000) 247-270

[8] Schoof, R.: Elliptic Curves over Finite Fields and the Computation of Square Roots mod p, Mathematics of Computation 44, (1985) 483–494

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

[10] Stevenhagen, P.: Hilbert’s 12th problem, complex multiplication and Shimura reciprocity Class field theory—its centenary and prospect (Tokyo, 1998) Adv. Stud. Pure Math., 30, Math. Soc. Japan, (2001) 161–176

[11] Weber, H.: Lehrbuch der Algebra, Vol. III. Chelsea reprint, original edition 1908