Open Access Published by De Gruyter November 17, 2020

Algebraic approaches for solving isogeny problems of prime power degrees

Yasushi Takahashi , Momonari Kudo , Ryoya Fukasaku , Yasuhiko Ikematsu , Masaya Yasuda and Kazuhiro Yokoyama

From the journal Journal of Mathematical Cryptology

https://doi.org/10.1515/jmc-2020-0072

Abstract

Recently, supersingular isogeny cryptosystems have received attention as a candidate of post-quantum cryptography (PQC). Their security relies on the hardness of solving isogeny problems over supersingular elliptic curves. The meet-in-the-middle approach seems the most practical to solve isogeny problems with classical computers. In this paper, we propose two algebraic approaches for isogeny problems of prime power degrees. Our strategy is to reduce isogeny problems to a system of algebraic equations, and to solve it by Gröbner basis computation. The first one uses modular polynomials, and the second one uses kernel polynomials of isogenies. We report running times for solving isogeny problems of 3-power degrees on supersingular elliptic curves over 𝔽_p² with 503-bit prime p, extracted from the NIST PQC candidate SIKE. Our experiments show that our first approach is faster than the meet-in-the-middle approach for isogeny degrees up to 3¹⁰.

Keywords: Elliptic curves; Isogenies; Vélu’s formulae; Gröbner basis computation

MSC 2010: Primary: 14G50; Secondary: 94A60

1 Introduction

Since Koblitz [25] and Miller [27] independently proposed elliptic curve cryptography (ECC) in 1985, elliptic curves have been used in cryptography. Since 2000, pairings over elliptic curves have been an indispensable tool to construct cryptographic protocols and functional encryption schemes. Since 2006, isogenies between elliptic curves have began to be used in [9, 29, 34] for constructing several cryptosystems and hash functions. In particular, supersingular isogeny graphs were first proposed in [9] for security, which introduced the supersingular isogeny graph path-finding problem as a hard problem in cryptography. Actually, there exists a subexponential quantum algorithm [10] for computing an isogeny between ordinary elliptic curves since it is reduced to solving the abelian hidden shift problem. (Recently, the key exchange protocol based on group action [8] is regarded as a credible post-quantum system despite of the existence of a sub-exponential attack.) In contrast, the supersingular case is non-abelian and seems to be a promising candidate for PQC. (cf., The security of ECC and pairing-based cryptography relies on the discrete logarithm problem over elliptic curves, which could be solved in polynomial-time by Shor’s algorithm [31] using quantum computers.) In 2011, Jao and De Feo [22] introduced the key exchange protocol using using supersingular isogenies as a post-quantum candidate, based on pseudo-random walks in supersingular isogeny graphs. (See [11] for the connection between the hard problems in [22] and the path-finding problem in [9].) Other cryptographic functions were subsequently developed in [14, 23, 35]. In 2017, Jao et al. [21] submitted algorithms of supersingular isogeny key encapsulation, called SIKE, for the NIST PQC standardization process. It remains as a second-round candidate [28].

The template for the security of isogeny cryptosystems is the general isogeny problem [20]; Given two elements j,ȷ˜ of a finite field 𝔽_q, to find an isogeny ϕ:E→E˜ if exists, such that j(E) = j and j(E˜)=ȷ˜ where j(E) denotes the j-invariant of an elliptic curve E. In the supersingular case, it is sufficient to take 𝔽_p² as the base field for a prime p since every supersingular curve over F‾p is isomorphic to one defined over 𝔽_p² [33]. A variant of this problem is when the degree of ϕ is known, and it arises from the cryptanalysis of the hash function of [9], which requires computing isogenies of degree ℓ0e for some small ℓ₀ and large e. The most cryptographically-interesting variant is the supersingular isogeny Diffie-Hellman (SIDH) problem [14], assuring the security of [21, 22], in which supersingular elliptic curves over 𝔽_p² are chosen for a special prime p=ℓ1e1ℓ2e2f±1 and a lot of auxiliary information are given (see also [19]). The basic algorithm to solve the general isogeny problem for ordinary curves is due to Galbraith [18]. His procedure is (1) to reduce the problem to the case of elliptic curves E^′ and E˜′ whose endomorphism ring is maximal, and then (2) to construct an isogeny between E^′ and E˜′ by building isogeny trees. The time complexity of step (1) for E is negligible when the conductor [𝒪_K : End(E)] is smooth, where 𝒪_K denotes the maximal order of the imaginary quadratic field K=End(E)⊗ZQ with endomorphism ring End(E). For step (2), a sub-exponential quantum algorithm was discovered in [10]. In contrast, the meet-in-the-middle approach in [18] is applicable to the supersingular isogeny graph by building isogeny trees from E to E˜ directly. The best known quantum algorithm is due to Biasse et al. [3] for the supersingular problem, and its time complexity is exponential. For the SIDH problem, a faster quantum algorithm is given by Tani’s claw finding algorithm [36], but its time complexity is still exponential. (See also [24].) In recent, several related computational problems have been discussed in [15] for supersingular elliptic curves, their isogeny graphs, and their endomorphism rings.

Although time complexities of fast quantum algorithms has been discussed actively as described above, we consider practical approaches for solving isogeny problems with classical computers. Here we focus on the supersingular isogeny graph path-finding problem in [9] of prime power degree ℓ=ℓ0efor smallℓ0. The meet-in-the-middle approach is a practical way to solve the isogeny problem with time complexity equal to the square root of ℓ. It seems the most practical way at least in the supersingular case. (It is reported in [1] that van Oorschot and Wiener’s parallel golden collision search [38] can outperform the meet-in-the-attack for solving isogeny problems with parallel computation. See also [12] for improvements of the parallel golden collision search.) While the meet-in-the-middle approach is a generic way for solving graph problems, we propose two algebraic approaches for solving the isogeny problem in this paper. (cf., Several algebraic approaches are considered in [26] to attack the SIDH/SIKE protocol.) Our basic strategy is to reduce the isogeny problem to a system of algebraic equations. (It does not depend on the supersingularity.) Specifically, our first approach uses the modular polynomial of prime level ℓ₀. We divide a system of equations of modular polynomial into two parts like the meet-in-the-middle approach, and compute their Gröbner bases to efficiently find j-invariants of intermediate curves between given two isogenous elliptic curves EandE˜. In contrast, our second approach takes an intermediate curve E ₀ between EandE˜, and consider kernel polynomials F (x) and F˜(x˜) of two isogenies φ:E→E0andφ˜:E˜→E0. Since the curve E₀ is unknown, we regard its Weierstrass coefficients as variables, and also represent kernel polynomials F(x)andF˜(x˜) as certain multivariate polynomials, based on Schoof’s work [30]. Furthermore, by using Vélu’s formulae [6, 39], we represent isogenies φ and φ˜ as rational functions over multivariate polynomial rings, and set up algebraic equations to determine the Weierstrass coefficients of E₀. Moreover, we report running times of our algebraic approaches for solving the isogeny problem of 3-power degrees on supersingular elliptic curves over 𝔽_p² with 503-bit prime p, extracted from SIKE-p503 parameters [21], to compare with the meet-in-the-middle approach in performance.

2 Mathematical background

In this section, we review some basic definitions and properties of elliptic curves and isogenies, which shall be required in Section 3 below.

2.1 Elliptic curves over finite fields

An elliptic curve E over a finite field 𝔽_q of characteristic p ≥ 5 is defined by the (short) Weierstrass form y² = x³ + ax + b with a, b ∈ 𝔽_q and discriminant Δ(E) = −16(4a³ + 27b²) ≠ 0. Its j-invariant is defined as j(E)=−1728(4a)3Δ(E). Two elliptic curves are isomorphic over F‾q if and only if they both have the same j-invariant, where F‾q denotes the algebraic closure of 𝔽_q. Moreover, given an element j₀ of 𝔽_q, there exists an elliptic curve over 𝔽_q with j-invariant equal to j₀. The set of 𝔽_q-rational points

E(Fq)=(x,y)∈Fq2:y2=x3+ax+b∪{OE}

forms an abelian group, where 𝒪_E denotes the infinity point on E. (See [33, Chap. III] for the group law.) The number of 𝔽_q-rational points on E, denoted by #E(𝔽_q), is finite, and it is represented as #E(𝔽_q) = q + 1 − t where t denotes the trace of the q^th-power Frobenius map (see [33, Chap. V] for the map). It satisfies |t|≤2q by Hasse’s theorem. An elliptic curve E over 𝔽_q is said supersingular if the characteristic p divides the trace t. Otherwise E is said ordinary. Every supersingular elliptic curve over F‾p has its j-invariant defined over 𝔽_p² [33, Thm. 3.1, Chap. V], and hence it is isomorphic (over F‾q ) to one defined over 𝔽_p². Furthermore, there are about p12 isomorphism classes of supersingular elliptic curves over F‾p.

2.2 Isogenies and Vélu’s formulae

Let EandE˜ be two elliptic curves over a finite field 𝔽_q. A morphism ϕ:E→E˜ satisfying ϕ(OE)=OE˜ is called an isogeny. Two elliptic curves EandE˜ are called isogenous if there is a non-zero isogeny between them. Tate’s theorem [37] states that EandE˜ are isogenous over 𝔽_q if and only if #E(Fqk)=#E˜(Fqk) for any positive integer k. Every non-zero isogeny ϕ:E→E˜ induces an injection of function fields ϕ∗:F‾q(E˜)→F‾q(E) (see [33, Chap. III]). The degree of a non-zero isogeny is defined as deg ϕ=[F‾q(E):ϕ∗F‾q(E˜)]. We say that ϕ is separable if the finite extension F‾q(E)/ϕ∗F‾q(E˜) is separable. A non-zero isogeny is separable if its degree is not divisible by the characteristic of the base field 𝔽_q.

A non-zero isogeny ϕ:E→E˜ induces a (surjective) group homomorphism from E(F‾q) to E˜(F‾q) Thm. 4.8, Chap. III], and its kernel is a finite subgroup of E(F‾q), denoted by E[ϕ]. It satisfies deg ϕ = #E[ϕ] if ϕ is separable. Conversely, given any finite subgroup SofE(F‾q), there are a unique elliptic curve E˜ and a separable isogeny ϕ:E→E˜withE[ϕ]=S [33, Prop. 4.12, Chap. III]. The curve E˜ is denoted by the quotient E/S. Vélu [39] showed how to explicitly represent the isogeny ϕ:E→E˜=E/S and the Weierstrass equation for E˜. A non-zero isogeny ϕ:E→E˜ is normalized if ϕ∗(ωE˜)=ωE, where ωEandωE˜ denote the invariant differentials of E$and$E˜, respectively (see [33, Chap. III] for the invariant differential of an elliptic curve). Vélu’s formulae give a normalized separable isogeny, and its modified form is shown in [6] as below (the form was given much earlier by Elkies [16]):

Modified Vélu’s formulae in [6]

Let E : y² = x³ +ax+b be an elliptic curve over a finite field 𝔽_q of characteristic p ≥ 5. Let ℓ be an odd number, and S a subgroup of E(F‾q) of order ℓ. Set S∗=S∖{OE}. Then a normalized separable isogeny ϕ:E→E˜=E/S of degree ℓ can be written as

(1) ϕ(x,y)=N(x)D(x), yN(x)D(x)′,

where D(x) is the polynomial defined as

D(x)=∏Q∈S∗(x−xQ)=xℓ−1−sxℓ−2+s2xℓ−3−s3xℓ−4+⋯

and N(x) is related to D(x) through the formula

(2) N(x)D(x)=ℓx−s−(3x2+a)D′(x)D(x)−2(x3+ax+b)D′(x)D(x)′.

Here T^′(x) denotes the derivative of a function T(x) and x_Q the x-coordinate of a point Q ∈ S^*. With coefficients s, s₂, s₃ of the polynomial D(x), set v=a(ℓ−1)+3(s2−2s2) and w = 3as +2b(ℓ−1)+5(s³ −3ss₂ +3s₃). Then the Weierstrass equation for E˜ is given by y2=x3+a˜x+b˜ with a˜=a−5v$and$b˜=b−7w. (Fast algorithms are shown also in [6] to compute the isogeny ϕ.)

Kernel polynomials

Partition the set S^* into two parts S⁺ and S⁻ such that S∗=S+∪S−andS−={−P:P∈S+}. We consider the kernel polynomial

(3) F(x)=∏Q∈S+(x−xQ)=xd+txd−1+t2xd−2+⋯+td

with d=ℓ−12, and it satisfies D(x) = F(x)². Schoof [30] studied the relation among coefficients t and t_i’s of F(x). Instead of working with E˜, he works with the isomorphic curve Eˆ:y2=x3+aˆx+bˆ with aˆ=ℓ4a˜ and bˆ=ℓ6b˜. To E, we associate the reduced Weierstrass ℘-function by ℘(z)=1z2+∑k=1∞ckz2k with c1=−a5 , c2=−b7 and ck=3(k−2)(2k+3)∑j=1k−2cjck−1−j for k ≥ 3. Similarly, the function ℘ˆ for Eˆ is defined using â and b̂. Then the polynomial F(x) satisfies

(4) zℓ−1F(℘(z))=exp−12tz2−∑k=1∞cˆk−ℓck(2k+1)(2k+2)z2k+2.

(This is by reduction of relations over C, and the characteristic p must be large enough to hold over 𝔽_q.) From this equation, we can represent every coefficient tiwitht,ck's andcˆk's, and with t,a,b,a˜andb˜ since c_k’s and cˆk's are obtained from a,b,a˜andb˜. (See [4, Chap. VII] for the first few coefficients of F(x).) In particular, every t_i can be represented as an element of the multivariate polynomial ring Fq[t,a,b,a˜,b˜] when we regard all t,a,b,a˜,b˜ as variables.

2.3 Modular polynomials

For every integer ℓ ≥ 2, there exists the modular polynomial Φ_ℓ(X, Y) ∈ ℤ[X, Y] to parameterize pairs of elliptic curves with a cyclic isogeny of degree ℓ between them. (See [32, Exercise 2.18, Chap. II].) For two elliptic curves EandE˜ over a finite field 𝔽_q, there is an isogeny of degree ℓ from EtoE˜ with cyclic kernel if and only if Φℓ(j(E),j(E˜))=0. Every modular polynomial Φ_ℓ(X, Y) is symmetric in each variable, and its integer coefficients become very large as ℓ increases. In particular, for a prime ℓ, the modular polynomial Φ_ℓ(X, Y) is equal to the form

(5) Xℓ+1−XℓYℓ+Yℓ+1+∑i,j≤ℓ,i+j<2ℓaijXiYj

with a_ij ∈ ℤ, since there are precisely ℓ + 1 subgroups of the ℓ-torsion group of an elliptic curve E. (Such each subgroup corresponds to an isogenous curve of degree ℓ with a j-invariant, which is a zero of the polynomial Φ_ℓ(X, j(E)).)

3 Algebraic approaches for isogeny problems

In this section, we propose several algebraic approaches. Before presenting our approaches, let us define our setting of isogeny problem as below:

Problem 1

Let ℓ=ℓ0e be a power of a small odd prime ℓ₀. Suppose that there exists an isogeny of degree ℓ between elliptic curves EandE˜ over a finite field 𝔽_q. Given ℓ,j=j(E)andȷ˜=j(E˜), find an isogeny ϕ:E→E˜ of degree ℓ.

Compared to the general isogeny problem in [20], the isogeny degree ℓ is restricted to a prime power and it is given in our setting.Note that a solution of ϕ is not unique in general for the above problem. (In particular, the kernel of ϕ is not necessarily cyclic.) There are also many (compact) representations of a solution ϕ such as a chain of isogenies of low degrees, a sequence of j-invariants of intermediate curves, and a path in an isogeny graph between EandE~. The meet-in-the-middle approach is a practical way to solve isogeny problems. For the above isogeny problem, it builds two trees of isogenies of prime degree ℓ₀ from both the sides of E and E˜, respectively, and it finds a collision between the two trees to find the shortest path from E to E˜. While the meet-in-the-middle approach is a generic way for solving graph problems, we shall propose two algebraic approaches for solving the above isogeny problem.

3.1 First approach using modular polynomials

In this subsection, we present our first approach. Consider a chain of isogenies ϕk of prime degree ℓ₀ from E to E˜ as

E⟶ϕ1E1⟶ϕ2E2ϕ3⟶⋯⟶ϕe−1Ee−1⟶ϕeE˜.

Let j_k denote the j-invariant of every elliptic curve E_k for 1 ≤ k ≤ e−1. In this approach, we regard j-invariants j_k’s as variables, and consider a system of equations using modular polynomials (cf., recall that j and ȷ˜ are elements of 𝔽_q.)

(6) Φℓ0(j,j1)=0,Φℓ0(jk,jk+1)=0(1≤k≤e−2),Φℓ0(je−1,ȷ˜)=0.

A solution of this system gives all j-invariants j_k’s of intermediate curves E_k.

Here we propose a method to solve the system (6) efficiently by Gröbner basis algorithms. (See textbooks [2, 13] for Gröbner basis computation.) We assume that the exponent e of the isogeny degree ℓ is even with e=2e0 for a positive integer e₀ for simple discussion. We divide the system (6) into two parts like the meet-in-the-middle approach. In terms of Gröbner basis computation, we consider two ideals in different multivariate polynomial rings

I=Φℓ0(j,j1),Φℓ0(j1,j2),…,Φℓ0(je0−1,je0)⊂Fq[j1,…,je0],

I˜=Φℓ0(je0,je0+1),Φℓ0(je0+1,je0+2),…,Φℓ0(je−1,ȷ˜)⊂Fq[je0,…,je−1].

Both ideals I and Ĩ are zero-dimensional since j,ȷ˜∈Fq. In particular, the dimensions of 𝔽_q-vector spaces Fq[j1,…,je0]/I and Fq[je0,…,je−1]/I˜ are both at most (ℓ0+1)e0 due to the form of the modular polynomial (see Equation (5)). Moreover, the above generators give a Gröbner basis for the ideal I (resp., the ideal Ĩ) with the lex term order with respect to je0≻⋯≻j2≻j1 (resp., je0≻⋯≻je−2≻je−1). Then we can efficiently compute minimal polynomials g and g˜ of the variable j_e0 with respect to ideals I and Ĩ, respectively, by using the FGLM algorithm [17]. (In this case, simple linear algebra might be more efficient since degrees of g and g˜ are known. See Remark 3.1 and Appendix A.1 below for details.) By the GCD computation over the univariate polynomial ring Fq[je0], we obtain a common root of two minimal polynomials gandg˜. Such a common root gives a solution of je0.

Once a solution of j_e0 is found, the isogeny problem is divided into two isogeny problems of smaller degree ℓ0e0=ℓ i.e., a divide-and-conquer strategy). By repeating this procedure, we can solve the whole isogeny problem.

Remark 3.1

Set g1=Φℓ0(j,j1), and let g_k denote the minimal polynomial of the variable j _k with respect to the ideal ⟨Φℓ0(j,j1),…,Φℓ0(jk−1,jk)⟩ for every 2 ≤ k ≤ e₀. Putting Gk(X):=Φℓ0k(j,X), we have generically

gk(jk)=Gk−2(jk)Gk(jk)

for 3 ≤ k ≤ e₀. (We verified by experiments that it holds in most cases.) For our target minimal polynomial g=ge0, it seems the best in performance to compute g_k from g_k₋₁ and Φℓ0(jk−1,jk) recursively for 2 ≤ k ≤ e₀, see Appendix A.1. In a similar way, we can compute another minimal polynomial ˜g efficiently.

The time complexity of the first approach shall be analyzed in Appendix A.1. The complexity depends on an algorithm for computing minimal polynomials, such as the FGLM algorithm and the GCD computation. Note that our first approach is not better than the meet-in-the-middle approach in time complexity.

Possible improvement (3-section method)

We introduce a possible improvement for the first approach. Our idea is to divide the system (6) into three parts. (cf., The original strategy is regarded as the “2-section method”.) With two parameters e₁ and e₂ satisfying 1 < e₁ < e₀ < e₂ < e and e₁ ≈ e − e₂, consider two ideals in different multivariate polynomial rings

I[1:e1]=Φℓ0(j,j1),Φℓ0(j1,j2),…,Φℓ0(je1−1,je1),

I˜[e2:e−1]=Φℓ0(je2,je2+1),Φℓ0(je2+1,je2+2),…,Φℓ0(je−1,ȷ˜).

As in the 2-section method, we use the lex term order with je1≻⋯≻j2≻j1 (resp., je2≻je2+1≻⋯≻ j_e₋₁) for the zero-dimensional ideal I[1:e1] (resp., I˜[e2:e−1]). Then the ideal I[1:e1] (resp., I˜[e2:e−1]) includes a polynomial g(je1) (resp., g˜(je2)) such that the set of its roots contains the roots of Φℓ′(j,je1) of level ℓ′=ℓ0e1 (resp., Φℓ˜′(je2,ȷ˜) of level ℓ˜′=ℓ0e−e2). Namely, the polynomials gandg˜ are minimal polynomials for zero-dimensional ideals I[1:e1] and I˜[e2:e−1], respectively. We then consider a new ideal

J[e1:e2]=g(je1),Φℓ0(je1,je1+1),…,Φℓ0(je2−1,je2),g˜(je2).

For this zero-dimensional ideal, we use the grevlex term order with je1≺je2≺je1+1≺je2−1≺⋯≺je0 to find intermediate j-invariants from je1 to je2. With these j-invariants, we obtain the other j-invariants j₁, . . . , je1−1 and je2+1,…,je−1 as in the 2-section method. The time complexity of this 3-section method might be improved so that it is better than that of the 2-section method.

3.2 Second approach using kernel polynomials

In this subsection, we present our second approach for solving the isogeny problem of prime power degree ℓ=ℓ0e. As in the first approach, we assume that the exponent e of the isogeny degree is even with e = 2e₀ for simple discussion. In this approach, we take an intermediate elliptic curve E₀ between EandE˜ such that there exist two normalized isogenies φandφ˜ of same degree ℓ′=ℓ0e0 as

E⟶φE0⟵φ˜E˜

Let a,b(resp.,a˜,b˜) denote Weierstrass coefficients in 𝔽_q defining the initial curve E (resp., the final curve E˜). These coefficients can be chosen with two j-invariants jandȷ˜ofEandE˜. In contrast, we use two variables a₀, b₀ as Weierstrass coefficients defining the intermediate curve E₀. (That is, .) E0:y2=x3+a0x+b0.)

Like Equation (3), consider the kernel polynomial

F(x)=xm+txm−1+t2xm−2+⋯+tm,

obtained from the isogeny φ:E→E0 with m=ℓ′−12. Now we regard t as an additional variable. Then we can represent the other coefficients t₂, . . . , t_m of F(x) as elements of the multivariate polynomial ring R = 𝔽_q[t, a₀, b₀] by Schoof’s work [30], described in Subsection 2.2 (recall a, b ∈ 𝔽_q). In other words, we regard the polynomial F(x) as an element of R[x]. Furthermore, like Equation (1), we can reconstruct the isogeny φ:E→E0 with rational functions over R[x] as

φ(x,y)=N(x)D(x),yN(x)D(x)′

for any point (x, y) ∈ E, by letting D(x) = F(x)². Here the rational function T(x)=N(x)D(x)∈R(x) is given by Equation (2). Since φ(x, y) ∈ E₀ for any point (x, y) ∈ E, we clearly have the relation

(7) (x3+ax+b)T′(x)2=T(x)3+a0T(x)+b0.

By expanding this relation with respect to x, we obtain a set of equations S defined over R. (Every coefficient of xⁱ in (7) corresponds to an equation in S.)

Similarly, we consider the kernel polynomial

F˜(x˜)=x˜m+t˜x˜m−1+t˜2x˜m−1+⋯+t˜m

obtained from another isogeny φ˜:E˜→E0. By regarding t˜ as another variable, we can also regard F˜(x˜) as an element of the polynomial ring R˜[x˜]withR˜=Fq[t˜,a0,b0]$. We can also reconstruct the isogeny φ˜ as φ˜(x˜,y˜)=(T˜(x˜),y˜T˜′(x˜)) for any point (x˜,y˜)∈E˜, for some rational function T˜(x˜)∈R˜(x˜). Similarly to the previous paragraph, the relation (x˜3+a˜x˜+b˜)T˜′(x˜)2=T˜(x˜)3+a0T˜(x˜)+b0 gives another set of equations S˜ defined over the ring R˜.

Finally, we solve the system of equations in the union set S∪S˜ by using Gröbner basis algorithms over the multivariate polynomial ring Fq[t,t˜,a0,b0]. (We used the grevlex term order with t≻t˜≻a0≻b0 in our experiments.) A solution of (t,t˜,a0,b0) determines the Weierstrass equation for the intermediate curve E₀, and also two kernels of isogenies φandφ~.

Remark 3.2

While the first approach requires many variables as the exponent e increases, this approach always requires four variables t,t˜,a0,b0. On the other hand, as e increases, total degrees of equations become large in this approach, while total degrees do not change in the first approach due to the use of the modular polynomial of same level ℓ₀. In contrast, one can consider several variants of this approach. As the simplest variant, we directly consider the kernel polynomial of the isogeny ϕ:E→E˜ (without taking any intermediate curve between EandE~). This variant requires only one variable, but total degrees of equations become very large as e increases. On the other hand, we consider multiple intermediate curves as a variant. This variant requires many variables as the number of intermediate curves increases. As a result, the meet-in-the-middle type described in this subsection seems the best in performance from our preliminary experiments.

As in the first approach, the time complexity of the second approach shall be analyzed in Appendix A.2. In particular, the analysis shows that the second approach can never be asymptotically faster than the first one.

4 Experiments

In this section, we report experimental results of our algebraic approaches for solving the isogeny problem of prime power degrees (see Problem 1 for the problem).

4.1 Experiment parameters

For Problem 1, we fix parameters p = 2²⁵⁰ · 3¹⁵⁹ − 1, q = p² and E : y² = x³ + x, extracted from SIKE-p503 parameters [21]. (The extension field 𝔽_p² is represented as 𝔽_p[z]/(z²+1).) The initial curve E is a supersingular elliptic curve defined over 𝔽_p² having #E(Fp2)=(p+1)2=(2250⋅3159)2 and j = j(E) = 1728. We also take 3-powers ℓ = 3^e (i.e., ℓ₀ = 3) as isogeny degrees for even exponents e=2e0. (cf., SIKE uses a combination of isogenies of degrees 2 and 3.) We follow [21] to generate the final supersingular curve E˜, isogenous to E of degree ℓ.

4.2 Implementation details

Here we describe details of implementation for our algebraic approaches and the meet-in-the-middle approach with Magma [5], a computational algebra system.

For our first approach by 2-section and 3-section methods, we used a combination of the modular polynomials Φ_N(X, Y) for N = 3, 3², 3³, which are pre-computed in Magma, in order to obtain the minimal polynomials g,g˜. For example, for ℓ = 310, we computed Gröbner bases for I=⟨Φ33(j,j3),Φ32(j3,j5)⟩, I˜=⟨Φ32(j5,j7),Φ33(j7,ȷ˜)⟩ with the Magma command GroebnerBasis to obtain g, g˜∈Fp2[j5], then computed the GCD of g, ˜g with the Magma commands GCD for the 2-section method.

For our second approach, we used Equation (4) to represent two kernel polynomials F(x)andF˜(x˜) of isogenies φ and φ˜, described in Subsection 3.2, as polynomials over Fq[t,t˜,a0,b0]. (As an alternative of (4), fast algorithms are introduced in [6] for computing the kernel polynomial of an isogeny, but they require very fast exponentiation.) We also used the grevlex term order with t≻t˜≻a0≻b0 to find a solution (t,t˜,a0,b0) from the union set of equations S∪S˜.

For the meet-in-the-middle approach, we construct two sets JandJ˜ of sequences (j, j₁, . . . , j_e0 ) and (ȷ˜,ȷ˜1,…,ȷ˜e0) of j-invariants of elliptic curves EkandE˜k, respectively. (Recall j = j(E) and ȷ˜=j(E˜).) Here E_k₋₁ and E_k (resp., E˜k−1andE˜k) are isogenous of degree 3 for every 1 ≤ k ≤ e₀, where we set E0=E(resp.,E˜0=E˜) for convenience. To construct sequences in JandJ˜, we use the modular polynomial of level 3. For example, we add each solution of Φ3(jk−1,x) to a sequence (j, j₁, . . . , j_k₋₁) of length k. (We used the Magma command Roots to find a solution.) In constructing such sequences, we remove a sequence whose ending point already appeared in the other sequences, in order to reduce the sizes of two sets JandJ˜. In our experiments, it terminates when we find a pair of two sequences of JandJ˜ satisfying je0=ȷ˜e0 (i.e., a collision).

4.3 Experimental results

In Table 1, we summarize average running times of our two approaches and the meet-in-the-middle approach for solving the isogeny problem of degrees ℓ = 3^e with even e from e = 6 up to 14. Specifically, we measured the running time of every approach until it finds the j-invariant or the Weierstrass coefficients of an intermediate curve between two isogenous curves EandE˜. We also experimented 5 times for every parameter set. All the experiments were performed using Magma 2.24-5 on 4.20 GHz Intel Core i7 CPU with 16 GByte memory. From Table 1, the first approach is the fastest for isogeny degrees up to ℓ = 310. In contrast, our second approach is costly and it did not terminate in one day for degrees larger than ℓ = 38. For degrees larger than ℓ = 310, the meet-in-the-middle approach is faster than our algebraic approaches. With respect to the memory usage, our first approach by the 2-section method requires about 64, 221 and 526MByte for ℓ = 310, 3¹² and 314, respectively, while the meet-in-the-middle approach requires about 24, 26 and 32 MByte for the same isogeny degrees.

Table 1

Average running times (seconds) of our two approaches and the meet-in-the-middle approach for the isogeny problem of degrees ℓ = 3^e on supersingular elliptic curves over 𝔽_p2 with 503-bit prime p = 2²⁵⁰ · 3¹⁵⁹ − 1 (These parameters are from SIKE-p503 [21], and the initial curve E is given by y² = x³ + x over 𝔽_p2)

Isogeny degrees	Our first approach		Our second approach	Meet-in-the-middle approach
ℓ = 3^e	2-section	3-section
3⁶ = 729	0.11	0.15	336.42	2.93
3⁸ = 6561	0.19	0.45	> 1 day	9.61
3¹⁰ = 59049	7.98	5.63	–	31.09
3¹² = 531441	287.77	159.39	–	96.75
3¹⁴ = 4782969	5071.16	2725.01	–	292.82

5 Concluding remarks

In this paper, we proposed two algebraic approaches for solving isogeny problems of prime power degrees. The first one is a straightforward way using modular polynomials of small prime levels. The second one is a more complex way, which uses kernel polynomials of isogenies based on [30] and reconstructs the isogenies with Vélu’s formulae [6, 39]. From the analysis in Appendix A, our approaches are not asymptotically faster than the meet-in-the-middle approach. However, our experiments showed that the first approach is faster than the meet-in-the-middle approach for isogeny degrees up to ℓ = 3¹⁰ over supersingular elliptic curves of SIKE-p503 [21] with a single classical computer. On the one hand, our first approach is applicable in the collision search step of the meet-in-the-middle approach. Such combination could make it faster in practice and reduce the memory size of the meet-in-the-middle approach. As future work, we would like to use the combination for solving isogeny problems of large degrees.

Acknowledgement

This work was supported by JSPS KAKENHI Grant Numbers 18H05836, 19K21026, and 19K22847, Japan

References

[1] Gora Adj, Daniel Cervantes-Vázquez, Jesús-Javier Chi-Domínguez, Alfred Menezes and Francisco Rodríguez-Henríquez, On the cost of computing isogenies between supersingular elliptic curves, in: Selected Areas in Cryptography–SAC 2018, Lecture Notes in Computer Science 11349, Springer, pp. 322–343, 2018.10.1007/978-3-030-10970-7_15Search in Google Scholar

[2] Thomas Becker and Volker Weispfenning, Gröbner bases, Graduate Texts in Mathematics 141, Springer, 1993.10.1007/978-1-4612-0913-3_5Search in Google Scholar

[3] Jean-François Biasse, David Jao and Anirudh Sankar, A quantum algorithm for computing isogenies between supersingular elliptic curves, in: Progress in Cryptology–INDOCRYPT 2014, Lecture Notes in Computer Science 8885, Springer, pp. 428–442, 2014.10.1007/978-3-319-13039-2_25Search in Google Scholar

[4] Ian F Blake, Gadiel Seroussi and Nigel Smart, Elliptic curves in cryptography, 265, Cambridge university press, 1999.10.1017/CBO9781107360211Search in Google Scholar

[5] Wieb Bosma, John Cannon and Catherine Playoust, The Magma algebra system. I. The user language, J. Symbolic Comput. 24 (1997), 235–265, Computational algebra and number theory (London, 1993).10.1006/jsco.1996.0125Search in Google Scholar

[6] Alin Bostan, François Morain, Bruno Salvy and Éric Schost, Fast algorithms for computing isogenies between elliptic curves, Mathematics of Computation 77 (2008), 1755–1778.10.1090/S0025-5718-08-02066-8Search in Google Scholar

[7] Alessio Caminata and Elisa Gorla, Solvingmultivariate polynomial systems and an invariant from commutative algebra, arXiv preprint arXiv:1706.06319 (2017).Search in Google Scholar

[8] Wouter Castryck, Tanja Lange, Chloe Martindale, Lorenz Panny and Joost Renes, CSIDH: An efficient post-quantum commutative group action, in: Advances in Cryptology–ASIACRYPT 2018, Lecture Notes in Computer Science 11274, Springer, pp. 395–427, 2018.10.1007/978-3-030-03332-3_15Search in Google Scholar

[9] Denis X Charles, Kristin E Lauter and Eyal Z Goren, Cryptographic hash functions from expander graphs, Journal of Cryptology 22 (2009), 93–113.10.1007/s00145-007-9002-xSearch in Google Scholar

[10] Andrew Childs, David Jao and Vladimir Soukharev, Constructing elliptic curve isogenies in quantum subexponential time, Journal of Mathematical Cryptology 8 (2014), 1–29.10.1515/jmc-2012-0016Search in Google Scholar

[11] Anamaria Costache, Brooke Feigon, Kristin Lauter, Maike Massierer and Anna Puskás, Ramanujan graphs in cryptography, IACR ePrint 2018/593 (2018).10.1007/978-3-030-19478-9_1Search in Google Scholar

[12] Craig Costello, Patrick Longa, Michael Naehrig, Joost Renes and Fernando Virdia, Improved Classical Cryptanalysis of the Computational Supersingular Isogeny Problem, IACR ePrint 2019/298 (2019).Search in Google Scholar

[13] David Cox, John Little and Donal O’Shea, Ideals, varieties, and algorithms: An introduction to computational algebraic geometry and commutative algebra, Springer Science & Business Media, 2013.Search in Google Scholar

[14] Luca De Feo, David Jao and Jérôme Plût, Towards quantum-resistant cryptosystems from supersingular elliptic curve isogenies, Journal of Mathematical Cryptology 8 (2014), 209–247.10.1515/jmc-2012-0015Search in Google Scholar

[15] Kirsten Eisenträger, Sean Hallgren, Kristin Lauter, Travis Morrison and Christophe Petit, Supersingular isogeny graphs and endomorphism rings: reductions and solutions, in: Advances in Cryptology–EUROCRYPT 2018, Lecture Notes in Computer Science 10822, Springer, pp. 329–368, 2018.10.1007/978-3-319-78372-7_11Search in Google Scholar

[16] Noam D Elkies et al., Elliptic and modular curves over finite fields and related computational issues, AMS IP STUDIES IN ADVANCED MATHEMATICS 7 (1998), 21–76.10.1090/amsip/007/03Search in Google Scholar

[17] Jean-Charles Faugere, Patrizia Gianni, Daniel Lazard and Teo Mora, Efficient computation of zero-dimensional Gröbner bases by change of ordering, Journal of Symbolic Computation 16 (1993), 329–344.10.1006/jsco.1993.1051Search in Google Scholar

[18] Steven D Galbraith, Constructing isogenies between elliptic curves over finite fields, LMS Journal of Computation and Mathematics 2 (1999), 118–138.10.1112/S1461157000000097Search in Google Scholar

[19] Steven D Galbraith, Christophe Petit, Barak Shani and Yan Bo Ti, On the security of supersingular isogeny cryptosystems, in: Advances in Cryptology–ASIACRYPT 2016, Lecture Notes in Computer Science 10031, Springer, pp. 63–91, 2016.10.1007/978-3-662-53887-6_3Search in Google Scholar

[20] Steven D Galbraith and Frederik Vercauteren, Computational problems in supersingular elliptic curve isogenies, Quantum Information Processing 17 (2018), 265.10.1007/s11128-018-2023-6Search in Google Scholar

[21] D Jao, R Azarderakhsh, M Campagna, C Costello, L DeFeo, B Hess, A Jalali, B Koziel, B LaMacchia, P Longa et al., SIKE: Supersingular isogeny key encapsulation. Submission to the NIST Standardization Process on Post-Quantum Cryptography, 2017.Search in Google Scholar

[22] David Jao and Luca De Feo, Towards quantum-resistant cryptosystems from supersingular elliptic curve isogenies, in: Post-Quantum Cryptography–PQCrypto 2011, Lecture Notes in Computer Science 7071, Springer, pp. 19–34, 2011.10.1007/978-3-642-25405-5_2Search in Google Scholar

[23] David Jao and Vladimir Soukharev, Isogeny-based quantum-resistant undeniable signatures, in: Post-Quantum Cryptography–PQCrypt 2014, Lecture Notes in Computer Science 8772, Springer, pp. 160–179, 2014.10.1007/978-3-319-11659-4_10Search in Google Scholar

[24] Samuel Jaques and John M Schanck, Quantum cryptanalysis in the RAM model: Claw-finding attacks on SIKE, IACR Cryptology ePrint Archive 2019/103 (2019).10.1007/978-3-030-26948-7_2Search in Google Scholar

[25] Neal Koblitz, Elliptic curve cryptosystems, Mathematics of Computation 48 (1987), 203–209.10.1090/S0025-5718-1987-0866109-5Search in Google Scholar

[26] Chloe Martindale and Lorenz Panny, How to not break SIDH, IACR ePrint 2019/558 (2019).Search in Google Scholar

[27] Victor S Miller, Use of elliptic curves in cryptography, in: Advances in Cryptology–CRYPTO 1985, Lecture Notes in Computer Science 218, pp. 417–426, Springer, 1985.10.1007/3-540-39799-X_31Search in Google Scholar

[28] National Institute of Standards and Technology (NIST), NISTIR 8240: Status Report on the First Round of the NIST Post-Quantum Cryptography Standardization Process.Search in Google Scholar

[29] Alexander Rostovtsev and Anton Stolbunov, Public-key cryptosystem based on isogenies, IACR Cryptology ePrint Archive 2006/145 (2006).Search in Google Scholar

[30] René Schoof, Counting points on elliptic curves over finite fields, Journal de théorie des nombres de Bordeaux 7 (1995), 219–254.10.5802/jtnb.142Search in Google Scholar

[31] Peter W. Shor, Algorithms for quantum computation: Discrete logarithms and factoring, in: Symposium on Foundations of Computer Science (FOCS 1994), pp. 124–134, IEEE, 1994.Search in Google Scholar

[32] Joseph H Silverman, Advanced topics in the arithmetic of elliptic curves, Graduate Texts in Mathematics 151, Springer-Verlag New York, 1994.10.1007/978-1-4612-0851-8Search in Google Scholar

[33] Joseph H Silverman, The arithmetic of elliptic curves, second ed, Graduate Texts in Mathematics 106, Springer Science & Business Media, 2009.10.1007/978-0-387-09494-6Search in Google Scholar

[34] Anton Stolbunov, Constructing public-key cryptographic schemes based on class group action on a set of isogenous elliptic curves, Advances in Mathematics of Communications 4 (2010), 215–235.10.3934/amc.2010.4.215Search in Google Scholar

[35] Xi Sun, Haibo Tian and Yumin Wang, Toward quantum-resistant strong designated verifier signature from isogenies, in: Intelligent Networking and Collaborative Systems–INCoS 2012, IEEE, pp. 292–296, 2012.10.1109/iNCoS.2012.70Search in Google Scholar

[36] Seiichiro Tani, Claw finding algorithms using quantum walk, Theoretical Computer Science 410 (2009), 5285–5297.10.1007/978-3-540-74456-6_48Search in Google Scholar

[37] John Tate, Endomorphisms of abelian varieties over finite fields, Inventiones mathematicae 2 (1966), 134–144.10.1007/BF01404549Search in Google Scholar

[38] Paul C Van Oorschot and Michael JWiener, Parallel collision search with cryptanalytic applications, Journal of Cryptology 12 (1999), 1–28.10.1007/PL00003816Search in Google Scholar

[39] Jacques Vélu, Isogénies entre courbes elliptiques, CR Acad. Sci. Paris Sér. AB 273 (1971), 238–241.Search in Google Scholar

A Time complexity analysis

In this appendix, we analyze time complexities of our algebraic approaches. Here we use the same notation as in Section 3. In this complexity analysis, we assume that ℓ and ℓ₀ are constants, and we set e and e₀ as asymptotic parameters to estimate time complexities.

A.1 For the first approach

The first approach proceeds with the below two steps:

Step 1. Compute minimal polynomials gandg˜ of the variable je0 with respect to ideals IandI˜, respectively.

Step 2. By computing the GCD of gandg˜ over the univariate polynomial ring Fq[je0], we obtain a common root in 𝔽_q of the minimal polynomials.

There are several methods to compute Step 1, such as the FGLM algorithm and the GCD computation. In the below, we estimate costs of the first approach using several different methods for Step 1.

Using the FGLM algorithm for Step 1

Applying the FGLM algorithm directly to the whole ideal

The cost of Step 1 by using the FGLM algorithm for the ideal I is estimated as O e0((ℓ0+1)e0)3=Oe0(ℓ0+1)3e0, since dimFqFq[j1,…,je0]/I≤(ℓ0+1)e0 and the number of variables is e₀. Similarly, the FGLM algorithm for another ideal Ĩ requires the same cost. Hence the complexity of Step 1 is

O(2e0(ℓ0+1)3e0)=O(e0(ℓ0+1)3e0)=O(e(ℓ0+1)3e/2)=O(eℓ032e).

Using the normal form computation and simple linear algebra only

Instead of directly using of the FGLM algorithm, we can compute minimal polynomials gandg˜ based on the below well-known fact:

Lemma A.1

Let K be a field, I a zero-dimensional ideal of S = K [X₁, . . . , X_n], and G a Gröbner basis of I with respect to a fixed term order ≻. Let f ∈ S and g(t) ∈ K [t] with leading coefficient 1. We denote by m_f (t) the minimal polynomial of f with respect to the ideal I. Then g(t) = m_f (t) if and only if g(t) satisfies the following two conditions:

NFI(g(f))=0, where NF_I(g(f )) denotes the normal form of g(f) ∈ S with respect to G.
For any h(t) ∈ K[t] with NFI(h(f))=0, we have deg(g) ≤ deg(h).

Recall from Remark 3.1 that we have generically g=Gk−2Gkfor3≤k≤e0, where we set Gk(jk):=Φℓk(j,jk) for 1 ≤ k ≤ e₀. (Until the end of this appendix, we assume this generic condition.) Since the degree of g equals to d:=deg(g)=deg(Ge0−2)+deg(Ge0)=(ℓ0+1)ℓ0e0−3+(ℓ0+1)ℓ0e0−1=O(ℓ0e0), it suffices to compute NF_I(X^k) for 0≤k≤d=O(ℓ0e0) with X:=je0, and solve the linear equation NFI(Xd)+∑k=0d−1ckNFI(Xk)=0 over 𝔽_q, and then we have g=∑k=0dckXk by Lemma A.1. The cost of computing each normal form NF_I(X^k) is O((ℓ0+1)2e0)=O(ℓ02e0). Hence the cost of computing NFI(Xk) for 0 ≤ k ≤ d is O((ℓ0e0ℓ02e0))=O(ℓ03e0). The cost of solving the linear equation is not higher than O(ℓ03e0) since it is O(d3)=O(ℓ03e0) if one uses the Gaussian elimination naively. Similarly, computing g˜ requires O(ℓ03e0). Here we estimate that the cost of Step 1 is O(2ℓ03e0)=O(ℓ03e0)=O(ℓ032e). In contrast, the cost of Step 2 mainly depends on the computation of the GCD of g and ̃g over 𝔽_q. Note that the complexity of computing the GCD of two univariate polynomials of degree at most d over 𝔽_q is O(d²) arithmetic over 𝔽_q (if one uses classical polynomial arithmetic). Since both degrees of gandg˜ are equal to (ℓ0+1)ℓ0e0−3+(ℓ0+1)ℓ0e0−1=O(ℓ0e0), we estimate that the cost of Step 2 is O(ℓ02e0)=O(ℓ0e). Summing up, we roughly estimate that the complexity of the first approach is

O(eℓ032e)+O(ℓ0e)=O(eℓ032e)

if the FGLM algorithm is applied. Moreover, the complexity becomes

O(ℓ032e)+O(ℓ0e)=O(ℓ032e)

if only the normal form computation and simple linear algebra is used. For the parameter in our experiments (i.e., ℓ₀ = 3), the cost is O(e033e0) or O(33e0).

Using sequential GCD for Step 1

As mentioned in Remark 3.1, it is practical to compute minimum polynomials g_k together with G_k(X) := Φℓ0k(j,X) sequentially for 2 ≤ k ≤ e₀. (Note the target minimal polynomial g is given by ge0.) The minimal polynomial g_k is the lowest degree polynomial in j _k belonging to the ideal generated by Gk−1(jk−1)=Φℓ0k−1(j,jk−1) and Φℓ0(jk−1,jk), that is

Gk−1(jk−1),Φℓ0(jk−1,jk)∩Fq[jk]=gk(jk).

Once g_k is computed, we obtain G_k as the quotient of g_k divided by G_k₋₂ since g_k(X) = G_k₋₂(X)G_k(X). Here we estimate the cost of computing G_k. Regarding G_k₋₁(j_k₋₁) and Φℓ0(jk−1,jk) as polynomials in j_k₋₁ over 𝔽_q[j_k], the cost of computing G _k can be upper-bounded by that of computing their greatest common divisor with respect to j _k₋₁ by the Euclidean algorithm with pseudo division. This is due to the following reason; In the below steps (1) and (2), the target minimal polynomial can be computed by eliminating j_k₋₁ from Gk−1′(jk−1,jk) and Φℓ0(jk−1,jk). The Euclidean algorithm with pseudo division corresponds to eliminating j_k₋₁ by using leading coefficients as polynomials over 𝔽_q[j_k], whereas the Gröbner basis computation corresponds to eliminating j_k₋₁ by using head terms as polynomials in 𝔽_q[j_k₋₁, j _k] with respect to the lexicographical term order with j_k₋₁ ≻ j_k. In each elimination step of the Euclidean algorithm with pseudo division, the elimination of a term is just the same as computing an S-polynomial (or its multiple by a monomial). To make our estimation simple, we can use so called the subresultant GCD algorithm which is the Euclidean algorithm with pseudo division and principal subresultant coefficient (PSC) and computes the resultant, that is, not the minimal polynomial but the characteristic polynomial (Gk−2(jk))ℓ0Gk(jk) of j _k with respect to ⟨Gk−1(jk−1),Φℓ0(jk−1,jk)⟩.

Now we estimate the cost of computing the resultant of G_k₋₁ and Φℓ0(jk−1,jk) with respect to j_k₋₁ by the subresultant GCD, where we may assume that the remainder sequence is normal. Note that degjk−1(Gk−1(jk−1))= (ℓ0+1)ℓ0k−2=O(ℓ0k−1), and degjk−1(Φℓ0(jk−1,jk))=ℓ0+1. The final target G_k is computed by conducting two procedures below:

Using the division algorithm, compute Gk−1′(jk−1,jk):=Gk−1(jk−1) mod Φℓ0(jk−1,jk), where Φℓ0(jk−1,jk) is a monic polynomial over 𝔽_q[j_k] of degree ℓ₀ + 1 in j_k₋₁. Note that Gk−1′ is a univariate polynomial over 𝔽_q[j_k] in variable j_k₋₁ of degree ≤ ℓ₀. More specifically, proceed with the following (we use simple notation A, B, C for polynomials and N for degrees):

(1-1) Set A←Φℓ0(jk−1,jk), and B ← G_k₋₁.

(1-2) Set N←degjk−1(B), (B), and C ←(the coefficient of jk−1N in B).

(1-3) Compute B−jk−1N−(ℓ0+1) CA, and set B←B−jk−1N−(ℓ0+1) CA.

(1-4) If N ≤ ℓ₀, stop the loop. Otherwise set N←degjk−1 (B) and go back to (1-2).

The resulting polynomial B coincides with the target Gk−1′.
Compute the resultant of Gk−1′(jk−1,jk) and Φℓ0(jk−1,jk) by the subresultant GCD, and then recover G_k from the resultant.

We estimate the cost of the procedure (1). Let n be the number of times the loop executes, and let (B_s , N_s , C_s) for 1 ≤ s ≤ n be tuples of the successive values B_s of B, N_s of N and C_s of C respectively. We denote by M_s the maximum degree of the coefficients of jk−1i for 0 ≤ i ≤ N_s − 1 in B_s. Note that B1=Gk−1∈Fq[jk−1], N1=degjk−1(Gk−1)=O(ℓ0k−1),C1∈Fq∖{0} and M₁ ≤ 0. Recall from (2.5) that

A=Φℓ0(jk−1,jk)=jk−1ℓ0+1−(jkℓ0)jk−1ℓ0+∑i,i′≤ℓ0,i+i′<2ℓ0aii′jk−1ijki′+jkℓ0+1,

and hence all the coefficients of jk−1i for 0 ≤ i ≤ N_s of jk−1Ns−(ℓ0+1)CsA are polynomials over 𝔽_q[j_k] of degrees at most deg(C_s)+ℓ₀+1. From this, the degree of the coefficient of each jki with 0 ≤ i ≤ N_s in Bs+1=Bs−jk−1Ns−(ℓ0+1)CsA is at most max {M_s, deg(C_s)+ℓ₀ +1}, and thus both deg(C_s₊₁) and M_s₊₁ are ≤ max {M_s, deg(C_s)+ℓ₀ +1}. Since M₁ ≤ 0 ≤ deg(C₁)+ℓ₀ +1, successively we have M_s ≤ (s −1)(ℓ₀ +1) and deg(Cs)≤(s−1)(ℓ0+1). Since A and C_s have at most (ℓ0+1)2+2=O((ℓ0+1)2) and deg(C_s) terms respectively, computing C _sA requires O((s−1)(ℓ0+1)3) arithmetic in 𝔽_q, which is clearly dominant in the n-th loop. It follows from n≤deg(B1)−degjk−1(A)+1= O(ℓ0k−1) and ∑s=1n(s−1)(ℓ0+1)3=O(ℓ03n2) that the total cost of the procedure (1) is

O(ℓ03(ℓ0k−1)2)=O(ℓ02k).

By the estimation of the cost of computing a resultant of two polynomials, we estimate the complexity of the procedure (2). The procedure (2) requires O(ℓ02) arithmetic over 𝔽_q[j_k] since both Gk−1′ and Φℓ0(jk−1,jk) have degree ≤ ℓ₀ + 1 = O(ℓ₀). Since the degree of the resultant gˆk=(Gk−2)ℓ0Gkisdˆk:=(ℓ0+1)2ℓ0k−2= O(ℓ0k), we may assume that all arithmetic over 𝔽_q[j_k] in the procedure (2) are arithmetic of two polynomials of degree O(ℓ0ℓ0k)=O(ℓ0k) by considering pseudo division and PSC. Hence the complexity of the procedure (2) is estimated as O(ℓ02(ℓ0k)2)=O(ℓ02k) if one uses classical polynomial arithmetic, and O(ℓ02(ℓ0k)1+ϵ)=O(ℓ0(1+ϵ)k) for some 0 < ϵ < 1 if one uses fast polynomial computation (e.g., fast Fourier transform). Dividing the degree O(ℓ0k) -polynomial gˆk by (Gk−2)ℓ0 provides G_k, and it is done in Oℓ02k (or O(ℓ0(1+ϵ)k)), which is equal to the cost of (2). Summing up, we estimate that the total cost is

O(ℓ04)+O(ℓ06)+⋯+O(ℓ02e0)=O(ℓ02e0)=O(ℓ0e)=O(ℓ)

or O(ℓ0(1+ϵ)e0)=O(ℓ(1+ϵ)/2). For the parameter in our experiments (i.e., ℓ₀ = 3), the cost is O(32e0).

Remark A.1

For the meet-in-the-middle approach, we factorize the univariate polynomial Φ₃(x, j) over 𝔽 p² for the j-invariant j = j(E) of an elliptic curve E to find j-invariants of all elliptic curves 3-isogeneous to E. It costs roughly O(log p), and hence the total cost of the meet-in-the-middle approach is O(3e0logp). On the other hand, the total cost of our first approach is O(32e0) as discussed above, and hence we expect that our first approach could be faster than the meet-in-the-middle approach when log p>3e0. This is a main reason why our first approach is faster than the meet-in-the-middle approach in practice for small exponents e₀.

A.2 For the second approach

We estimate the complexity of our second approach for solving the isogeny problem of prime degree ℓ=ℓ0e with e=2e0. We use the same notation as in Subsection 3.2. To estimate the complexity, we calculate the number of equations in the system, that is, #(S∪S˜), and the maximum value of the total degrees of the equations with respect to t,t˜,a0 and b₀.

First, we bound the total degrees of tiandt˜i with respect to t,t˜,a0 and b₀ for 2 ≤ i ≤ m. The kernel polynomial F (x) for E is written as F=xm+txm−1+∑k=0m−2tm−kxk. For each 2 ≤ i ≤ m, the total degree of t_i with respect to t, a₀ and b₀ is i, and thus it is bounded by m. Similarly, the total degree of t˜i with respect to t˜ a₀ and b₀ is bounded by m. Next, we explicitly write down a polynomial G=∑kckxk∈R[x] such that S={ck:0≤k≤degxG}, where the coefficient ring is R = 𝔽_q[t, a₀, b₀]. Recall that we obtain S by expanding (3.2) with respect to x. It follows from T=ND and D = F ² that we have T′=N′F−2NF′F3, and thus (3.2) is rearranged as

(8) (x3+ax+b)N′F−2NF′F32=N3F6+a0NF2+b0.

Note that

N=(ℓx+2t)F2−2(3x2+a)F′F−4(x3+ax+b)(F′′F−(F′)2)

by (2.2), and hence N is of degree 2m + 1 with respect to x. We also note that each coefficient of N is a polynomial in R = 𝔽_q[t, a₀, b₀] of total degree ≤ 2m + 1 = O(m). Multiplying by F ⁶ the both side of (8), we have

(x3+ax+b)N′F−2NF′2=N3+a0NF4+b0F6

and set

G:=(x3+ax+b)N′F−2NF′2−N3−a0NF4−b0F6,

which is of degree 6m+3 with respect to x. Writing G=∑k=06m+3ckxk with c_k ∈ R for 0 ≤ k ≤ 6m+3, we clearly have S = {c_k : 0 ≤ k ≤ 6m+3}, and the total degree of each cis equal to or less than 6m+3 = O(m). Similarly, we can construct a polynomial G˜=∑kc˜kxk∈R˜[x] with R˜=Fq[t˜,a0,b0] such that S˜={c˜k:0≤k≤degxG˜} with degxG˜=6m+3. The total degree of each ̃c_k is equal to or less than 6m + 3 = O(m).

Here we assume that we use the algorithm F ₄, which is default for computing Gröbner bases in Magma, to solve a system of (non-homogeneous) multivariate equations. Let IS∪S˜ denote the ideal defined by S∪S˜. Let solv. deg(IS∪S˜) denote the solving degree of the ideal IS∪S˜⊂Fq[t,t˜,a0,b0] with respect to the grevlex order. To estimate the complexity of computing a Gröbner basis of IS∪S˜, we use the following proposition proved by Caminata and Gorla:

Proposition A.1 ([7])

Let K be a field, and R = K [X₁, . . . , X_n] the polynomial ring of n variables. Let f₁, . . . , f_ℓ be (not necessarily homogeneous) polynomials in R, and I the ideal generated by f₁, . . . , f_ℓ. We denote by s = solv. deg(I) the solving degree of I, that is, the highest degree of the polynomials involved in the computation of a Gröbner basis of I (we fix the grevlex order for Gröbner bases). Let d_i be the total degree of f_i for 1 ≤ i ≤ ℓ, and

m:=∑i=1ℓn+s−dis−di.

Then the number of operations in K required to compute a Gröbner basis of I is

On+ssmω−1ifm≤n+ss,

Omn+ssω−1otherwise,

where 2 ≤ ω ≤ 3 is the exponent of matrix multiplication.

Since the total degree of each c∈S∪S˜is≈6m+3=O(m), we have

solv. deg(IS∪S˜)≥6m+3=O(m),

where the solving degree is not less than the maximum degree of input polynomials. It also follows from #(S∪S˜)≈2(6m+3)=12m+6=O(m) and

4+solv. deg(IS∪S˜)solv. deg(IS∪S˜)=4+solv. deg(IS∪S˜)4≥O(m4)

that the number of operations in 𝔽_q required to compute a Gröbner basis of IS∪S˜ in the second approach is ≥Om5=O(ℓ′)5=Oℓ05e0 by Proposition A.1, and so is the total complexity of the second approach. For the parameter in our experiments (i.e., ℓ₀ = 3), the complexity is not less than 35e0.

Received: 2019-06-05

Accepted: 2019-07-01

Published Online: 2020-11-17

This work is licensed under the Creative Commons Attribution 4.0 International License.