Antonin Leroux
email: antonin.leroux@polytechnique.org

Computation of Hilbert class polynomials and modular polynomials from supersingular elliptic curves

We present several new heuristic algorithms to compute class polynomials and modular polynomials modulo a prime p by revisiting the idea of working with supersingular elliptic curves. The best known algorithms to this date are based on ordinary curves, due to the supposed inefficiency of the supersingular case. While this was true a decade ago, the recent advances in the study of supersingular curves through the Deuring correspondence motivated by isogeny-based cryptography has provided all the tools to perform the necessary tasks efficiently. Our main ingredients are two new heuristic algorithms to compute the j-invariants of supersingular curves having an endomorphism ring contained in some set of isomorphism class of maximal orders. The first one is derived easily from the existing tools of isogeny-based cryptography, while the second introduces new ideas to perform that task efficiently for a big number of maximal orders at the same time. For each of the polynomials (Hilbert and modular), we obtain two algorithms. The first one, that we will qualify as direct, is based on the computation of a set of well-chosen supersingular j-invariants defined over F p 2 and uses the aforementioned algorithm to translate maximal orders to j-invariants as its main building block. The second one is a CRT algorithm that applies the direct algorithm on a set of small primes and reconstruct the result modulo p with the chinese remainder theorem. In both cases, the direct algorithm achieves the best known complexity for primes p that are relatively small compared to the discriminant (for the Hilbert case) and to the level (for the modular case). Our CRT algorithms matches the complexities of the state-of-the-art CRT approach based on ordinary curves, while improving some of the steps, thus opening the possibility to a better practical efficiency.

Introduction

Hilbert class polynomials and modular polynomials are central objects in number theory, and their computation have numerous applications. One field where these computations are of particular interest is cryptography. The main applications are to be found in elliptic curve cryptography and pairing-based cryptography, but we can also mention, more marginally, the recent field of isogeny-based cryptography.

Class polynomials, for instance, play a central role in the CM method, which is the main approach to find ordinary curves with a prescribed number of points over a given finite field (see [AM93,BS07]). This has applications to primality proving with the ECPP method and finding pairing friendly-curves with the Cocks-Pinch method.

Modular polynomials are related to isogenies between elliptic curves. Historically, they play a very important role in the SEA point counting algorithm [E + 98,Sch95] which remains one of the main algorithms used in elliptic-curve cryptography to generate cryptograhic curves. Moreover, the interest in isogenies have been renewed with the rise of isogeny-based cryptography. While most applications tend to use the more efficient Vélu formulas [Vél71], we can cite a few instances where modular polynomials have been considered. For example, it is used in the CRS key exchange [Cou06,RS06], the very first isogeny-based protocol, and we can also mention the OSIDH construction [CK19].

The goal of this work is to explore theoretical and practical improvements to the best-known algorithms to compute class polynomials and modular polynomials modulo prime numbers through the use of supersingular curves.

Related work. One of the main problems behind the computation of class polynomials and modular polynomials is the huge size of their coefficients over Z. There exists several algorithms of quasi-linear complexity [Eng09,CH02,Sut11,BLS12], but more often than not, memory is the real bottleneck in the concrete computations of those polynomials. In theory, size is less an issue when the result is needed modulo some prime p, but this is only true in practice if we have a way to skip entirely the computation over Z, which is not so easy to get.

Nonetheless, Sutherland [Sut11] proved that this could be done for class polynomials by a careful application of the CRT method. The result was later applied to modular polynomials by Bröker, Lauter and Sutherland (BLS) [BLS12]. The main advantage of the CRT method compared to other approaches is the low memory requirement (almost optimal in the size of the final output), and this is why this method has achieved the best practical results.

For both class and modular polynomials, the main tools used in the CRT algorithms from [Sut11,Sut12,BLS12] are /ordinary elliptic curves.

Supersingular curves had been considered in the context of the Hilbert polynomial computation in [BBEL08] and in the context of the modular polynomials computation in [CL05], but these methods were discarded over the time because they were slower than those based on ordinary curves.

The situation has changed with the recent interest on the connection between supersingular curves and quaternion algebras sparked by isogeny-based cryptography.

Since the work of Deuring [Deu41], it is known that endomorphism rings of supersingular curves in characteristic p are isomorphic to maximal orders in the quaternion algebra B p,∞ ramified at p and infinity, and that, conversely, every such maximal order types arises in this way. This is the first result of what is now called the Deuring correspondence. In this work, we are particularly interested in the task of computing the j-invariants of the (at most 2) supersingular elliptic curves over F p 2 having a given maximal order type as endomorphism ring.

The first concrete effort to realize that task is an algorithm of Cervino [Cer04] to compute the endomorphism rings of all supersingular curves in characteristic p. The complexity of this algorithm is O(p 2+ε ) and it becomes rapidly impractical. This algorithm was more recently improved by Chevyrev and Galbraith in [CG14] but the complexity is still O(p 1,5+ε ). As part of cryptanalytic efforts to understand the difficulty of various problems related to the Deuring correspondence, a heuristic algorithm of polynomial complexity in log(p) was introduced by Eisenträger, Hallgren, Lauter, Morrison and Petit [EHL + 18] . This algorithm builds upon the previous works of Kohel,Lauter,Petit and Tignol [KLPT14] and Galbraith,Petit and Silva [GPS17].

More concretely, these works prove that an isogeny can be efficiently computed between two supersingular curves of known endomorphism ring by translating the problem over the quaternions with the Deuring correspondence, solving the translated problems over quaternions, before translating back the solution as an isogeny. This can be applied directly to compute the j-invariants of all curves with an endomorphism ring contained in a given maximal order type by using one starting curve E 0 of known endomorphism ring (such a curve can always be computed efficiently with the CM method).

Recently, in [Rob22], Robert introduced new algorithms to compute modular polynomials from supersingular curves, this time using the high-dimensional isogeny techniques recently introduced in the context of isogeny-based cryptography cryptanalysis. Like us, Robert has a direct algorithm and a CRT algorithm. The direct algorithm has a quadratic complexity in the level ℓ, just as our direct algorithm, but the asymptotic complexity is somewhat differently balanced. The CRT algorithm reaches the same complexity has Sutherland's. Robert algorithm have the very nice advantage to have a proven complexity without any heuristic, but they require more memory, in particular the CRT algorithm.

Contributions. Our main contribution is to revisit previous methods to compute Hilbert and modular polynomials from supersingular curves with the recent algorithmic progress on the Deuring correspondence.

Our main sub-routine aims at translating a set of isomorphism class of maximal orders into their corresponding supersingular j-invariants under the Deuring correspondence. We introduce two algorithms, with different performance profiles, to perform that task.

With these new algorithms, we obtain an improvement over the asymptotic complexity of the class and modular polynomials computation in a wide range of primes below some upper-bounds that depend either on the discriminant of the class polynomial or the level of the modular polynomial. Moreover, we show that our new algorithm can also be used in the CRT method to reach the same complexity as ordinary curves, but with possibly better practical efficiency.

Technical overview

We start by looking at our main-subroutine that consists in the computation of the j-invariants of supersingular elliptic curves corresponding to some set of maximal order isomorphism classes (called maximal order types, see Definition 1). In the rest of this article, unless specified otherwise, a curve is considered to be a supersingular elliptic curve.

Maximal orders to j-invariants. We propose two algorithms dedicated to that task. Let consider that a set S of types is given in input, together with some prime p.

Our first algorithm is called OrdersTojInvariantSmallSet and it consists merely in a sequential execution of a sub-algorithm from [EHL + 18] that performs the desired translation for one type of maximal order. When everything is done carefully, it can be executed in O(log(p) 5+ε ) under experimentally verified heuristics detailed in [KLPT14] and related to the probability for numbers represented by some quadratic form to be prime. Thus, since OrdersTojInvariantSmallSet consists in #S executions of this sub-algorithm, the total heuristic complexity of Orders-TojInvariantSmallSet is O(#S log(p) 5+ε ). For a generic p and set of maximal order type S, we do not know how to do better than that. However, when S contains almost all maximal order types (the maximal size being upper-bounded by the number of supersingular curves), it becomes sub-optimal due to the amount of redundant computation performed along the way. In that case, it becomes much more practical to use an algorithm designed to sieve through the entire set of types, only focusing on the ones in S when they're met along the way. It requires a bit of care to perform this task in the most efficient manner but it can be done, and this leads to the algorithm OrdersTojInvariantBigSet of complexity O(#S log(p) 2+ε + p log(p) 1+ε ). This algorithm requires one heuristic that we detail in Section 3, as Claim 1. It is related to the expansion property of the supersingular isogeny graphs.

We stress that both algorithms are designed to work (and analyzed) for a generic prime p which is why they are so interesting.

The direct algorithms. Our heuristic algorithms OrdersTojInvariantSmallSet and OrdersTojInvariantBigSet can be used directly to compute the roots of class and modular polynomials modulo p (under the assumption that these roots are supersingular). The method is pretty straightforward: find the maximal order types corresponding to the desired roots, then, compute them with either Orders-TojInvariantSmallSet or OrdersTojInvariantBigSet. With the complexity we have stated, this is already enough to obtain an asymptotic improvement over existing generic methods when p is not too big (compared to the discriminant or level of the associated class or modular polynomial).

If we write S the "degree" of the polynomial (it is h(D) = O( |D| log(D) ε ) for Hilbert polynomials of discriminant |D| and O(ℓ 2 ) for modular polynomials of level ℓ), then we obtain the following complexity with OrdersTojInvariantSmall-Set:

O(S log 5+ε p + S log 2+ε S log p).

With OrdersTojInvariantBigSet, the complexity becomes O(S log 2+ε p + p log 1+ε p + S log 2+ε S log p).

In both cases, the latter term comes from the polynomial reconstruction step that must be performed to recover the polynomial from its roots. Note that the size of the output is O(S log p). In terms of space, the requirement is quasioptimal in both cases: so O(S 1+ε log p).

It is clear that the second algorithm will be better when p = O(S log(S)). However, whenever S = o(p) (which is often the case in applications), it will be better to use the variant with OrdersTojInvariantSmallSet.

The CRT for class polynomials. Let us take a prime p and a discriminant D < 0. We want an efficient algorithm to compute H D (X) mod p. Our main algorithm is essentially the same as the one introduced by Belding, Bröker, Enge and Lauter in [BBEL08], and later improved by Sutherland in [Sut11] .

Let us write O for the quadratic imaginary order of discriminant D < 0. We may assume that the factorization of D is known as computing it is negligible compared to the rest of the computation. We define P O to be a set of primes. We write B D for the bound on the bit-size of the coefficients of H D over Z.

Here is how the algorithm works:

1. Select some primes p 1 , . The main difference between the several variants of the CRT approach (including ours) lies is in the choice of the set P O and in the way to compute H D mod p i . In [BBEL08], different algorithms were proposed to handle the distinct cases of split and inert p i . In both cases, the H D mod p i are constructed from their roots. These roots are always j-invariants of elliptic curves in characteristic p i , but this is where the similarity ends. In the former case, the elliptic curves are ordinary and are defined over F pi , whereas in the latter case, we obtain supersingular elliptic curves defined over F p 2 i . The ordinary and supersingular case are very different and the resulting algorithms are also very different.

Sutherland [Sut11] improved the method from [BBEL08] for split primes by a careful choice of the primes and other tricks to improve the computation of H D mod p i .

We improve the method presented in [BBEL08] for non-split primes by making use of the Deuring correspondence with our OrdersTojInvariantBigSet algorithm to replace Cervino's algorithm [Cer04].

The CRT for modular polynomials. Let ℓ be a prime distinct from p. We want a CRT approach to compute Φ ℓ (X, Y ) mod p. It can be done in a very similar fashion to class polynomials by computing Φ ℓ mod p i for some p i in a set P ℓ and reconstruct the final polynomial via the CRT.

Bröker, Lauter and Sutherland (BLS in the rest of this article) in [BLS12] proposed to use primes of the form (t 2 -4v 2 Dℓ 2 ) with t, v, D ∈ N for which there is a very specific volcano structure involving ℓ-isogenies. This structure implies the existence of two distinct sets of ordinary curves defined over F pi : the curves with endomorphism ring isomorphic to O for some quadratic imaginary order O of discriminant D and class number bigger than ℓ + 2 and the curves with endomorphism ring isomorphic to Z + ℓO. Since the latter are ℓ-isogenous to the former, it is possible to recover the full Φ ℓ mod p i by computing the jinvariants corresponding to these two sets of curves. The volcano structure allows for efficient computation by minimizing the number of ℓ-isogeny computations.

For supersingular curves, the choice of primes is even easier than for class polynomials: we can use any primes p i that is big enough. As long as the number of supersingular curves is bigger than ℓ + 2 we will be able to recover the full modular polynomial. This idea has already been considered by Charles and Lauter in 2005 [CL05] but in a rather direct way (where each of the ℓ-isogeny involved is computed using the Vélu formulas).

We prove that using the Deuring correspondence and OrdersTojInvariantBigSet, we can avoid entirely any ℓ-isogeny computation and minimize the cost of elliptic curve operations.

Generic improvements to the CRT method. There are several ways to improve the CRT method in practical applications. First, alternative class polynomials and modular functions (with smaller height bounds) can be used instead of the standard Hilbert class polynomial and modular polynomials for the same practical purpose. This was the focus of the paper [ES10].

Second, for a number of applications such as the CM method and the SEA point counting algorithm, computing these polynomials is actually not necessary. What is really needed is the ability to evaluate them. Sutherland showed [Sut12,Sut13] that it was possible to do better than compute-then-evaluate in both cases, providing, in particular, an additional improvement in terms of memory requirement for those applications.

Using supersingular curves rather than ordinary ones should not prevent from applying all these practical improvements. For clarity's sake we focus on the simpler computation of the standard polynomials and leave to the reader the task of adapting these improvements to our new setting which should not be too daunting.

Organisation of the article. The rest of this paper is organized as follows: in Section 2, we introduce some background on isogenies, quaternion algebras and the Deuring correspondence. Then, in Section 3, we introduce our main new algorithm to compute efficiently j-invariants corresponding to maximal order types. In Section 4, we explain in details how this algorithm can be applied to the computation of class polynomial with the CRT method. In Section 5, we do the same for modular polynomials.

Acknowledgement. We thank Andrew Sutherland and anonymous reviewers for very useful feedback on this work.

Background material

Notations

Basic complexities. We write M Z (b) for the cost of multiplying two integers of less than b bits. For asymptotic complexities we consider M Z (b) = O(b 1+ε ). For instance, this covers the complexity of all arithmetic operations in a finite field F p of characteristic p of less than b bits.

Similarly, we write M P (b) for the cost (in terms of arithmetic operations over k) of multiplying two polynomials of degree smaller than b over a base field k. Depending on the size of b we will either use

M P (b) = O(b log(b) 1+ε or O(b 1+ε ).
Finally, the cost of fast interpolation algorithm for a polynomial of degree b is O(M ( P (b) log(b)).

Elliptic curves, quaternion algebras and the Deuring correspondence

More precise references on the topics covered in this section are: the book of Silverman [Sho94] for elliptic curves and isogenies, the book of John Voight [Voi18] on quaternion algebras and theoretical aspects of the Deuring correspondence, the thesis of Antonin Leroux [Ler22] for the algorithmic aspects of the Deuring correspondence.

Supersingular elliptic curves and isogenies. An isogeny ϕ : E 1 → E 2 is a nonconstant morphism sending the identity of E 1 to that of E 2 . The degree of an isogeny is its degree as a rational map (see [HS09] for more details). When the degree deg(ϕ) = d is coprime to p, the isogeny is necessarily separable and d = # ker ϕ. An isogeny is said to be cyclic when its kernel is a cyclic group. The Vélu formulas [Vél71] can be used to compute any cyclic isogeny from its kernel. For any ϕ : E 1 → E 2 , there exists a unique dual isogeny φ :

E 2 → E 1 , satisfying ϕ • φ = [deg(ϕ)].
Endomorphism ring. An isogeny from a curve E to itself is an endomorphism. The set End(E) of all endomorphisms of E forms a ring under addition and composition. For elliptic curves defined over a finite field F q , End(E) is isomorphic either to an order of a quadratic imaginary field or a maximal order in a quaternion algebra. In the first case, the curve is said to be ordinary and otherwise supersingular. We focus on the supersingular case in this article. Every supersingular elliptic curve defined over a field of characteristic p admits an isomorphic model over F p 2 . It implies that there only a finite number of isomorphism class of supersingular elliptic curves. The Frobenius over F p is the only inseparable isogeny between supersingular curves and it has degree p. We write π : E → E p . For any supersingular curve E, the property End(E) ∼ = End(E p ) is satisfied but we have E ∼ = E p if and only if E has an isomorphic model over F p .

Quaternion algebras. For a, b ∈ Q ⋆ we denote by H(a, b) = Q + iQ + jQ + kQ the quaternion algebra over Q with basis 1, i, j, k such that i 2 = a, j 2 = b and k = ij = -ji. Every quaternion algebra has a canonical involution that sends an element α = a 1 + a 2 i + a 3 j + a 4 k to its conjugate α = a 1 -a 2 i -a 3 j -a 4 k. We define the reduced trace and the reduced norm by tr(α) = α+α and n(α) = αα.

Orders and ideals. A fractional ideal I of a quaternion algebra B is a Z-lattice of rank four contained in B. We denote by n(I) the norm of I, defined as the Z-module generated by the reduced norms of the elements of I.

An order O is a subring of B that is also a fractional ideal. Elements of an order O have reduced norm and trace in Z. An order is called maximal when it is not contained in any other larger order. A suborder O of O is an order of rank 4 contained in O.

In this work, we will work with isomorphism classes of maximal orders in some quaternion algebra B and this is why we introduce the notion of type. The left order of a fractional ideal is defined as O L (I) = {α ∈ B p,∞ | αI ⊂ I} and similarly for the right order O R (I). A fractional ideal is integral if it is contained in its left order, or equivalently in its right order. An integral ideal is primitive if it is not the scalar multiple of another integral ideal. We refer to integral primitive ideals hereafter as ideals.

The product IJ of ideals I and J satisfying O R (I) = O L (J) is the ideal generated by the products of pairs in I × J. It follows that IJ is also an (integral) ideal and O L (IJ) = O L (I) and O R (IJ) = O R (J). The ideal norm is multiplicative with respect to ideal products. An ideal I is invertible if there exists another ideal The Deuring correspondence is an equivalence of categories between isogenies of supersingular elliptic curves and the left ideals over maximal order O of B p,∞ , the unique quaternion algebra ramified at p and ∞, inducing a bijection between conjugacy classes of supersingular j-invariants and maximal orders (up to equivalence) [Koh96]. Moreover, this bijection is explicitly constructed as E → End(E). Hence, given a supersingular curve E 0 with endomorphism ring O 0 , the pair (E 1 , ϕ), where E 1 is another supersingular elliptic curve and ϕ : E 0 → E 1 is an isogeny, is sent to a left integral O 0 -ideal. The right order of this ideal is isomorphic to End(E 1 ). One way of realizing this correspondence is obtained through the kernel ideals defined in [Wat69]. Given an integral left-O 0 -ideal I, we define the kernel of I as the subgroup

I -1 verifying II -1 = O L (I) = O R (I -1 ) and I -1 I = O R (I) = O L (I -1
E 0 [I] = {P ∈ E 0 (F p 2 ) : α(P ) = 0 for all α ∈ I}.
To I, we associate the isogeny

ϕ I : E 0 → E 0 /E 0 [I].
Conversely, given an isogeny ϕ, the corresponding kernel ideal is

I ϕ = {α ∈ O 0 : α(P ) = 0 for all P ∈ ker(ϕ)}.
In Table 1, we recall the main features of the Deuring correspondence.

Supersingular j-invariants over F p 2 Maximal orders in Bp,∞ j(E) (up to Galois conjugacy) O ∼ = End(E) (up to isomorpshim) (E1, ϕ) with ϕ : E → E1 Iϕ integral left O-ideal and right O1-ideal θ ∈ End(E0) Principal ideal Oθ deg(ϕ) n(Iϕ) φ Iϕ ϕ : E → E1, ψ : E → E1 Equivalent Ideals Iϕ ∼ I ψ Supersingular j-invariants over F p 2 Cl(O) τ • ρ : E → E1 → E2 Iτ•ρ = Iρ • Iτ
Table 1: The Deuring correspondence, a summary from [DFKL + 20].

Effective Deuring correspondence After establishing the nice theoretical results of the Deuring correspondence, it is natural to ask if we can obtain efficient algorithms to perform the translation between the two sides of our correspondence. This trend of work was started by Kohel,Lauter,Petit and Tignol in [KLPT14], and developed By Galbraith,Petit and Silva in [GPS17]. In [EHL + 18], Eisentrager, Haller, Lauter, Petit and Morrison provided the first complete picture of the situation (at least heuristically). It turns out that if we start from the quaternion side (either as a maximal order or an ideal), there are polynomialtime algorithms to compute the corresponding element (j-invariant, or isogeny).

In particular, Eisentrager et al. introduced a heuristic polynomial-time algorithm that computes the j-invariant corresponding to a maximal order type given in input. Henceforth, we call this algorithm SingleOrderTojInvariant(). It will be a crucial building block in one of our algorithm.

3 Computing j-invariants corresponding to maximal orders.

Let us fix some prime number p.

In this section, we introduce two algorithms to compute the j-invariants of supersingular curves over F p 2 corresponding to a set S of maximal order types in B p,∞ . By the Deuring correspondence, we know that each maximal order type in B p,∞ corresponds to one or two j-invariants of supersingular curve over F p 2 .

We will explain in Section 3.1 how to represent efficiently maximal order types as elements in some set H. Concretely, the input S to our algorithms will be given as some subset of H.

Our two algorithms presented in Sections 3.2 and 3.3 target two opposite situations with respect to the relative size of p and S. The first algorithm is called OrdersTojInvariantSmallSet. As the name suggests, it targets the case where #S/p is "small". It is a quite direct application of standard results on the effective Deuring correspondence, and handles each maximal order independently. The asymptotic complexity is O(#S log p 5+ε ) and works for any prime p and set S.

The second algorithm is called OrdersTojInvariantBigSet and it is more involved in both design and analysis. It targets the case where S is made of a significant portion of all O(p) possible types and is based on the idea that since S is big enough, the strategy that consists in going through the entire supersingular isogeny graph, collecting the j-invariants we want along the way, is quite optimal. Its complexity is O(S log p 2+ε + p log p 1+ε ). Hence, the cutoff between the two methods will be for some S with #S = Θ(p/ log p 2+ε ). Note that our second algorithm will be optimal when p/#S = Θ(log p 1+ε ).

Hashing to maximal order types

One of the important point for making our algorithms practical is to have a good way to handle sets of maximal order types and test if a type belong in some set of types.

For any maximal order O, we will represent Typ O by an invariant H(O). The purpose of this section is to introduce an efficiently computable invariant H(O) and the corresponding function H.

To derive an invariant for an isomorphism classes of lattices, it is quite natural to look at the smallest elements of that lattice. This idea was introduced by Chevyrev and Galbraith [CG14] in a related context, and their result was recently strengthened by Goren and Love [GL23]. Let us take a maximal order O. It can be shown (see [CG14] and [GL23, Theorem 1.4]) that if x 1 , x 2 , x 3 realize the successive minimas of the Gross lattice O T = {2x-tr(x)|x ∈ O}, then the values n 1 , n 2 , n 3 where n i = n(x i ) for i = 1, 2, 3 determine uniquely the maximal order O.

This is enough to obtain an invariant of size O(log p) as it can be shown that log(n i ) = O(log p) for all i ∈ {1, 2, 3}. If needed, for compactness, one can then apply some kind of hash function h : {0, 1} * → H where H is big enough to make the probability of a collision negligible over all maximal order types in B p,∞ .

For a generic statement, we take an arbitrary function h (which might be the identity) and assume its computational cost is negligible. Then, we define H as:

1. Compute O T . 2. Compute the three successive minimas x 1 , x 2 , x 3 of O T . 3. Output H(O) = h (n(x 1 ), n(x 2 ), n(x 3 )).
Proposition 1. The hash function H presented above can be computed in O(log(X) 1+ε ) when all the coefficients in the decomposition of the basis of O over 1, i, j, k are smaller than X.

Proof. This can be achieved by computing the successive minimas of O T with the algorithm to reduce ternary quadratic form of [ER01], or the L1 algorithm from [NSV11] to perform lattice reduction in small dimension.

3.2 Matching maximal orders and j-invariant: a generic algorithm.

Let us fix a prime p. We assume that a function H as introduced in Section 3.1 is defined and we assume that the underlying hash function h is such that there is no collisions over all maximal order types in B p,∞ .

The algorithm OrdersTojInvariantSmallSet is easily derived from the ideal-toisogeny algorithm described as [EHL + 18, Algorithm 12] and whose main building blocks where first introduced in [GPS17]. We provide a somewhat detailed description of this algorithm below. We refer to [GPS17,EPSV23] for a description of the main algorithmic building blocks.

1. Compute a curve E 0 and a maximal order

O 0 with End(E 0 ) ∼ = O 0 . 2. Set M = {}. 3. For each O ∈ S do: 4.
Compute an ideal I of powersmooth norm connecting O 0 and O 5.

Translate that ideal into the corresponding isogeny ϕ I . 6.

Compute the j-invariant j O of the codomain of ϕ I . 7.

M = M ∪ {H(O), j O }. 8. Return M .
Under some heuristic assumptions regarding the distribution of numbers represented by quadratic forms (see [KLPT14]), the complexity of this algorithm is O(#S log p 5+ε ).

For a concrete implementation, we propose to use the approach described in [EPSV23].

3.3 Maximal orders to j-inv: a faster algorithm when #S ≈ p.

OrdersTojInvariantSmallSet treats independently each order in S. This is fine when S is quite small compared to the number of maximal order types in B p,∞ , but it becomes less and less optimal as #S/p increases. Indeed, we end up redoing a lot of computations in that manner.

To overcome this problem when the ratio #S/p increases, we can try to mutualize as much computation as possible by using an approach that explores the entire isogeny graph and only targets the specific elements of S along the way. Since, our exploration of the isogeny graph is completely generic, we obtain an algorithm with a much better amortized cost per maximal order type.

More concretely, our idea is the following: take a smooth degree L such that all supersingular curves are L-isogenous to some starting curve E 0 of endomorphism ring in Typ O 0 for some maximal order O 0 . Compute all the O 0 -ideals of norm L and their right orders and select the ones contained in the set S. Then, enumerate efficiently through all the corresponding isogenies of degree L and collect the j-invariants of the codomains. Since quaternion operations cost less, we will do the exhaustive part over the quaternions, while minimizing the cost of elliptic curve operations by "selecting" the L-isogenies than we cannot avoid to compute.

Note that we do the ideal and isogeny phases in a simultaneous manner in OrdersTojInvariantBigSet. The good complexity we obtain will come from the care we take in computing all the required isogenies in the most efficient way possible, and in particular to avoid as many useless isogeny computations as possible. For that, the choice of L will be very important. In particular, if L = L 1 L 2 , by factorization of isogenies, we can compute all relevant L-isogenies by computing all L 1 -isogenies and only a subset of all the O(L 1 L 2 ) L 2 -isogenies. This subset obviously depends on S and the L 1 L 2 -isogeny computations account for the S log p 2 terms in the complexity (because we will choose L 2 to be smooth). The p log p covers the costs of the quaternion operations (which does not depend on S since we cover all maximal order types) and L 1 isogeny computations.

Before describing the algorithm in itself, we need to provide several properties and one heuristic claim that are going to be crucial for analyzing the algorithm complexity and proving its correctness.

Preliminary results and a heuristic assumption. Our first results target the degree L of the isogenies that we will use. As we explained above, this degree is crucial for optimizing the algorithm. To enable the fast computation of a lot of L-isogenies at the same time, we need it to be power-smooth (for efficient application of the Vélu formulas), but with coprime factors that are not too small.

Let us write Φ(N ) the number of cyclic isogenies of degree N for any N ∈ N.

Lemma 1. There exists a bound B and constants B 0 , C 0 , 1 < C 1 < C 2 such that for every number N > B, there exists a value n < B 0 log(N ) and integers (e i ) 1≤i≤n such that, if we define ℓ i to be the (n -i + 1)-th smallest prime and

L i = ℓ ei i then we have C 0 log(N ) ≤ Φ(L i ) < C 0 log(N ) for all 1 ≤ i ≤ n -1, (C 0 /2) log(N ) ≤ Φ(L n ) < C 0 log(N ) and C 1 N ≤ n i=1 Φ(L i ) < C 2 N .
Proof. We remind that when ℓ is a prime number and e ≥ 1 is an integer, Φ(ℓ e ) = ℓ e-1 (1 + ℓ).

Let us take C 0 , 1 < C 1 < C 2 three constants. We leave these constants unspecified for now. We start our reasoning without assuming anything on those constants, and we will encounter conditions on them during the proof. At the end, we will verify that these constraints can all be satisfied.

We also take a bound B that we assume to be "big enough" for various asymptotic inequalities to be verified. Let us finally take any integer N > B.

Let us construct two integers n, Λ n recursively from Λ 0 = 1 in the following manner: if

Λ i ≥ C 1 N , then set n = i, otherwise let Λ i+1 = Φ(λ i+1 )Λ i with λ i+1 = ℓ en-i n-i
where each ℓ i is the n -i + 1-th prime and e i is defined to be 1 if

Φ(ℓ i ) ≥ C 0 log(N ) or the biggest exponent such that Φ(λ i+1 ) < C 0 log(N ).
It is clear that this recursive algorithm always terminates and so we can get two integers n, Λ n in that manner for any N . Moreover, we have

C 0 /2 log(N ) ≤ Φ(λ 1 ) since ℓ 1 = 2 and C 0 log(N ) ≤ Φ(λ i ) for all other i (because if ℓ i ≤ Φ(λ i ) < C 0 log(N ), then ℓ n-i Φ(λ i ) = Φ(ℓ n-i λ i ) ≤ C 0 log(N )
and this is impossible by construction of the exponent e n-i ).

Then, there are two possibilities: either

Λ n < C 2 N or Λ n ≥ C 2 N . First case: Λ n < C 2 N . Let us set L i = λ n-i+1
for 1 ≤ i ≤ n, we will show that L 1 , . . . , L n is satisfies all the required properties.

In fact, the only thing that remains to be verified is that Φ(L i ) < C 0 log(N ) for 1 ≤ i ≤ n -1. By construction of the L i , it is clear that is suffices to prove that the n-th prime (the divisor of L 1 ) is smaller than C 0 log(N ) -1.

By taking logarithms of the inequality

C 0 log(N ) n ≤ Λ n < C 2 N , we otbain (n/2)(log(C 0 ) + log log(N )) < Λ n < log(C 2 ) + log(N ).
From there, it follows (since B is big enough) that n < 3 log(N )/ log log(N ).

From asymptotic results on the size of the n-th prime (see for instance [RS62, 3.13]), we get that we can assume ℓ 1 < 2n log(n) < 6 log(N ). Thus, if C 0 ≥ 7, the result is proven.

Second case: Λ n ≥ C 2 N . Since Λ n is too big, we will try to divide Λ n by a factor δ. Since N is big enough, we can assume n ≥ 4. We don't want to modify λ 1 since the bound is tight, but we can divide λ 2 , λ 3 by some powers of ℓ n-1 = 3, ℓ n-2 = 5. The idea is that we can remove up to O( log(N )) from λ 2 , λ 3 .

Since we have Λ n-1 < C 1 N , by the same reasoning as before we get n < 1 + 3 log(N )/ log log(N ). Hence, with ℓ 1 < 2n log(n) we get ℓ 1 < 7 log(N ).

Since, Φ(λ n ) is either smaller than C 0 log(N ) or equal to ℓ 1 + 1, we deduce that Φ(λ n ) < max(8, C 0 ) log(N ) and it suffices to choose C 0 ≥ 8 to ensure that Φ(λ n ) < C 0 log(N ).

From there, we deduce

Λ n = Λ n-1 Φ(λ n ) < C 0 C 1 N log(N ), from which we get C 1 N ≤ Λ n < C 0 C 1 N log(N ). (3.1)
Now, we want to find δ|Λ n such that C 1 N ≤ Λ n /δ < C 2 N . We start by using λ 2 . Let us take δ 2 as the biggest divisor of Φ(λ 2 ) such that

δ 2 < min(Λ n /(C 1 N ), C 0 log(N )/3).
If we assume C 2 > 3C 1 , then we can clearly take δ 2 ≥ 3 without having

Λ n /δ 2 ≤ C 1 N . If Λ n /δ 2 < C 2 N , then we are done with L n-1 = λ 2 /δ 2 and L i = λ n-i+1 for all other 1 ≤ i ≤ n, as we have C 1 N ≤ Λ n /δ 2 < C 2 N and the bound Φ(λ 2 )/δ 2 > C 0 log(N )
is satisfied by choice of δ 2 and the fact that Φ(λ 2 ) > C 0 log(N )/3.

If not and Λ n /δ 2 ≥ C 2 N , we will now try to remove a divisor of λ 3 . Before that, it is useful to try to upper-bound Λ n /δ 2 more precisely.

We have δ 2 > min(Λ n /(C 1 N ), C 0 log(N )/3)/3 (otherwise we could have multiplied δ 2 by 3). If the inequality Λ n /(C 1 N ) < C 0 log(N )/3 were to be true, then we would have δ 2 > Λ n /(3C 1 N ) and so Λ n /δ 2 < 3C 1 N . This is smaller than C 2 N by our assumption on C 1 and C 2 , and so this is a contradiction. Thus we must have Λ n /(C 1 N ) ≥ C 0 log(N )/3 and we can get the lower-bound δ 2 ≤ C 0 log(N )/9. If we plug that into Eq. (3.1), we obtain the following

C 1 N ≤ Λ/δ 2 < 9 C 0 C 1 N log(N ). (3.2)
Now, we define δ 3 as the biggest divisor of Φ(λ 3 ) such that

δ 3 < min(Λ n /(C 1 δ 2 N ), C 0 log(N )/ℓ 3 ).
Then, following the same reasoning as for δ 2 , we get that if we assume further 5C 1 < C 2 , then we can take a non-trivial δ 3 . Furthermore, we have that either Λ n /(δ 2 δ 3 ) < C 2 N and we are done (for the same reasons as before), or we must have

Λ n /(δ 2 C 1 N ) ≥ C 0 log(N )/ℓ 3 .
In that case, we can plug this in Eq. (3.2) to get

C 1 N ≤ Λ n δ 2 δ 3 < 225C 1 N. (3.3)
If we assume that C 2 > 225C 1 this is impossible and so we are done. In summary, suitable integers L 1 , . . . , L n can be found for any big enough N assuming that the constants C 0 , C 1 , C 2 satisfies:

-C 0 > 7. -C 2 > 3C 1 . -C 2 > 5C 1 . -C 2 > 225C 1 .
This system of inequalities can clearly be satisfied and so this proves the result.

Definition 2. Given any integer N > B, the integers L 1 , . . . , L n as defined in Lemma 1 is called a good-degree basis for N .

We derive the following result that will be useful in our analysis.

Proposition 2. Take an integer N > B with good degree-basis L 1 , L 2 , • • • , L n with n > 7. Then, there exists a constant C 3 such that

n-6 i=1 (n -i) i j=1 Φ(L j ) < C 3 N log(N ) 2 (3.4)
Proof. To prove the result, we are going to use the equalities

m i=0 X i = X m+1 -1 X-1 and m i=1 iX i = X( mX m+1 -(m+1)X m +1 (X-1) 2
) for any m and X. By our definition of a degree-basis in Lemma 1, we have that Φ(L j ) ≥ C 0 log(N ) for all 1 ≤ j ≤ n. Thus,

i j=1 Φ(L j ) < C 2 N/ n j=i+1 Φ(L j ) < C 2 N/( C 0 log(N )) n-i . Let us take X = 1/ C 0 log(N ). We have n-6 i=1 (n-i) i j=1 Φ(L i ) < C 2 N n-6 i=1 (n-i)X n-i . Then, we have n-6 i=1 (n -i)X n-i = n-1 i=6 iX i = X 6 n-7
i=0 (i + 6)X i . Then, we apply our two equalities to get

n-6 i=1 (n -i)X n-i =X 6 (6 X n-6 -1 X -1 ) + (n -7)X n-5 -(n -6)X n-6 + X (X -1) 2 ) =X 6 (n -1)X n-5 -(n -12)X n-6 -5X + 6 (X -1) 2 .
Since, asymptotically, we will have X < 1 while n will increase, we have that the leading term in the numerator of the fraction is 1, and so there exists

C ′ 3 such that n-6 i=1 (n -i)X n-i ≤ C ′ 3 X 4 . Thus, n-6 i=1 (n -i) i j=1 Φ(L i ) < C2C ′ 3 C 2 0 N log(N ) 2 . ⊓ ⊔
Our heuristic claim is about the size of L required to meet the condition that all supersingular curves are L-isogenous to some curve E 0 . We rewrite this under the Deuring correspondence as a condition on maximal orders and ideals.

Claim 1 There exists a constant C 4 such that for any prime p and maximal order O 0 in B p,∞ , given any number N > pC 4 , every maximal order type in B p,∞ is obtained as the type of the right order of a left integral O 0 -ideal of norm N .

Remark 1. Claim 1 is consistent with experiments regarding the diameter of the graph of supersingular 2-isogenies made in [ACNL + 21]. We also made some small experiments that seems to be consistent with that idea.

In case our claim fails, it is possible to use several starting curve E 0 to decrease the probability of missing some curve. If this is still not enough, and a few types are not obtained in this manner, it is always possible to apply SingleOrderTojInvariant to compute the remaining j-invariants without damaging too much the complexity.

For a given input p to OrdersTojInvariantBigSet, we assume the knowledge of a good degree-basis L 1 , • • • , L n (see Definition 2) of C 4 p (where C 4 is the constant in Claim 1). In OrdersTojInvariantBigSet, we will use the L i torsion points for all i. For supersingular curves, these points can always be defined over an extension F p m i . We have the following lemma to bound the value of all the m i from L i .

Lemma 2. Let p be a prime number and E 0 a supersingular curve over F p 2 . For any integer N , coprime with p, there exists Ê0 ∼ = E 0 over F p such that the torsion subgroup Ẽ0 [N ] is defined over an extension F p m of even degree m ≤ N over F p . Proposition 3. Assuming Claim 1, OrdersTojInvariantBigSet is correct.

Proof. By Claim 1 and our choice of L, we see that each maximal order types in B p,∞ is obtained as the right order of a L-ideal. Thus, we need to prove that our algorithm goes through every possible L-ideal and that it computes correctly the j-invariant associated with the right orders of those ideals. We will make our reasoning over isogenies and the Deuring correspondence will allow us to conclude the result over ideals.

Every cyclic L-isogeny can be factored as ϕ n • . . . • ϕ 1 where ϕ i is an isogeny of degree L i . When R i , S i is a basis of E[L i ], then all cyclic subgroups of order L i are generated by an element [C]R i + [D]S i . There is a 1-to-1 correspondence between cyclic subgroups of order L i and cyclic isogenies of degree L i . Since R 0,i , S 0,i is a basis of E 0 [L i ] and ϕ i-1 • • • • • ϕ 1 has degree coprime with L i , the two points R i , S i are a basis of E[L i ] and so this proves that our enumeration covers all possible isogenies of degree L i at each iteration of index i of the loop in line [START_REF]for all cyclic subgroups CRi + [D]Si of order Li in Ri[END_REF].

Thus, at the end of the loop we have covered all L-isogenies, and this means that if our ideal computation is correct, then we have covered all maximal order types and our algorithm is correct.

Remains to prove that the ideal I 1 ∩ . . . ∩ I i where each

I j = O 0 α 0 (C j + D j α), L j is the ideal corresponding to the isogeny ϕ i • . . . ϕ 1 where each ϕ j has kernel ϕ j-1 • • • • • ϕ 1 ([C j ]R 0,j + [D j ]S 0,j ) or equivalently that the kernel of ϕ i • . . . • ϕ i = i j=1 ([C j ]R 0,j + [D j ]S 0,j ).
For that, it suffices to prove the result for each coprime factor of the degree, so we need to prove that E 0 [I j ] = [C j ]R 0,j + [D j ]S 0,j . Let us go back to the definition of an ideal kernel given in Section 2. We have E 0 [I j ] = {P, β(I J ) = 0∀β ∈ I j {. Since I j contains L j O 0 , it is clear that the kernel must be a subgroup of E 0 [L j ]. Since multiplication in B p,∞ amounts to composition of the corresponding isogenies, it suffices to verify that ker α 0

(C j + D j α) ∩ E 0 [L j ] = [C j ]R 0,j + [D j ]S 0,j .
First, note that we have ker α 0 ∩ E 0 [J j ] = R 0,j by definition of α 0 and R 0,j . Then, we have (

[C j ] + [D j ]α)([C j ]R 0,j + [D j ]S 0,j ) = ([C j ] + [D j ](α)([C j ] + [D j ]α)(R 0,j ) = [n(C j +D j α)]R 0,j . This proves that we have [C j ]R 0,j +[D j ]S 0,j ⊂ ker α 0 (C j + D j α) ∩ E 0 [L j ].
By definition of α the scalar n(C j + D j α) is coprime with L j and so the endomorphism C j + D j α is a bijection on E 0 [L j ]. Thus, there cannot be another subgroup than [C j ]R 0,j + [D j ]S 0,j that is sent to R 0,j and this concludes the proof that ker α

0 (C j + D j α) ∩ E 0 [L j ] = [C j ]R 0,j + [D j ]S 0,j .
With that last fact, we have proven the point.

Complexity analysis. Below, we give as Theorem 1, a complexity statement for Algorithm 1. We derive this result from smaller statements for all the main Steps of Algorithm 1. The proofs of these statements include a more detailed description of the steps when needed. When a step has a complexity that will end-up being negligible before the total cost, we will sometimes not bother with a precise statement.

Note that all operations involving manipulations of the list M can be done efficiently in O(1) using adequate data structure so we do not analyze this part of the computation.

Since we look for an asymptotic statement, we assume that the size n of the good-degree basis used is bigger than 6 so we can apply Proposition 2.

Proposition 4. Under GRH, Step 2 can be executed O(poly(log p)).

Proof. The first step can be performed using an algorithm that was described as part of the proof of [EHL + 18, Proposition 3]. The idea is the following. Select the smallest fundamental discriminant d such that p is inert in the ring of integer

R d of Q( √ d).
Note that it can be proven that d = O(log p 2 ) under GRH. Then, we know that there are supersigular curves that admit an embedding of R d in their endomorphism ring. The j-invariants of these curves are the roots to the Hilbert class polynomial H d . It suffices to find one root of H d over F p to get a supersingular curve E 0 whose endomorphism ring will contain the Frobenius π and an endomorphism ι of norm d. This endomorphism can be found by computing all the isogenies of degree d with the Vélu formulae. Since the suborder 1, ι, π, ι • π has an index in O(d) = O(poly(log p)) inside End(E 0 ), and recovering the full endomorphism ring can be done in O(poly(log p)).

This proves the result for Step 2. The Steps 3 and 4 can be solved in a similar manner despite a final goal that is quite different. For Step 3, it is sufficient to find a β 0 such that that the matrix of the action of β 0 on a basis of the L-torsion has two distinct eigenvalues. Then, we can take α 0 = β 0 -λ where λ is one of the two eigenvalues. Note that this is equivalent to saying that β 0 needs to have two distinct eigenvalues mod L for all 1 ≤ i ≤ n.

For Step 4, a sufficient condition to obtain α such that gcd(n(x + yα), L) = 1 for all x, y with gcd(x, y, L) = 1 is to have that the matrix of α over a basis of the L i -torsion with no eigenvalues for all 1 ≤ i ≤ n.

The existence and value of these eigenvalues mod L i for a quaternion element β can be verified directly by computing the roots of the polynomial

X 2 + tr(β)X + n(β) mod L i .
Thus, to solve the two steps, we can apply the following method. First, for each L i dividing L, find one element β 0,i (resp. α i ) in O 0 /L i O 0 with two distinct (resp. no) eigenvalues mod L i . Since the ring O 0 /L i O 0 is isomorphic to M 2 (Z/L i Z) it is clear that a solution can be found by enumerating over the

L 4 i elements of O 0 /L i O 0 .
Then, the final element β 0 (resp. α) can be obtained by CRT (doing this coefficient-wise over the basis of O 0 ). The coefficients of β 0 (resp. α) will have size O(L) over the basis of O 0 and so we get the desired result by taking into account the coefficients of this basis over the canonical basis of B p,∞ . For each L i , we have to enumerate at most L 4 i quaternion elements, so the final complexity statements follows from the complexity of computing modular squareroot and the complexity of CRT. Proposition 6. The FOR loop in line 5 can be executed in O(poly(log(p))).

Proof. The loop is repeated O(log p) times. Computing a basis is a very standard task and it can be done in O(poly(log(p))) since m i = O(log p). The generator R 0,i can be computed in O(poly(log(p))) using the algorithm described in [GPS17] to find the kernel of an ideal. Evaluating the endomorphism can be done by evaluating a basis of End(E 0 ) and then performing the scalar multiplications corresponding to the coefficients of α in this basis. The first part can be done in O(poly(log(p))) with the ideas of the proof of Proposition 4 by choice of E 0 and the second can also be done in O(poly(log(p))) by the size bound given on the coefficients of α in Proposition 5.

The main computational task of OrdersTojInvariantBigSet is quite clearly performed during the loop in line 15. This FOR loop contains two inner loops. We will incrementally provide complexity statements for each of these loops in orders to clearly decompose the cost of each operations.

Proposition 7. At index i < n -6,
Step 26 produces a polynomial over F p 2 of degree smaller than L i in O(M P (L i ) log(L i )) operations over F p m i . This polynomial defines uniquely the isogeny ϕ. At index n -6 ≤ i ≤ n, L i = ℓ ei i for some prime ℓ i ≤ 13 (the 6-th prime) and Step 26 produces e i polynomials of degree ℓ i over F p 2 in O(e 2 i ) operations over F p m i . These polynomials uniquely defines the isogeny ϕ.

Proof. It can be shown with the Vélu formulaes [Vél71], that to represent an isogeny it suffices to compute its kernel polynomial, i.e., the polynomial whose roots are the x-coordinate of the points of the kernel and that this polynomial is always defined over F p 2 even if the kernel points are not.

From a single kernel point, the O(L i ) points of the kernel can be generated in O(L i ) operations over F p m i . Then, the kernel polynomial can be constructed with complexity O(M P (L i ) log(L i )) from its roots.

When n -6 ≤ i ≤ n, we can write L i = ℓ ei i for some prime ℓ i = O( 1) and then, we can factor our isogeny of degree L i as e i isogenies of degree ℓ i . All these isogenies can be computed in time O(e 2 i ) from a kernel generator (see [JDF11] for more on this topic).

Proposition 8. There exists a constant C such that, at any index i ≤ n -6, the number of F p -operations executed in Steps 27,28 of the FOR loop in line 19 is upper bounded by C(n -i) log pM P (log p). When i > n -6, we have the upperbound C(n -i)llog(p) 2 M P (log p) Proof. For each i < j ≤ n, we need to evaluate the polynomial produced by

Step 26 on the points R i+1 , . . . , R n and S i+1,...Rn . By Proposition 7, when i ≤ n -6, each evaluation costs O(L i ) operations over F p m j . The cost of arithmetic over F p m j in operations over F p is upper-bound by C ′ M P (m j ) for some constant C ′ . So we get the result from L i = O(log p) and m j = O(log p) by Lemma 2.

When i > n -6, there are e i = O(llog(p)) polynomials of degree ℓ i = O(1) and we get the desired result.

Remark 2. The Vélusqrt algorithm from [BFLS20] cannot be applied here because the kernel of the isogenies and the point on which the evaluation is performed do not live in the same extension.

Proposition 9. There exists a constant C such that, at any index i ≤ n -6, the number of binary operations executed in each execution of the FOR loop in line 19 is upper bounded by CΦ(L i )(n -i) log pllog(p)M P (log p)M Z (log p). When n -6 < i < n, the number of binary operations is smaller than

CΦ(L i )llog(p) 2 M P (log p)M Z (log p).
Proof. For each execution of this loop, the number of iteration is exactly Φ(L i ). Arithmetic over quaternion orders and ideals such as intersection and right order computation can be performed in O(M Z (log p)) (because the coefficients have size in O(log p) and these operations are simple linear algebra in dimension 4). Then, the hash can be computed in O(log p 1+ε ). Thus, the total cost of the loop is O(Φ(L i ) log p 1+ε ) for quaternion operations for any i. This is negligible compared to other operations when i < n.

The verification that H(O) ∈ S and insertion in the hash map operation can be done in O( 1) with the appropriate hash-map structure.

Then, there are the cost of operations over F p for the isogeny computation. To derive the total complexity, we apply Propositions 7 and 8. Arithmetic over F p takes O(M Z (log p)) binary operations. For any i, the kernel computation is negligible. When i ≤ n -6, the isogeny computation takes

O(Φ(L i ) log pllog(p)M P (log p)M Z (log p)),
then the evaluations take O((n -i)Φ(L i ) log pM P (log p)M Z (log p)).

Similarly, when n -6 < i < n, this cost is replaced by

O(Φ(L i )llog(p) 2 M P (log p)M Z (log p)).
Proposition 10. There exists a constant C such that, at any index i ≤ n -6, the number of binary operations executed in the FOR loop in line 17 is upper bounded by

C i j=1 Φ(L j )(n -i) log p 3+ε . When n -6 < i < n, it is upper bounded by C i j=1 Φ(L j )(n -i) log p 2+ε . When i = n, it is upper-bound by C(#S log p 2+ε + p log p 1+ε ).
Proof. There are Φ(L i ) cyclic subgroups of order L i . Thus at index i < n, the size of List is i-1 j=1 Φ(L j ) and the result follows directly from Proposition 9. When i = n, we perform the quaternion computations (intersection, right order and computation of the hash value) for all n j=1 Φ(L i ) subgroups. Thus, since we have n j=1 Φ(L j ) = O(p) by Lemma 1 and Claim 1, and the cost of quaternion operations is O(log p 1+ε ) as in the proof of Proposition 9, we get that the cost for quaternion operations is O(p log p 1+ε ). The isogeny computation is only performed when the right maximal order is contained in S, thus, we can upper-bound the number of times where an isogeny is computed by #S. As for all the i ≥ n -6, the cost of each L n -isogeny computation is O(llog(p) 2 M P (log p)M Z (log p)) and this proves the result.

Proposition 11. The loop in line 15 can be executed in O(#S log p 2+ε +p log p 1+ε ) binary operations.

Proof. The total cost of the loop is directly obtained by summing over all i the bounds in Proposition 10. We start by summing over all i ≤ n -6. We get an upper-bound on the number of binary operation by

C n-6 i=1 (n - i) i
j=1 Φ(L j ) log p 3+ε . So we get O(p log p 1+ε ) after applying Proposition 2. For n -6 < i < n, we get the upper-bound

C n-1 i=n-5 (n -i) i j=1 Φ(L j ) log p 2+ε .
Since by Lemma 1, Φ(L n ) > C 0 /2 log p, we have i j=1 Φ(L j ) ≤ C ′ p/ log p for some constant C ′ and any i < n. Thus, since there is a constant number of summands, we get the cost is in O(p log p 1+ε ) for those indices.

Finally, at i = n, we can apply directly the bound from Proposition 10. The final cost is O(#S log p 2+ε + p log p 1+ε ).

All the results above lead directly to the following theorem. Proof. The cost of Step 1 is negligible before all the other computations. The result follows from all the propositions proven in this section.

For the space requirement, the proof is quite similar. After each iteration of the loop for an index i ≤ n, the list contains i j=1 Φ(L i ) entries of the form E, K, P i+1 , . . . , P n , Q i+1 , . . . , Q n where E is an elliptic curve, K is an ideal of norm i j=1 L i and P j , Q j are points of order L i . Thus, this can be represented in O(

i j=1 Φ(L i )(log p + (n -i) log(p 2 )). Using Proposition 2 for i ≤ n -6 we see that this takes O(p log p) space. With Φ(L n ) = Θ(log(p)), we get that i j=1 Φ(L i ) = O(
p/ log p) and so the memory required for indices n -6 ≤ i < n is also in O(p log p). Finally, at index i = n the second term is zero and so we get O(p log p).

Finally, the space required by the list M is O(#S log p) = O(p log p), and this concludes the proof.

Good choice of primes. The analysis we provided above for both OrdersTo-jInvariantSmallSet and OrdersTojInvariantBigSet does not assume anything on the prime p. There are some "nice" choices of primes p for which we could basically gain a factor log p over all elliptic curve operations by having all the required torsion point defined over F p 2 (thus saving the cost of operations over big F pextensions). Since we are interested in a generic statement, we do not bother with these marginal gains.

In the context of applying OrdersTojInvariantBigSet to the CRT method, this idea will have its importance in the concrete choice of primes p i . However, due to the linear dependency in p, it does not appear possible to select all CRT primes among these "nice" primes. And so the CRT complexity will depend on the worst-case complexity of our algorithm OrdersTojInvariantBigSet.

Computation of the Hilbert class polynomial

In this section, we introduce our two algorithms to compute Hilbert class polynomial. In Section 4.1, we present the direct algorithm, and in Section 4.2, we introduce the CRT algorithm and we also provide a comparison with Sutherland's CRT algorithm from [Sut11].

Computing the class polynomial modulo a non-split prime.

Let us fix some negative discriminant D and a non-split prime p in Q( √ D) (and we assume further that p is coprime to the conductor of O). We write h for the class number of O. For simplicity, we assume below that the factorization of D is known so this step is not part of our estimates.

Our goal in this section is to explain how to compute the Hilbert class polynomial H D (X) mod p. This polynomial can be reconstructed from its roots that are j-invariants of some supersingular curves over F p 2 .

When p is split, the roots correspond to ordinary curves. In Sutherland's algorithm they are obtained in two main steps algorithm [Sut11]: start by identifying one root, and then enumerate through all the roots using the action of the class group Cl(O) through isogenies of small degree.

When p is non-split, we have three main steps. We start to do something very similar to what is done for ordinary curves, but instead of working directly with the roots, we first use the Deuring correspondence, and identify the roots from the the corresponding maximal orders of the quaternion algebra B p,∞ . The class group action of Cl(O) can be realized much more efficiently in that case because all operations in B p,∞ are essentially linear algebra.

The Deuring correspondence states that the set of maximal order we obtain in this manner are isomorphic to the endomorphism rings of the elliptic curves we want to compute. Thus, we constitue a set S D (p) of maximal order type and we can apply OrdersTojInvariantBigSet on this set to compute the j-invariants we need. This execution constitutes our third step.

The fourth and final step in our algorithm is common with the third step of the ordinary case: recover the polynomial H D from its roots. Note that this is done using standard polynomial arithmetic. This is described more precisely in Algorithm 2.

Algorithm 2 SupersingularHilbertComputation(D, p)

Input: A discriminant D < 0 and a prime p such that p is non-split in, and coprime to the conductor of,O the quadratic order of discriminant D. Output: HD(X) mod p. Proof. Correctness follows from [Onu21,Theorem 3.4]. We will now explain how each step can be performed and what is the best known complexity.

Step 1. This task has already been solved in the context of generating backdoor curves to the SIDH scheme [QKL + 21] and generating keys for the Séta encryption scheme [DFFdSG + 21]. First, we need to solve a quadratic equation over Q to find a, b, c, d ∈ Q such that Z[a+ ib + jc+ kd] is the quadratic order of discriminant D. This can be done using Simon's algorithm [Sim05] to solve quadratic forms in dimension 4. The complexity of Simon's algorithm is polynomial in the logarithm of the determinant once the factorization of the determinant is known and the size of the output is also logarithmic in the determinant. In our case, the quadratic form we consider is basically b, c, d, e → (qb 2 + p(c 2 + qd 2 ) -e 2 D and its determinant is equal to p 2 q 2 D2 f for some small integer f . The computation of the full factorization is sub-exponential in log(pD). Now that θ = a + ib + jc + kd has been computed, we need to find a maximal order O containing it. Let us take A as the smallest common denominator of a, b, c, d, we have A = O(poly(log(pD))). Then Aθ ∈ O 0 where O 0 is any maximal order containing the sub-order 1, i, j, k . Since Aθ ∈ O 0 , the right order of the ideal I = O 0 Aθ+O 0 C contains θ. We can set O = O R (I) and O can be computed in O(poly(log(p|D|))). Hence, this step can be performed in sub-exponential in log p|D| and is negligible compared to the rest of the computation.

Step 2. We go from one maximal order type to all maximal order types of interest by using the group action of the class group in a manner similar to what is used by Sutherland in [Sut11]. For that, we use a polycyclic representation of Cl(D) as introduced in [Sut11, Section 5]. But, in our case, instead of isogeny computation, we can simply use arithmetic over quaternions through the action of ideals of the form O(θ -λ) + Oℓ on the set of maximal orders containing θ which cover all maximal order types we need. Any group action computation for an ideal of norm ℓ takes O(log(ℓ)). Thus, using the same estimates than in [Sut11,Lemma 7] where we see that under GRH the biggest norm of an ideal involved a polycyclic representation is in O(log 2 (D)), we see that this part can be performed in

O(h log(|D|) ε ) = O( |D| log(|D|) ε ).
We can hash (with the function introduced in Section 3.1) all the maximal order types obtained in this manner to create the set

S D (p) in O( |D| log(|D|) ε - log(p) 1+ε ).
Step 3. This step consists simply in the execution of OrdersTojInvariantBigSet on the set S D (p) computed in Step 2. Thus, by Theorem 1 and the estimates on h, the complexity of this step is O( |D| log(|D|) ε log(p) 2+ε + p log(p) 1+ε ). Alternatively, it is possible to use the OrdersTojInvariantSmallSet algorithm and obtain a complexity of O( |D| log(|D|) ε log(p) 5+ε ). The reconstruction step. The complexity of this step is O( |D| log(|D|) 2+εlog(p) 1+ε ) as was proven in [Sut11]. The total complexity. Putting together all the results above, we get the complexity result. Space complexity. In terms of memory requirement, the polynomial reconstruction requires O(h(D) log p) space (since we only need to store two levels of the product tree at the same time). The set of maximal orders takes O(h(D) log p) space. For the variant with OrdersTojInvariantBigSet, the results follows from Theorem 1. When OrdersTojInvariantSmallSet is used, the computation of each j-invariant is done sequentially, and so the space requirement is optimal and is O(h(D) log p).

⊓ ⊔

Which version is optimal? We obtain two versions of SupersingularHilbertComputation depending on weither we use OrdersTojInvariantSmallSet or OrdersTojInvariantBigSet.

For a given value of D, the former has a better asymptotic complexity. However, for primes p in O( |D| log(|D|) 4+ε ), the latter has a better complexity. For a range of medium-sized primes, this algorithm will have the best known complexity. The cut-off with the CRT method (whose complexity is O(|D| 1+ε )) (as we will see below) happens for p = O(2 |D| 1/10 ).

4.2

Application to the CRT method and comparison with existing method.

In this section, we present and analyze the algorithm yielded by applying the CRT method from [Sut11] on top of the SupersingularHilbertComputation algorithm introduced in the previous section. The algorithm is pretty straightforward, and the only constraint on the selection of the small primes.

Algorithm 3 HilbertComputation(D, p)

Input: A discriminant D < 0 and a prime p. Output: HD(X) mod p.

1: Set BD as the upper-bound on the coefficients of HD(X) over Z.

2: Set P O = {}, P = 1, q = 3. 3: while P < BD do 4: if q is prime and q is non-split in and coprime to the conductor of O then 5:

P ← q • P , P ℓ ← P ℓ ∪ {q}.

6:

end if 7:

q ← q + 2 8: end while 9: Perform the precomputations for the explicit CRT mod p using P O . 10: for q ∈ P O do 11:

H O,q (X) ← SupersingularHilbertComputation(D, q). 12:

Update the CRT sums for each coefficient of H O,q (X). 13: end for 14: Perform the postcomputation for the explicit CRT to obtain HD(X) ∈ Fp[X]. 15: return Return HD(X).

Proposition 13. Under GRH and Claim 1, the complexity of HilbertComputation is O(|D| log 3+ε |D|). The space requirement is O( |D|(log 2+ε |D| + log p))

Proof. Under GRH, we can lower bound the probability that a given prime is non-split in O by 1/2. Then, following the same reasoning as in [Sut11], it can be shown that we have #P O = O(log B D /llog(B D )) and max q∈P O q = O(log(B D )). With the usual log B D = O( |D| log 1+ε |D|) that holds under GRH, we get that we can take O( |D|) primes with max q∈P O q = O( |D| log(|D|)). By Proposition 12, the dominant complexity of each CRT computation is the cost of the polynomial reconstruction. Thus, following the complexity estimates from [Sut11], we obtain that the overall complexity estimate is O(|D| log(|D|) 3+ε ). The space complexity follows from the space requirements computed in [Sut11] and Proposition 12.

Practical comparison with the CRT algorithm based on ordinary curves from [Sut11]. We can see that the complexity reported in Proposition 15 is the same as Sutherland's algorithm from [Sut11]. Even if the two algorithms have the same asymptotic complexity, they might not have the same practical efficiency. Below, we try to see which one could be faster in practice.

First, note that the set of CRT primes are very different between the two methods. Since the probability that a given prime can be selected in P O in our algorithm HilbertComputation is 1/2, the primes that we select are going to be much smaller than the ones used in Sutherland's algorithm. Now, we look at the concrete algorithmic steps. Let us start with the polynomial reconstruction step as it is the asymptotically dominant step. We remind that the concrete complexity of this step is O( |D| log(|D|) 2+ε log(p). It is pretty similar in both algorithms and we argue that the practical cost should be roughly the same. This is not completely obvious since the primes will not have the same size in the two cases and the roots are defined over F p for ordinary curves against F p 2 over supersingular curves. First, the size of the primes does not really matter because the product n i=1 p i have roughly the same size in the two cases and the complexity of the reconstruction is linear in log(p i ) for all i. Second, since in the supersingular case, the Galois conjugate (by the action of the Frobenius) of a root of H D mod p i is also a root, by building the remainder tree from polynomials of the form (X -j)(X -j pi ) ∈ F pi [X], we see that we can make the entire computation over F pi as in the ordinary case (and thus avoid the constant overhead brought by multiplications over F p 2 i ). We conclude from this brief reasoning that the reconstruction cost will be essentially the same in the two cases. Now, if we forget the reconstruction step, we see that using supersingular curve offers an asymptotic advantages. Indeed, in Steps 1, 2 and 3 of our algorithm, the dominant step is the execution of OrdersTojInvariantBigSet in Step 3, which has a O(|D| log(|D|) 2+ε ) complexity (if we consider the executions over all primes p i ∈ P D and we use log(p i ) = O(log(|D|))). In particular, this is smaller than the O(|D| log(|D|) 5/2+ǫ ) that dominates that part of the computation in Sutherland's algorithm (corresponding to the computation of one curve with the correct endomorphism ring). This is the first reason that suggests that the supersingular case might be more efficient than the ordinary one, but this is not the main one. The main reason behind the practical speed-up we hope to obtain is that we can use smaller primes. Indeed, the expected maximum of our primes is in O( |D| log(|D|)) (against O(|D| log(|D|) 1+ε for ordinary curves). Moreover, we can take all the small primes that satisfy the reduosity condition. In particular, we will be able to use a good portion of primes significantly smaller than |D|.

We hope that the very small primes will give a nice improvement in practice because for these primes, some of the roots will have big multiplicities, which should help perform every steps more efficiently in practice (for example there will be less than O( |D|) j-invariants to compute in that case) and it should improve the practical efficiency.

Note that there is also a good potential for practical improvement by carefully selecting the primes in P D and choosing the good-degree basis used for each of those primes in order to minimize the degree of the extension required to compute the isogenies in OrdersTojInvariantBigSet. A selection is also performed in the algorithm of Sutherland to help improve the cost of finding one curve with the good cardinal, so it would be natural to do the same thing in our case.

Even if OrdersTojInvariantBigSet proves to be too slow to beat the version of Sutherland by using only supersingular primes, it is clear that it is worth considering an hybrid set of primes P D containing a mix of supersingular and ordinary primes to obtain the best efficiency as the computation will be definitely very fast for a lot of small non-split primes.

Further improvement: batching class polynomial computation. OrdersTojInvariantBigSet can be easily modified to handle several sets S 1 , . . . , S k more efficiently than k executions of OrdersTojInvariantBigSet for each S i .

Thus, if we have several discriminants D 1 , . . . , D k , and a prime p in 1≤i≤k P D k . A good part of the computations performed to compute H D1 , . . . , H D k mod p can be done at the same time at a reduced cost.

Moreover, if some H Di have some common roots, the common divisors could be constructed once and for all.

Thus, our new method could be used to batch efficiently the computation of several class polynomial at the same time.

Computation of the modular polynomials

In this section, we introduce our new algorithms to compute modular polynomials. In Section 5.1, we present the direct algorithm, and in Section 5.2, we describe the CRT algorithm and compare the various existing algorithms of the literature.

Direct modular polynomial computation with supersingular curves

Let us take a prime level ℓ and a prime p. To make a direct computation of Φ ℓ (X, Y ) mod p it suffices to identify ℓ + 1 distinct j-invariants and for each of those, to interpolate Φ ℓ (j, Y ) from its roots that are the ℓ + 1 ell-isogenous j-invariants.

The idea of our algorithm for modular polynomials follows the same principle as the class polynomials algorithm. There is a slight difference because modular polynomials are bivariate but it does not change the generic principle of the algorithm. Indeed the full polynomial Φ ℓ mod p is interpolated from the Φ ℓ (j i , Y ) for ℓ + 1 j-invariants (j i ) 0≤i≤ℓ of elliptic curves defined in F p . Each univariate polynomial Φ ℓ (j i , Y ) is reconstructed from its roots, that are the j-invariants of curves ℓ-isogenous to j i .

In SupersingularModularComputation, we compute the maximal orders that correspond to all the j-invariants that we need under the Deuring correspondence, and then, we apply OrdersTojInvariantBigSet (or OrdersTojInvariantSmall-Set) to find the required j-invariants.

Algorithm 4 SupersingularModularComputation(p, ℓ)

Input: A prime p and a prime ℓ such that ⌈p/12⌉ + 1 > ℓ. Output: Φ ℓ (X, Y ) mod p. P (ji, Y ) ← j∈J i (Y -j). 12: end for 13: Reconstruct P (X, Y ) mod p from the P (j, Y ) for j ∈ J0. 14: return Return P (X, Y ). Proof. With the condition on the respective size of ℓ and p, we know there are enough supersingular j-invariants defined over F p 2 . We now briefly explain how each of the steps can be performed and what are their complexities.

By Claim 1, there is a value L = O(p) such that every supersingular jinvariant are L-isogenous to a given j-invariant j 0 . Thus, under the Deuring correspondence, the ℓ + 1 maximal order types can be computed in O(ℓ log(p) 1+ε ) by starting from a canonical maximal order type O 0 ⊂ B p,∞ (for which there exists formulas in [KLPT14] for any p). Thus, the complexity of Step 1 is O(ℓ log(p) 1+ε ).

The computation of an ℓ-ideal and of its right order can be done in O(log(ℓp)). Thus, the set of maximal order types S can be computed in O(ℓ 2 (log p + log ℓ)).

As show in [BLS12], the interpolation of all the Φ ℓ (j i , Y ), and the final interpolation of Φ ℓ (X, Y ) take O(ℓ 2 log ℓ 2+ε log p 1+ε ) .

The final complexity statement follows from the results shown previously on the complexities of OrdersTojInvariantBigSet and OrdersTojInvariantSmallSet.

The space complexity of the polynomial reconstruction is O(ℓ 2 log p) and the results follow from the space complexites of OrdersTojInvariantBigSet and OrdersTojInvariantSmallSet.

5.1.1 Which version is better ? Thus, we see that the first algorithm based on OrdersTojInvariantBigSet will be better for small values of p. We can estimate a cut-off for a value of p in O(ℓ 2+ε ). In that range of primes, our algorithm with OrdersTojInvariantBigSet has the best known generic complexity to compute Φ ℓ mod p. For primes bigger than that, it is better to use the variant with Orders-TojInvariantSmallSet to avoid the quasi-linear dependency in p. For a range of medium-sized primes, this algorithm will have best known complexity.

Space complexity. In terms of memory requirement our two algorithms are optimal and require O(ℓ 2 log p).

5.2 CRT methods to compute modular polynomials from supersingular curves.

In this section, we analyse the benefit of using supersingular curves in the CRT method and compare it with the algorithm described by Bröker,Lauter and Sutherland in [BLS12]. We refer the reader to the algorithm outlined in Section 1.1. Since the CRT is typically based on a lot of very small primes, we use the variant with OrdersTojInvariantBigSet.

Algorithm 5 ModularComputation(p, ℓ)

Input: A prime p, j ∈ Fp, a prime ℓ.

Output: Φ ℓ (X, Y ) mod p.

1: Let j be the integer in [0, p -1] equal to j mod p. 2: Set B ℓ = 2 6ℓ log ℓ+18ℓ+log(ℓ+2) . 3: Set P ℓ = {}, P = 1, q = 12ℓ. 4: while P < B ℓ do 5: if q is prime then 6:

P ← q • P , P ℓ ← P ℓ ∪ {q}.

7:

end if 8:

q ← q + 2 9: end while 10: Perform the precomputations for the explicit CRT mod p using P ℓ . 11: for q ∈ P ℓ do 12:

Pq(X, Y ) ← SupersingularModularComputation(q, ℓ).

13:

Update the CRT sums for each coefficient of Pq(X, Y ). 14: end for 15: Perform the postcomputation for the explicit CRT to obtain P (X, Y ) ∈ Fp[X, X]. 16: return Return P (X, Y ).

3. The CRT algorithm from BLS [BLS12]. 4. Our CRT algorithm ModularComputation. 5. The CRT algorithm from Robert [Rob22].

The two direct algorithms are the only known algorithms for generic primes whose complexity is only quadratic in ℓ. We have stated the complexity of the two variants of SupersingularModularComputation in Proposition 14. The complexity of Robert's direct algorithm is O(ℓ 2 log p log 2+3u ℓ) for some parameter u that might be 2 or 4 (it's not really clear which one).

On the other hand, the complexity of the three CRT algorithms is the same: it is cubic in ℓ and independand of p. BLS and ours are optimal in memory, while Robert's is not.

We will compare these various complexities by makin the value of p vary while ℓ remains fixed.

For small values of p, SupersingularModularComputation with OrdersTojInvariantBigSet will be the best algorithm as noted in the end of Section 5.1. However, due to the linear complexity of p, this algorithm will quickly be outperformed by the other algorithms as p grows. As soon as ℓ = o( √ p), it will be better to use SupersingularModularComputation with OrdersTojInvariantSmallSet. The log ℓ factor is only quadratic in the complexity of SupersingularModularComputation so it will remain better than Robert's direct algorithm when p remains not too big. However, since the log p factor has an exponent of 5 in SupersingularModularComputation while Robert's algorithm is linear in log p, the latter will eventually outperform our direct algorithm. The exact cross-point will depend on the value of u.

Finally, due to the fact that the CRT methods have a running time independant of p, they will end up being more efficient than the direct algorithms when ℓ = o(log p).

The practical efficiency of Robert's algorithm is hard to estimate, but it might not be competitive with the other two due to its much bigger memory requirement (which ends being quite problematic as modular polynomial are huge). We have already argued that we expect our method to outperform the method from BLS in practice.

Conclusion

We have introduced several new algorithms to compute modular polynomials of level ℓ and Hilbert polynomials of discriminant D modulo a generic prime number p from supersingular curves. The direct version of our algorithmsn have complexity in ℓ 2 and |D| respectively for generic primes. Depending on the relative size of ℓ,|D| and p, we exhibit improvements over the best known asymptotic complexities for a significant range of primes.

Moreover, when applied to the CRT method, we obtain an algorithm whose complexity is the same as previously known CRT method, but with the potential to give a practical improvement (in particular in the case of modular polynomials).

It remains to see how efficient our new algorithms are in practice. There are several practical challenges to overcome before providing an implementation of the proposed algorithms (in particular related to the field extensions involved in the computations of some isogenies), and this is why we leave the concrete implementation to future work.

Definition 1 .

 1 The type of an order O written Typ O is the isomorphism class of O.

  ). The conjugate of an ideal I is the set of conjugates of elements of I, which is an ideal satisfying II = n(I)O L (I) and II = n(I)O R (I). We define an equivalence on orders by conjugacy and on left O-ideals by right scalar multiplication. Two orders O 1 and O 2 are equivalent if there is an element β ∈ B ⋆ such that βO 1 = O 2 β. Two left O-ideals I and J are equivalent if there exists β ∈ B ⋆ , such that I = Jβ. If the latter holds, then it follows that O R (I) and O R (J) are equivalent since βO R (I) = O R (J)β. For a given O, this defines equivalences classes of left O-ideals, and we denote the set of such classes by Cl(O).

Proposition 5 .

 5 Step 3 and Step 4 can be executed in O(poly(log(p))) and the output α 0 and α have coefficients in O(poly(p)) over the canonical basis of B p,∞ .Proof. First, note that we can fix a basis of O 0 with coefficients (over the canonical basis of B p,∞ ) in O(p) (see [KLPT14, Section 2.3] for instance).

Theorem 1 .

 1 Under Claim 1, on input p and S, OrdersTojInvariantBigSet can be executed in O(#S log p 2+ε + p log p 1+ε ) binary operations. The space requirement of OrdersTojInvariantBigSet is in O(p log p).

1 :

 1 Find a maximal order O in Bp,∞ with O ֒→ O. 2: Use the action of Cl(D) to find SD(p), the set of types of maximal orders in Bp,∞ with an optimal embedding of O. To each type O ∈ SD(p) is associated a number nO ≥ 1 counting the multiplicity of the type O in this enumeration. The set of these multiplicites is write ND(p). 3: Compute JD(p) = OrdersTojInvariantSmallSet(p, SD(p)) (or use OrdersTojInvariantBigSet if this is more efficient). 4: Reconstruct H(X) = j D (p)∈J (X -j) mod p from SD(p), JD(p) and ND(p). 5: return Return H(X).Proposition 12. SupersingularHilbertComputation is correct and it can be executed in O( |D|(log(|D|) 2+ε log(p) 1+ε + log(|D|) ε log(p) 2+ε ) + p log(p) 1+ε ), binary operations and requires O((p+h(D)) log p) memory when OrdersTojInvariantBigSet is used in Step 1 (assuming GRH and Claim 1), and O( |D|(log(|D|) ε log(p) 5+ε + log(|D|) 2+ε log(p) 1+ε )), binary operations and requires O(h(D) log(p)) space with OrdersTojInvariantSmall-Set (under GRH and the heuristics from [KLPT14]).

Proposition 14 .

 14 SupersingularModularComputation can be executed in O(ℓ 2 (log ℓ 2+ε log p 1+ε + log p 2+ε ) + p log p 1+ε ) binary operations and requires O((ℓ 2 + p) log p) when OrdersTojInvariantBigSet is use in Step 8 (assuming GRH and Claim 1), and O(ℓ 2 (log p 5+ε + log ℓ 2+ε log p 1+ε )) binary operations and requires O(ℓ 2 log p) space with OrdersTojInvariantSmallSet (under GRH and the heuristics from [KLPT14]).

  1: Compute a set S0 = {O1, . . . , O ℓ+1 } of distinct maximal order types in Bp,∞. 2: Set S1, . . . , S ℓ+1 = {}, . . . , {}. S) if this is more efficient). 9: Divide J as J0, J1, . . . J ℓ+1 corresponding to S0, S1, . . . , S ℓ+1 . 10: for i = 1 to ℓ + 1 do 11:

	3: for i = 1 to ℓ + 1 do
	4:	Compute Ii,1, . . . , I i,ℓ+1 , the ℓ + 1 left Oi-ideals of norm ℓ.
	5:	Si ← {OR(I i,k )|1 ≤ k ≤ ℓ + 1}.
	6: end for
	7: S = 0≤i≤ℓ+1 S0. 8: Compute J	=	OrdersTojInvariantSmallSet(p, S)	(or	use
		OrdersTojInvariantBigSet(p,

Proposition 15. The complexity of ModularComputation is O(ℓ 3 log ℓ 3+ε ). The space requirement is O(ℓ 2 log p).

Proof. Following the same reasoning as in [Sut11], it can be shown that we have #P ℓ = O(B ℓ / log(B ℓ )) and max q∈P ℓ q = O(log(B ℓ )). The choice of B ℓ for the bound on the coefficients of Φ ℓ holds under GRH. We get that we can take O(ℓ) primes with max q∈P ℓ q = O(ℓ log(ℓ)).

By Proposition 14, for the complexity of SupersingularModularComputation with OrdersTojInvariantBigSet, and the complexity of the CRT steps described in [BLS12], we see that the complexity of each CRT computation takes O(ℓ 2 log(ℓ) 3+ε ). Thus, we obtain that the overall complexity estimate is O(ℓ 3 log(ℓ) 3+ε ). The space complexity follows in a similar manner.

On the cost of elliptic curve operations. The complexity of our algorithm is dominated by the polynomial reconstruction step. Note that this is also the case for the BLS algorithm.

However, unlike BLS, the rest of our computation is much more efficient. Indeed, if we look at the cost of the elliptic curve operation (the cost of all the elliptic curve operations inside the calls to SupersingularModularComputation), we see that the global complexity is only quadratic in ℓ. Indeed, we use primes q of size O(ℓ 1+ε ) for which there exists only O(ℓ 1+ε ) supersingular j-invariants. This means that we can simply compute them all with OrdersTojInvariantBigSet in quasi-linear time in ℓ, thus yielding a quadratic complexity over all CRT primes. This is why we expect our algorithm to produce a practical speed-up over the BLS method.

Batching the computation. Similarly to the Hilbert polynomial case, the set of small primes can be reused in the computations over various ℓ. In fact, the situation is even better in the modular case, because, apart from size, there are no restrictions on the primes. Thus, if we want to compute Φ ℓi for primes ℓ 1 , ℓ 2 , . . . , ℓ k of the same size, we will be able to use the same exact set P D for all the computations. Thus, only the polynomial reconstruction phase will be specific to each ℓ i , and the rest needs to be done only once.

Comparison between different existing methods

Here, we compare the complexity of various known methods to compute modular polynomials. Note that some of these methods are proven under some quite ad-hoc heuristics (ours in particular) while some others are proven rigorously. However, all the heuristics used in the various method are plausile and have been verified experimentally. Thus, we don't expect that these heuristics should be problematic for a practical comparison.

Here is the list of algorithms to consider:

1. Our direct algorithm SupersingularModularComputation.

2. The direct algorithm from Robert [Rob22].