Radical Isogenies

Isogeny-based cryptography is one of the more promising candidates for post-quantum cryptography and although it is slower than lattice-based cryptography, it has the advantage of smaller key and ciphertext sizes. Isogeny-based protocols can be broadly categorized into two families: SIDH and CRS/CSIDH.

SIDH is a key agreement protocol introduced by Jao and De Feo in 2011 [16]. This protocol is based on random walks in isogeny graphs of supersingular elliptic curves E over \(\mathbb {F}_{p^2}\), and is reminiscent of the CGL hash function due to Charles, Goren and Lauter from 2009 [8]. The prime p is chosen such that the torsion subgroups \(E[2^n]\) and \(E[3^m]\) are defined over \(\mathbb {F}_{p^2}\), for large exponents \(n,\,m\). The random walks then correspond to choosing a random point P in \(E[2^n]\) or \(E[3^m]\) and constructing the isogeny with kernel \(\langle P \rangle \), as a composition of isogenies of degree 2 respectively 3.

CRS/CSIDH [7] takes a different approach and computes an action of the ideal-class group \({\text {cl}}(\mathcal {O})\) of some order \(\mathcal {O}\) in an imaginary quadratic field on the set \(\mathcal {E}\ell \ell _p(\mathcal {O}, t)\) of elliptic curves over a prime field \(\mathbb {F}_p\) with \(\mathbb {F}_p\)-rational endomorphism ring \(\mathcal {O}\) and trace of Frobenius t. The idea of using this class group action in cryptography was independently proposed by Couveignes [11] and Rostovtsev-Stolbunov [22] for ordinary elliptic curves. In [7] this idea was ported to the supersingular case, resulting in a speed-up of several orders of magnitude. The computation of the class group action boils down to computing chains of \(\ell \)-isogenies for many small primes \(\ell \), e.g., for CSIDH-512, \(\ell \) ranges from 3 to 587. This is in stark contrast with SIDH where only 2- and 3-isogenies are used.

In the CSIDH setting, computing an \(\ell \)-isogeny \(\varphi \) from an elliptic curve \(E/\mathbb {F}_p\) consists of two steps: first, a generator P of the kernel of \(\varphi \) is computed, i.e. an \(\mathbb {F}_p\)-rational point of order \(\ell \), and secondly, given P, an equation for the isogenous curve \(E/\langle P \rangle \) is determined.

The most basic approach to solve the first step is to generate a random point \(Q \in E(\mathbb {F}_p)\) and to multiply this by the cofactor \(\#E(\mathbb {F}_p) / \ell \). Generating a random point is essentially a square root computation at a cost of about \(1.5 \log p\) multiplications in \(\mathbb {F}_p\), and the multiplication by the cofactor can be done using the Montgomery ladder [2] and takes roughly \(11\log p\) multiplications in \(\mathbb {F}_p\). Generating a point of order \(\ell \) is thus a costly operation, even further exacerbated by the fact that multiplication by the cofactor results in the point at infinity \(\mathcal {O}_E\) with probability \(1/\ell \), which is non-negligible for small \(\ell \). Note that this also makes the algorithm non-deterministic, negatively affecting constant time implementations. The cost of generating \(\ell \)-torsion points from scratch can be mitigated somewhat by considering a chain of \(\ell _i\)-isogenies for many different primes \(\ell _i\). Instead of sampling an \(\ell _i\)-torsion point for every \(\ell _i\)-isogeny separately, it is cheaper to sample an \(\prod _{i=1}^k \ell _i\)-torsion point and push it through the isogeny to create a chain of isogenies of respective degrees \(\ell _1,\ell _2,\ldots ,\ell _k\), multiplying this point with a cofactor that gets smaller in each iteration.

The second step is typically carried out using some form of Vélu’s formulae [28], which compute the coefficients of \(E/\langle P \rangle \) from the coefficients of E and the coordinates of the scalar multiples of P. Vélu’s formulae can also be used to compute the image \(\varphi (Q)\) of any point Q under the isogeny. The original implementation of CSIDH uses these formulae on elliptic curves in Montgomery form [7, 21], and requires \(O(\ell )\) arithmetic operations in \(\mathbb {F}_p\) per \(\ell \)-isogeny. Since then many optimizations to CSIDH have been proposed, such as:

• using different forms of elliptic curves, e.g. twisted Edwards curves [18, 19] and Hessian curves [12, 14];

• adapting Vélu’s formulae to only require \(\widetilde{O}(\sqrt{\ell })\) operations in \(\mathbb {F}_p\) [1] instead of \(O(\ell )\);

• changing CSIDH into CSURF to allow the use of very efficient 2-isogenies [6],

• lowering the number of \(\ell \)-isogenies that has to be computed for each \(\ell \) [9, 20].

A number of alternative approaches have been considered that avoid the generation of \(\ell \)-torsion points altogether, e.g. by using modular polynomials [3, 13] or division polynomials [3]. This leads to deterministic algorithms which can outperform the above method using Vélu’s formulae for small \(\ell \). Highly optimized approaches exist for 2-isogenies [6] and 3-isogenies [12, 14], where the speed-up stems from two ingredients: firstly, an elliptic curve model is chosen that is nicely adapted to 2-torsion (a variant of Montgomery curves) resp. 3-torsion (Hessian curves). The second and main ingredient however is that the coefficients of \(E/\langle P \rangle \) can be expressed in terms of the coefficients of E and a single radical of a simple algebraic expression in the coefficients of E. This radical is a square root for 2-isogenies and a cube root for 3-isogenies.

Contributions

The main contribution of this paper is the generalization of the aforementioned special cases of 2- and 3-isogenies to all isogenies of any degree \(N \ge 2\).

Concretely, given an elliptic curve E with a point P of order N, one can use Vélu’s formulae to compute a defining equation for \(E' = E /\langle P \rangle \). We present accompanying formulae which produce a point \(P'\) on \(E'\) again of order N, such that the composition

is a cyclic isogeny of degree \(N^2\). These formulae are algebraic expressions in the coefficients of E and the coordinates of P, and one radical (an Nth root) of another algebraic expression in the coefficients of E and the coordinates of P. An important implication of this construction is that the same formulae now apply to \(E'\) and \(P'\), which allows us to compute chains of N-isogenies of arbitrary length without needing to generate an N-torsion point in every step. In practice, we assume \(P = (0,0)\), thereby suppressing its coordinates from the formulae.

More in detail, we proceed as follows: an elliptic curve E over a field K together with a K-rational point P of order \(N\ge 4\) can be represented by the Tate normal form

We then compute the curve \(E' = E/\langle P \rangle \) using Vélu’s formulae. The point \(P'\) on \(E'\) can be constructed as a pre-image of P under the dual isogeny \(\hat{\varphi } : E' \rightarrow E\), which guarantees that the composition of \(\varphi \) with \(E' \rightarrow E' / \langle P' \rangle \) is cyclic of order \(N^2\). Our central observation is that \(P'\) is defined over \(K(b, c,\! \root N \of {\rho } \, )\) for some \(\rho \in K(b,c)\) and we prove that one can take \(\rho = t_N(P,-P)\) where \(t_N\) denotes the Tate pairing. Indeed, since \(\hat{\varphi }(P') = P\) and using the compatibility of the Tate pairing with isogenies, we have

which shows that the field of definition of \(P'\) must contain \(\root N \of {t_N(P,-P)}\), and we show that this is also sufficient.

The fact that we only require one Nth root explains the name “radical isogenies”. By rewriting \((E', P')\) again in Tate normal form with coefficients \(b'\) and \(c'\), we are ready for another iteration. The formulae we derive in fact express \(b'\) and \(c'\) directly as elements of \(K(b, c,\! \root N \of {\rho } \, )\).

By specializing to finite fields \(\mathbb {F}_q\) with \(\gcd (q-1, N) = 1\), we immediately obtain that the radical \( \root N \of {\rho } \, \) is again defined over \(\mathbb {F}_q\), since Nth powering is a field automorphism in this case. We implemented our formulae and considered two application scenarios: firstly, we show that using our formulae, chains of N-isogenies can be computed much faster than using the state-of-the-art methods: for \(N=3, 5, 7\) the best previous approach was to use modular polynomials and we obtain speed-ups of factors 9, 18 and 27. For \(N = 11, 13\), the best previous approach was to generate N-torsion points in combination with Vélu’s formulae and our radical isogenies outperform this by factors 12 and 5 respectively. Secondly, we implemented a version of CSIDH using radical isogenies for all primes \(\le \)13 and obtain a speedup of \(19\%\) over the state of the art implementation [1].

Paper organization

Section 2 briefly recaps the necessary background on isogenies, division polynomials, the Tate normal form, the Tate pairing, simple radical extensions, and isogeny-based protocols. Section 3 proves the existence of radical isogeny formulae, while Sect. 4 works out these formulae explicitly for small values of N. Section 5 discusses how our formulae perform when computing chains of N-isogenies, while Sect. 6 reports on an improved implementation of CSIDH using radical isogenies. Finally, Sect. 7 concludes the paper and lists a number of open problems.

Throughout this section we let K denote an arbitrary field.

2.1 Isogenies and Vélu’s Formulae

Let E and \(E'\) be elliptic curves over K. An isogeny \(\varphi :E\rightarrow E'\) is a non-constant morphism such that \(\varphi (\mathcal {O}_E)=\mathcal {O}_{E'}\), where \(\mathcal {O}_E, \mathcal {O}_{E'}\) denote the respective points at infinity. The degree of \(\varphi \) is its degree as a morphism and there always exists a dual isogeny \(\hat{\varphi }:E' \rightarrow E\) such that \(\hat{\varphi } \circ \varphi = [\deg (\varphi )]\), where as usual \([ \cdot ]\) denotes scalar multiplication. The kernel of \(\varphi \) is a finite subgroup of E, more precisely its size is a divisor of \(\deg (\varphi )\), where equality holds if and only if \(\varphi \) is separable (which is automatic if \({{\,\mathrm{char}\,}}K \not \mid \deg (\varphi )\)). Conversely, given a finite subgroup \(C \subset E\), there exists a unique^{Footnote 1} separable isogeny \(\varphi \) having C as its kernel. Concrete formulae for this isogeny were given by Vélu:

Theorem 1

Let C be a finite subgroup of the elliptic curve

over K. Fix a partition \(C = \{\mathcal {O}_E\}\cup C_2 \cup C^+ \cup C^-\), where \(C_2\) are the order 2 points of C, and \(C^+\) and \(C^-\) are such that for any \(P\in C^+\) it holds that \(-P\in C^-\). Write \(S = C^+\cup C_2\), and for \(Q\in S\) define

Then the separable isogeny \(\varphi \) with domain E and kernel C has codomain \(E' = E/C\) with Weierstrass equation

over \(\overline{K}\). Furthermore, for \(P\in E\) we can compute the image of P as

Proof

See [28].   \(\blacksquare \)

2.2 Division Polynomials

Let E/K be defined by \(y^2 + a_1 xy + a_3 y = x^3 + a_2 x^2 + a_4 x + a_6\), and let \(b_2 = a_1^2 + 4a_2\), \(b_4 = 2a_4 + a_1 a_3\), \(b_6 = a_3^2 + 4a_6\), \(b_8 = a_1^2 a_6 + 4a_2a_6 - a_1a_3a_4 + a_2a_3^2 - a_4^2\). For all integers \(N \ge 0\), the N-division polynomial is given by

where \(t = N\) if N is odd and \(t = \frac{N}{2} \cdot \varPsi _{E,2}\) if N is even. By definition, we have that for any non-trivial \(P \in E[N]\), \(\varPsi _{E,N}(P) = 0\). The division polynomials satisfy the following recurrence relation which allows them to be computed efficiently:

Note that \(\varPsi _{E,2}^2 = 4x^3 + (a_1^2 + 4a_2)x^2 + (2a_1a_3 + 4a_4)x + a_3^2 + 4 a_6\), i.e. a univariate polynomial in x.

If one is interested in points of exact order N (so not just in E[N]), then one can use the reduced N-division polynomial \(\psi _{E,N}\) defined as

For all primes \(\ell \), we have that \(\varPsi _{E,\ell }=\psi _{E,\ell }\). Note that for \(N > 2\), the reduced N-division polynomial of an elliptic curve E is a univariate polynomial in x.

The multiplication by N-map can be expressed explicitly using division polynomials as follows [23, Exercise 3.6]:

with \(\phi _{E,N} = x \varPsi _{E,N}^2 - \varPsi _{E,N+1} \varPsi _{E,N-1}\) and \(\omega _{E,N} = \frac{1}{2 \varPsi _{E,N}} (\varPsi _{E, 2N} - \varPsi _{E,N} (a_1 \phi _{E,N} + a_3 \varPsi _{E,N}^2))\).

2.3 The Tate Normal Form

We will be interested in elliptic curves E over K with a distinguished point \(P \in E(K)\) of some finite order N. By translating this point to (0, 0) and requiring that the tangent line is horizontal, and with proper scaling, one can easily prove the following lemma; we refer to [25, Lem. 2.1] for further details.

Lemma 2

Let E be an elliptic curve over K and let \(P \in E(K)\) be a point of order \(N \ge 4\), then (E, P) is isomorphic to a unique pair of the form

with \(b,c \in K\) and

The resulting curve-point pair is said to be in Tate normal form.

Given a Tate normal form, the first few scalar multiples of \(P = (0,0)\) are given by simple expressions in b and c, e.g.

Higher multiples can be computed using (3). Using these multiples, for each \(N \ge 4\) one can write down an irreducible polynomial \(F_N(b,c) \in \mathbb {Z}[b,c]\) whose vanishing, along with the non-vanishing of \(\varDelta (b,c)\) and of \(F_m(b,c)\) for \(4 \le m < N\), expresses that P has exact order N. For instance, for \(N = 4\) we find the equation \(F_4(b,c) = c = 0\), by imposing that \(3P = -P\). Similarly, for \(N=5\) we find \(F_5(b,c) = c - b = 0\) and for \(N = 6\) we find \(F_6(b,c) = c^2 + c - b = 0\). Further examples can be found in Table 1 below. Alternatively, the polynomial \(F_N(b,c)\) can be recovered as a factor of the constant term of the N-division polynomial of the curve (4), when considered over the rational function field \(\mathbb {Q}(b,c)\). This is the approach taken in [25, §2], to which we refer for more details.

Remark 3

Up to birational equivalence, \(F_N(b,c)\) is a defining polynomial for the modular curve \(X_1(N)\). See again [25] for more background.

2.4 The Tate Pairing

Given an elliptic curve E/K and an integer \(N \ge 2\), the Tate pairing is a bilinear map

which can be computed as follows. Consider a Miller function \(f_{N,P_1}\), i.e., a function on E with divisor \(N(P_1) - N(\mathcal {O}_E)\). Let D be a K-rational divisor on E that is linearly equivalent with \((P_2) - (\mathcal {O}_E)\) and whose support is disjoint from \(\{ P_1, \mathcal {O}_E \}\). Then \(t_N(P_1,P_2) = f_{N,P_1}(D)\). If \(P_1 \ne P_2\) and the Miller function is normalized, i.e., the leading coefficient of its expansion around \(\mathcal {O}_E\) with respect to the uniformizer x/y equals 1 (we are assuming that E is in Weierstrass form), then one can simply compute \(t_N(P_1,P_2)\) as \(f_{N,P_1}(P_2)\).

For certain instances of K, the Tate pairing is known to be non-degenerate, meaning that for each \(P_1 \in E(K)[N] \setminus \{ \mathcal {O}_E \}\) there exists a \(P_2 \in E(K) / N E(K)\) such that \(t_N(P_1,P_2) \ne 1\), and vice versa. Most notably, this is true if \(K = \mathbb {F}_q\) is a finite field containing a primitive Nth root of unity \(\zeta _N\) [15], i.e., for which \(N \mid q-1\).

Another important feature is that the Tate pairing is compatible with isogenies, in the following sense: if \(\varphi : E \rightarrow E'\) is an isogeny over K then the rule \(t_N(\varphi (P_1), P_2') = t_N(P_1,\hat{\varphi }(P_2'))\) applies. In particular we have

for all \(P_1 \in E(K)[N]\) and \(P_2 \in E(K)/NE(K)\). For a proof of this compatibility we refer to [4, Thm. IX.9], which assumes \(\zeta _N \in K\), but this condition can be discarded (it is not used in the proof).

2.5 Simple Radical Extensions

Following [10], we say that a field extension \(K \subset L\) is simple radical of degree \(N \ge 2\) if there exists an \(\alpha \in L\) such that (i) \(L = K(\alpha )\), (ii) \(\rho := \alpha ^N \in K\), and (iii) \(x^N - \rho \in K[x]\) is irreducible. Property (iii) can be verified easily using the following theorem.

Theorem 4

Let K be a field, consider an integer \(N \ge 2\), and let \(\rho \in K^*\). Assume that for all primes \(m \mid N\) we have \(\rho \notin K^m\). If \(4 \mid N\), assume moreover that \(\rho \notin -4K^4\). Then the polynomial \(x^N - \rho \in K[x]\) is irreducible.

Proof

See [17, Thm. VI.9.1].    \(\blacksquare \)

We will usually write \(L = K(\! \root N \of {\rho } \, )\), although it should be noted that \( \root N \of {\rho } \, \) is only well-defined up to multiplication by \(\zeta _N^i\) for some \(i \in \{0, 1, \ldots , N-1\}\). Apart from this subtlety, we note that the field \(K(\! \root N \of {\rho } \, )\) does not change if we multiply \(\rho \) with the Nth power of an element of \(K^*\), or if we raise \(\rho \) to some power that is coprime with N.

Remark 1 If \(K \subset L\) is simple radical of degree N and if \({{\,\mathrm{char}\,}}K \not \mid N\), then the Galois closure of L over K is obtained by adjoining a primitive Nth root of unity \(\zeta _N\), and

where the first factor is cyclic of order N. In particular, if \(\zeta _N \in L\) then L is Galois over K with cyclic Galois group. Kummer theory provides a converse statement [24, Lem. 9.13.1].

2.6 CSIDH

We briefly review the CSIDH key agreement protocol, which is our main application of radical isogenies. Let \(\mathbb {F}_p\) be a large finite field with \(p = c \ell _1 \ell _2 \cdots \ell _r - 1\), where the \(\ell _i\) are small distinct primes and where c is some small cofactor. Alice and Bob agree on an order \(\mathcal {O}\subset \mathbb {Q}(\sqrt{-p})\) containing \(\mathbb {Z}[\sqrt{-p}]\), and they consider the set \(\mathcal {E}\ell \ell _p(\mathcal {O}) = \mathcal {E}\ell \ell _p(\mathcal {O}, 0)\) of elliptic curves \(E / \mathbb {F}_p\) whose endomorphism ring \(\mathrm {End}_{\mathbb {F}_p} E\) is isomorphic to \(\mathcal {O}\). Such curves are necessarily supersingular, and without loss of generality it can be assumed that the isomorphism \(\mathrm {End}_{\mathbb {F}_p} E \cong \mathcal {O}\) identifies the Frobenius endomorphism \(\pi _p\) on E with \(\sqrt{-p}\).

To any \(E \in \mathcal {E}\ell \ell _p(\mathcal {O})\) and any invertible ideal \(\mathfrak {a}\subset \mathcal {O}\) one can, using the above isomorphism, associate the finite subgroup

It turns out that the isogenous curve \(E / E[\mathfrak {a}]\) is again contained in \(\mathcal {E}\ell \ell _p(\mathcal {O})\) and that it depends on the class \([\mathfrak {a}]\) of \(\mathfrak {a}\) only; furthermore, this defines a free and transitive action of the ideal-class group \({\text {cl}}(\mathcal {O})\) on \(\mathcal {E}\ell \ell _p(\mathcal {O})\). The key agreement then works as follows: Alice and Bob agree on a starting curve \(E \in \mathcal {E}\ell \ell _p(\mathcal {O})\), then both sample a secret ideal-class \([\mathfrak {a}]\) resp. \([\mathfrak {b}]\), compute the isogenous curves \(E / E[\mathfrak {a}]\) resp. \(E / E[\mathfrak {b}]\), and exchange the outcomes. Both parties can now compute \(E / E[\mathfrak {a}\mathfrak {b}]\) by acting with their own secret ideal-class on the other party’s curve.

In order for this to be practical, Alice and Bob should sample \(\mathfrak {a}, \mathfrak {b}\) as products of ideals of the form \((\ell _i, \sqrt{-p} - 1)^{e_i}\), whose action corresponds to a chain of \(|e_i|\) easy-to-compute \(\ell _i\)-isogenies; this is also true if \(e_i < 0\), in which case one considers the equivalent ideal \((\ell _i, \sqrt{-p} + 1)^{|e_i|}\). The prime \(\ell _i = 2\) requires special treatment: it should be skipped unless \(p \equiv 7 \bmod 8\) and \(\mathcal {O}\) is the maximal order, in which case one considers \((2, (\sqrt{-p} - 1)/2)\) resp. \((2, (\sqrt{-p} + 1)/2)\) instead of the principal ideals \((2, \sqrt{-p} - 1), (2, \sqrt{-p} + 1)\).

In this section we prove the existence of radical isogeny formulae, without deriving these formulae explicitly. The explicit derivation for small N, including the cases \(N = 2,3\), is given in the next section. As such, we assume \(N \ge 4\) and consider the ‘universal’ Tate normal curve

over the field

so that the base point \(P = (0,0)\) has order N. Note that \(\mathbb {Q}_N(b,c)\) is simply the function field of \(X_1(N)\) over \(\mathbb {Q}\). Let \(\varphi : E \rightarrow E'\) be the isogeny with kernel \(\langle P \rangle \); for concreteness it can be assumed that the codomain curve \(E'\) is given by Eq. (2) provided by Vélu’s formulae, although this is not needed for what follows.

Recall that we are interested in those points \(P' \in E'\) for which the composition

is a cyclic \(N^2\)-isogeny. It is easy to check that these points are characterized by the condition

(5)

with \(\hat{\varphi } : E' \rightarrow E\) the dual of \(\varphi \). In particular, there are \(N \phi (N)\) such points, generating N distinct subgroups of \(E'\), where \(\phi \) denotes Euler’s totient function. The points corresponding to \(\lambda = 1\) will be called P-distinguished; they can be viewed as a set of canonical generators for these subgroups.

Define

where the Miller function \(f_{N,P}\) on E is assumed to be normalized, so that \(\rho \) is just \(t_N(P, -P)\) when considered modulo Nth powers in \(\mathbb {Q}_N(b,c)^*\). The main result of this section is:

Theorem 5

Let \(P' \in E'\) be a point satisfying (5). Then the field extension \(\mathbb {Q}_N(b,c) \subset \mathbb {Q}_N(b,c)(P')\), obtained by adjoining the coordinates of \(P'\), is simple radical of degree N. More precisely, \(\mathbb {Q}_N(b,c)(P') = \mathbb {Q}_N(b,c)(\! \root N \of {\rho } \, )\) for an appropriately chosen Nth root \( \root N \of {\rho } \, \) of \(\rho = f_{N,P}(-P)\).

Proof

The fibre \(\hat{\varphi }^{-1} \{ \lambda P \}\) decomposes as a union of orbits under the action of the absolute Galois group of \(\mathbb {Q}_N(b,c)\), together containing N elements. One of these orbits contains \(P'\). Its number of elements equals the degree of the corresponding closed point, which in turn equals the degree of the extension \(\mathbb {Q}_N(b,c) \subset \mathbb {Q}_N(b,c)(P')\). In particular, this extension has degree at most N. On the other hand, by Lemma 6 below, the extension \(\mathbb {Q}_N(b,c) \subset \mathbb {Q}_N(b,c)(\! \root N \of {\rho } \, )\) is of degree precisely N. Therefore, it suffices to prove that \(\mathbb {Q}_N(b,c)(P')\) contains an Nth root of \(\rho \).

To this end we consider \(\alpha := f_{N, P'}(-P') \in \mathbb {Q}_N(b,c)(P')\), where the Miller function \(f_{N, P'}\) is again assumed normalized, and we let \(\mu \) be such that \(\lambda ^2 \mu \equiv 1 \bmod N\). Modulo Nth powers in \(\mathbb {Q}_N(b,c)(P')^*\) we have

showing that \(\rho \) is indeed the Nth power of some element of \(\mathbb {Q}_N(b,c)(P')\).    \(\blacksquare \)

Lemma 6

The polynomial \(x^N - \rho \in \mathbb {Q}_N(b,c)[x]\) is irreducible.

Proof

According to Theorem 4 it suffices to prove:

1. (i)

   for all primes \(m \mid N\) we have \(\rho \notin \mathbb {Q}_N(b,c)^m\),

2. (ii)

   if \(4 \mid N\) then \(\rho \notin -4 \mathbb {Q}_N(b,c)^4\).

Let \(p \equiv 1 \bmod 2N\) be a prime number such that \(4\sqrt{p} > N^2\). Then the Hasse interval \([p + 1 - 2 \sqrt{p}, p + 1 + 2\sqrt{p} ]\) contains the integers \(\lambda N\) for N consecutive values of \(\lambda \). At least one of these values satisfies \(\gcd (\lambda , N) = 1\). By [27, Thm. 2.4.31] there exists an elliptic curve \(\overline{E} / \mathbb {F}_p\) such that \(\overline{E}(\mathbb {F}_p) \cong \mathbb {Z}/(\lambda N)\), so in particular \(\overline{E}( \mathbb {F}_p)[N^\infty ] \cong \mathbb {Z}/ (N)\). Without loss of generality we can assume that \(\overline{E}\) is in Tate normal form, say with coefficients \(\overline{b}, \overline{c} \in \mathbb {F}_p\), and that \(\overline{P} = (0,0)\) is a point of order N on \(\overline{E}\).

Then, in order to prove (i), assume that \(\rho \in \mathbb {Q}_N(b,c)^m\) for some prime divisor \(m \mid N\). Since Miller functions are compatible with reduction mod p and with specialization at \(\overline{b},\overline{c} \in \mathbb {F}_p\) (this follows, for instance, from Miller’s algorithm), we find that

in turn implying that \(t_N(\overline{Q}, [-N/m] \overline{P}) = 1\) for all \(\overline{Q} \in \overline{E}(\mathbb {F}_p)[N]\). This contradicts the non-degeneracy of the Tate pairing over \(\mathbb {F}_p\) (which contains all Nth roots of unity by our choice of p). Indeed, \([-N/m]\overline{P}\) is a non-trivial element of \(\overline{E}(\mathbb {F}_p) / N \overline{E}(\mathbb {F}_p)\).

As for (ii): if \(4 \mid N\) then \(p \equiv 1 \bmod 8\), from which it follows that \(-1\) and 4 are 4th powers in \(\mathbb {F}_p\), in particular the same holds for \(-4\). As above, if \(\rho \in -4 \mathbb {Q}_N(b,c)^4\) then we can conclude that

again contradicting the non-degeneracy of the Tate pairing.    \(\blacksquare \)

An immediate consequence of Theorem 5 is that for each point \(P' = (x_0', y_0')\) satisfying (5) there exist concrete algebraic formulae

for its coordinates: these are the radical isogeny formulae we are after. Note that, in order to find these formulae explicitly, it suffices to consider the cases where \(P'\) is P-distinguished, i.e., where \(\lambda = 1\). Indeed, all other cases are then dealt with by feeding these formulae to the multiplication-by-\(\lambda \) map from (3). Experimentally, it seems that the P-distinguished case yields the simplest formulae.

Remark 2 Our choice of radicand \(\rho = f_{N,P}(-P)\) is somewhat arbitrary: any representant of \(t_N(P, \mu P)\) for any \(\mu \in (\mathbb {Z}/ N)^*\) would have worked equally well, with the same proofs. This reflects the fact that scaling \(\rho \) by Nth powers, or raising \(\rho \) to an exponent that is coprime with N, results in the same simple radical extension.

Given the coordinates of a P-distinguished point \(P'\), all other P-distinguished points are found by varying the choice of \( \root N \of {\rho } \, \):

Lemma 7

Let \(\lambda \in (\mathbb {Z}/ N)^*\) and consider formulae of the form (7) expressing the coordinates of a point \(P'\) such that \(\hat{\varphi }(P') = \lambda P\). Then, by varying the choice of the Nth root \( \root N \of {\rho } \, \), i.e., by scaling it with \(\zeta _N^i\) for \(i = 0, 1, \ldots , N-1\), these formulae compute the coordinates of all points \(P'\) for which \(\hat{\varphi }(P') = \lambda P\).

Proof

From the proof of Theorem 5 it follows that \(\hat{\varphi }^{-1}\{ \lambda P \}\) consists of a single Galois orbit, which implies our claim.    \(\blacksquare \)

For the applications we have in mind, we want to interpret the formulae (7) in some concrete field K, with the indeterminates b, c replaced by concrete elements \(\overline{b}, \overline{c} \in K\). It follows from general principles in algebraic geometry that these specialized formulae continue to produce the coordinates of a point \(P'\) defining a cyclic \(N^2\)-isogeny, with the possible exception of finitely many field characteristics \(p > 0\) and finitely many \((\overline{b}, \overline{c}) \in K^2\). Loosely based on good reduction arguments from the theory of modular curves, we actually believe:

Conjecture 1

The formulae (7) are compatible with specialization to all fields K satisfying \({{\,\mathrm{char}\,}}K \not \mid N\) and to all elements \(\overline{b}, \overline{c} \in K\) satisfying \(F_N(\overline{b}, \overline{c}) = 0\), \(\varDelta (\overline{b}, \overline{c}) \ne 0\) and \(F_m(\overline{b}, \overline{c}) \ne 0\) for all \(4 \le m < N\) (in other words, to all \(\overline{b}, \overline{c}\) for which \(y^2 + (1 - \overline{c})xy - \overline{b} x = x^3 - \overline{b}x^2\) is an elliptic curve on which \(\overline{P} = (0,0)\) has exact order N).

It is easy to confirm this conjecture for small values of N, by explicitly factoring the N-division polynomial of \(E'\): this is the approach followed in the next section, leading to explicit expressions for the formulae (7). In particular, the above conjecture does not affect any of our conclusions in Sects. 5 and 6, which are based on radical N-isogenies for these small values of N only. But from a purely mathematical point of view, we leave the validity of Conjecture 1 as an interesting open question.

We conclude by recalling that by rewriting \((E', P')\) in Tate normal form, one obtains a curve equation

where now

are certain algebraic expressions in \(b, c, \! \root N \of {\rho } \, \). The formulae (8) can be applied iteratively, effectively allowing to compute a cyclic \(N^k\)-isogeny for arbitrary k without needing to explicitly generate points of order N in each step.

In this section, we explain how to find concrete formulae of the forms (7) and (8) for small values of N, by factoring the reduced N-division polynomial of \(E'\) with the help of Magma [5]. As a by-product, we get a confirmation of Conjecture 1 in these cases. In particular, throughout this section, we work over an arbitrary field K with \({{\,\mathrm{char}\,}}K \not \mid N\).

We first deal with the cases \(N = 2, 3\), which require to use a different curve model. We note however that the same principles, in particular using the Tate pairing, also applies in these cases.

Case \(\varvec{N=2}\). Since \({{\,\mathrm{char}\,}}K \ne 2\), we can assume that \(E : y^2 = x^3 + a_2x^2 + a_4x\) for \(a_2, a_4 \in K\) and \(P = (0,0)\). A simple calculation shows that the isogenous curve \(E / \langle P \rangle \) can be given by

The dual isogeny corresponds to quotienting out (0, 0) on \(E'\), so any other point of order 2 on \(E'\) is a suitable instance of \(P'\); note that it is automatically P-distinguished. If we define \(\rho = a_4\) and \(\alpha = \sqrt{\rho }\), then these points are of the form

and by translating \(P'\) to (0, 0), we find the isomorphic model \(E' : y^2=x^3+a'_2x^2+a'_4x\), where

We are now ready to repeat the whole process, since we can divide out by (0, 0) again.

Remark 3 We cannot use \(f_{2,P}(-P)\) as an instance of \(\rho \) in this case, since \(P = -P\). Nevertheless, the reader can check that \(\rho = a_4\) is a representant of \(t_2(P,-P)\).

Case \(\varvec{N=3}\). By requiring that the inflexion point \(P = (0,0)\) has a horizontal tangent line, we can assume that \(E : y^2 + a_1xy + a_3y = x^3\) for certain \(a_1, a_3 \in K\). Vélu’s formulae yield

as a defining equation for \(E / \langle P \rangle \). The 3-division polynomial of \(E'\) splits as

and one checks through explicit computation that the linear factor is the kernel polynomial of the dual isogeny. Therefore, any root of the cubic factor is the x-coordinate of a P-distinguished point \(P'\). Letting \(\rho = f_{3, P}(-P) = -a_3\) and writing \(\alpha = \root 3 \of {\rho }\), this cubic factor splits as

(note that it splits completely over \(K(\zeta _3)\) in view of Remark 1 and/or Lemma 7). Thus we can take \(x_0' = -a_1\alpha + 3\alpha ^2\) and then one checks that \(y_0' = 4a_3\) is the y-coordinate of the corresponding P-distinguished point \(P' = (x_0', y_0')\). Translating \(P'\) to (0, 0) yields a model

with \(a'_1=-6\alpha +a_1\) and \(a'_3= 3a_1\alpha ^2 - a_1^2\alpha + 9a_3 \), and we can repeat. We recall that the simple radical nature of iterated 3-isogenies is not a new observation, see [12, 14].

Case \(\varvec{N=4}\). For \(N \ge 4\) we switch to the Tate normal form as in Sect. 3. Concretely, for \(N=4\) we have \(F_4(b,c) = c =0\) so we obtain the defining equation \(E: y^2+xy-by=x^3-bx^2\). From Vélu’s formulae we find

as a defining equation for \(E/\langle P \rangle \), with reduced 4-division polynomial

The first linear factor corresponds to the x-coordinate of a generator of the dual isogeny. The second linear factor corresponds to the x-coordinate of a 4-torsion point Q such that 2Q is in the kernel of the dual isogeny. Any root of the quartic factor is the x-coordinate of a P-distinguished point \(P'\). Letting \(\rho = f_{4,P}(-P) = -b\) and writing \(\alpha = \root 4 \of {\rho }\), one can verify that

is such a P-distinguished point. Translating \(P'\) to (0, 0) we find an isomorphic model of \(E'\) given by

with

This formula can be applied iteratively.

Case \(\varvec{N=5}\). For \(N = 5\) we have \(F_5(b,c) = b - c = 0\), so we obtain the defining equation \(E : y^2 + (1 -b)xy -by = x^3 -bx^2\). Vélu’s formulae yield

as a defining equation for the codomain of \(\varphi : E \rightarrow E / \langle P \rangle \). The 5-division polynomial of \(E'\) can be verified to split as

where the quadratic polynomial factor is the kernel polynomial of the dual isogeny. The roots of the first quintic factor are the x-coordinates of the P-distinguished points. Those of the second quintic factor are the x-coordinates of the points \(P'\) for which \(\hat{\varphi }(P') = 2 P\) (i.e., the doubles of the P-distinguished points). Concretely, letting \(\rho = f_{5,P}(-P) = b\) and writing \(\alpha = \root 5 \of {\rho }\), the first quintic factor admits the root

(with all other roots obtained by scaling \(\alpha \) with powers of \(\zeta _5\)) and then one can check that

is the y-coordinate of the corresponding P-distinguished point \(P'\). Translating \(P'\) to (0, 0), we obtain the isomorphic form

where

and again we can repeat.

Case \(\varvec{N=6}\). For \(N = 6\) we have \(F_6(b,c) = c^2 + c - b = 0\), so we work with \(E : y^2 + (1-c)xy - (c^2 + c)y = x^3 - (c^2 + c)x^2\). Vélu’s formulae yield

as a model for \(E' = E / \langle P \rangle \). Its reduced 6-division polynomial \(\psi _{E', 6}(x)\) behaves much like in the degree 4 case: there is a unique interesting factor

whose roots are the x-coordinates of the P-distinguished points \(P' \in E'\). Letting \(\rho = f_{6, P}(-P) = -b^2/c = -c(c+1)^2\) and writing \(\alpha = \root 6 \of {\rho }\), one checks that

is such a root; all other roots are found by scaling \(\alpha \) with some power of \(\zeta _6\). One then verifies that

is the y-coordinate of the corresponding P-distinguished point \(P'\). When writing \((E',P')\) in Tate normal form, we find

with

Once again, this formula can be applied iteratively.

Radical Isogenies of Degree \(\varvec{N \ge 7}\). A similar reasoning can be made for \(N \ge 7\), but a direct factorization of the reduced N-division polynomial of \(E'\) over \(\mathbb {Q}_N(b,c)(\! \root N \of {\rho } \, )\) quickly becomes unwieldy, for several reasons: the coefficients of \(E'\) become more involved, the degree of \(\psi _{E',N}\) grows quadratically, and both \(\rho \) and the base field \(\mathbb {Q}_N(b,c)\) become increasingly complicated, see Table 1. For instance, from \(N = 7\) onwards it is no longer possible to eliminate one of the variables b, c using the relation \(F_N(b,c) = 0\). As long as the modular curve \(X_1(N)\) has genus 0, it is possible to get around this by using a different parametrization, see Table 2, but for \(N=11\) and \(N \ge 13\) this is no longer the case.

An approach that already works much better is to use number fields, i.e. assign a large enough integer value to b, construct the number field defined by \(F_N(b,c)=0\) and the degree N extension by adjoining \( \root N \of {\rho } \, \). The root of \(\psi _{E',N}(x)\) is an expression in c and \( \root N \of {\rho } \, \) with rational coefficients. We know that each such coefficient is a rational function in b, so if b is large enough, this function can be found using lattice reduction. The most effective method is similar to the previous method, but uses p-adic fields instead of number fields. Again we need to choose a “large enough” value for b and a large enough precision with which we represent the p-adic field, to be able to reconstruct the rational function in b. We followed this approach for \(N = 13\), since Magma struggles to find the formulae using direct root finding. All formulae for \(N=2,\ldots ,13\) can be found online at https://github.com/KULeuven-COSIC/Radical-Isogenies.

In this section we use our iterable radical isogeny formulae of the form (8) to compute chains of N-isogenies between elliptic curves over finite fields \(\mathbb {F}_q\) with \({{\,\mathrm{char}\,}}{\mathbb {F}_q} \not \mid N\); the application to CSIDH is given in Sect. 6. Here we just concentrate on the computation of long chains of N-isogenies for some fixed \(N \ge 2\), and address the following two issues. Firstly, the radicand \(\rho \) might not admit an Nth root over \(\mathbb {F}_q\): in the worst case, this could mean that at every iteration we need to replace the base field with a degree N extension. Secondly, over \(\overline{\mathbb {F}}_q\) there are N choices for \( \root N \of {\rho } \, \), hence the question arises which root to take if we want to navigate the N-isogeny graph in a controlled way. We discuss three special cases given by \(\gcd (q-1,N) = 1\), \(\gcd (q-1,N) = N\) and \(\gcd (q-1,N) = 2\).

5.1 The Case \(\gcd (q - 1, N) = 1\)

The most straightforward case is \(\gcd (q-1, N)=1\), where there is a very natural choice for \( \root N \of {\rho } \, \). Indeed, in this case the map \(\mathbb {F}_q \rightarrow \mathbb {F}_q : a \mapsto a^N\) is a bijection, so if the starting curve \(E : y^2 + (1-c)xy -by = x^3 - bx^2\) is defined over \(\mathbb {F}_q\), then so is \(\rho (b,c)\) and it admits a unique Nth root which is again defined over \(\mathbb {F}_q\). Choosing this instance of \( \root N \of {\rho } \, \) results in new coefficients \(b', c' \in \mathbb {F}_q\) and the argument repeats. Moreover, the Nth root can be computed as \(\rho ^\mu \) where \(\mu \) is such that \(\mu N \equiv 1 \bmod (q-1)\). Thus, the condition \(\gcd (q-1, N)\) naturally pulls out a chain of N-isogenies whose cost, at least for small N, is dominated by a single \(\mathbb {F}_q\)-exponentiation at each step.

Lemma 8

Assume that \({{\,\mathrm{char}\,}}{\mathbb {F}_q} \not \mid N\) and \(\gcd (q-1,N)=1\), then \(\mathrm {End}_{\mathbb {F}_q}E\) is an imaginary quadratic order which is locally maximal at all primes dividing N, and our chain of N-isogenies corresponds to the repeated action of the ideal class \([(N, \pi _q - 1)]\).

Proof

Observe that

where the last equality follows from \(\gcd (q-1, N) = 1\) along with the fact that \(P = (0,0)\) is an \(\mathbb {F}_q\)-rational point of order N. These properties also imply that

with t the trace of Frobenius, showing that \(\mathrm {End}_{\mathbb {F}_q}E\) is indeed an imaginary quadratic order which is locally maximal at all primes dividing N; see [29, §4]. Thus the isogeny \(E \rightarrow E' = E / \langle P \rangle \) is the horizontal isogeny corresponding to the invertible ideal \((N, \pi _q - 1) \subset \mathrm {End}_{\mathbb {F}_q}E\). Since such isogenies do not change the structure of \(E(\mathbb {F}_q)\), and since choosing the unique \(\mathbb {F}_q\)-rational Nth root of \(\rho \) clearly produces an \(\mathbb {F}_q\)-rational point of order N, the reasoning can be repeated and the lemma follows.    \(\blacksquare \)

Estimating the rough cost of an exponentiation as \(1.5\log q\) multiplications in \(\mathbb {F}_q\), our method should be compared with:

1. (i)

   generating an \(\mathbb {F}_q\)-rational N-torsion point and applying (some form of) Vélu’s formulae; the main cost in this approach is the generation of the N-torsion point, which consists of generating a random point and multiplying by the cofactor \(\#E(\mathbb {F}_q) / N\), taking roughly \(11\log q\) multiplications in \(\mathbb {F}_q\); furthermore this procedure has to be repeated with probability 1/N, which is non-negligible for small N,

2. (ii)

   finding an \(\mathbb {F}_q\)-rational root of \(\Phi _N(x, j(E))\), with \(\Phi _N\) the classical modular polynomial of level N; this roughly amounts to computing \(x^q\) modulo the polynomial \(\Phi _N(x, j(E))\), whose degree is at least \(N+1\), so we estimate this cost as \(1.5(N+1)^2\log q\) multiplications in \(\mathbb {F}_q\).

However, for growing N it becomes unfair to measure the cost of a radical isogeny by merely an exponentiation in \(\mathbb {F}_q\): the algebraic expressions for \(b'\) and \(c'\) in terms of \(b, c, \! \root N \of {\rho } \, \) become increasingly complicated, and the cost of evaluating these expressions quickly overtakes the cost of the exponentiation as shown in Table 3. We also remark that the majority of the multiplications are with small constants coming from the explicit formulae as illustrated in Sect. 4. The size of these constants also grows with N, e.g. for \(N=13\) the constants have a size of up to 14 bits.

A similar overhead is present in approach (ii) using modular polynomials (where moreover one is left with the task of determining the correct twist), which seems consistently outperformed by our radical isogeny formulae. As for the basic approach (i) using Vélu’s formulae, it is shown in Table 4 that for small N, radical isogenies are up to 50 times faster, the main reason being that radical isogenies can be chained without explicitly generating a new N-torsion point on each curve. From \(N \approx 15\) onwards, the overhead becomes so large that radical isogenies become less efficient.

5.2 The Case \(\gcd (q-1, N) = N\)

At the other extreme, if \(N \mid q-1\) then \(\mathbb {F}_q\) contains a primitive Nth root of unity \(\zeta _N\). As a consequence, if \(\rho \in \mathbb {F}_q^*\) admits an Nth root \( \root N \of {\rho } \, \in \mathbb {F}_q\), then all Nth roots are defined over \(\mathbb {F}_q\). But the probability that a random \(\rho \in \mathbb {F}_q^*\) admits an Nth root in \(\mathbb {F}_q\) is 1/N only, so one would expect that the base field needs to be extended at most steps of the iteration.

The situation is much better in the following special case: let \(q = p^2\) for some prime \(p \equiv -1 \bmod N\), so that indeed \(N \mid q-1\), and let \(E / \mathbb {F}_{q}\) be a supersingular elliptic curve, say with \(|E(\mathbb {F}_{q}) |= (p+1)^2\). Such curves are used in the CGL hash function and in SIDH, but since these rely exclusively on 2 and 3 isogenies which are already heavily optimized, we do not expect any real improvement for these applications. On these curves we have \(\pi _q = [-p]\), from which it follows that \(E[N] \subset E(\mathbb {F}_q)\). Let \(P \in E\) be any point of order N, then we claim that \(\rho = f_{N,P}(-P) \in \mathbb {F}_q^*\) is an Nth power, i.e. \(t_N(P, -P) = 1\).

To see this, note that the codomain of \(\varphi : E \rightarrow E' = E / \langle P \rangle \) again satisfies \(|E'(\mathbb {F}_q) |= (p+1)^2\) and therefore \(E'[N] \subset E'(\mathbb {F}_q)\). In particular, any P-distinguished point \(P'\) takes coordinates in \(\mathbb {F}_q\) and we conclude

The argument of course repeats, so in this case one can keep applying our radical isogeny formulae, choosing an Nth root of \(\rho \) at each iteration, without ever leaving \(\mathbb {F}_q\). A performance comparison with the modular polynomial method (ii) from the previous section can be found in Table 5.

5.3 The Case \(\gcd (q-1, N ) = 2\)

An interesting intermediate case is \(\gcd (q-1, N) = 2\), where an element \(\rho \in \mathbb {F}_q^*\) is an Nth power if and only if it is a square. If it is, then it has exactly two Nth roots \(\pm \! \root N \of {\rho } \, \). If \(q \equiv 3 \bmod 4\) then one of these Nth roots is a square and one of them is not; they can be computed as \(\rho ^\mu \) resp. \(-\rho ^\mu \), where \(\mu \) is such that \(\mu N \equiv 1 \bmod (q-1)/2\).

For \(N = 2\), it was observed in [6] that this distinction allows for a controlled navigation of the 2-isogeny graph of supersingular elliptic curves E over a finite prime field \(\mathbb {F}_p\) with \(p \equiv 7 \bmod 8\). Concretely, such curves come in two types: curves on ‘the floor’ have endomorphism ring \(\mathbb {Z}[\sqrt{-p}]\) and admit a unique \(\mathbb {F}_p\)-rational point of order 2, while curves on ‘the surface’ have endomorphism ring \(\mathbb {Z}[(1+\sqrt{-p})/2]\) and have three distinguished \(\mathbb {F}_p\)-rational points of order 2:

• \(P^-\), whose halves have x-coordinates that are not defined over \(\mathbb {F}_p\),

• \(P_1^+\), whose halves are not defined over \(\mathbb {F}_p\), but their x-coordinates are,

• \(P_2^+\), whose halves are defined over \(\mathbb {F}_p\)

(see Fig. 1). Quotienting out \(P^-\) takes us from the surface to the floor, while quotienting out \(P^+_1\) and \(P^+_2\) amounts to traveling along the surface, using the horizontal isogenies corresponding to the respective ideals \((2, (\sqrt{-p} + 1)/2)\), \((2, (\sqrt{-p} - 1)/2)\) of \(\mathbb {Z}[(1 + \sqrt{-p})/2]\), see [6, Lem. 5].

Lemma 9

If the curve point pair \((E, P_1^+)\) resp. \((E, P_2^+)\) is in the form (E, P) with

as in Sect. 4, then \(\rho = a_4\) is a square. Applying the iterative formulae (9) corresponds to the repeated action of \([(2, (\sqrt{-p} + 1)/2)]\) resp. \([(2, (\sqrt{-p} - 1)/2)]\) if one consistently computes \(\sqrt{\rho }\) as \(-\rho ^\mu \) resp. \(\rho ^\mu \).

Proof

The fact that \(\rho = a_4\) is a square follows from the proof of [6, Lem. 3]. From [6, Lem. 4] it follows that selecting \(-\rho ^\mu \) resp. \(\rho ^\mu \) corresponds to selecting \(P_1'^+\) resp. \(P_2'^+\) on \(E'\), which implies the lemma. Note that the other square root of \(\rho \) corresponds to \(P'^-\) in both cases, taking us to the floor.    \(\blacksquare \)

A connected component of the 2-isogeny graph of supersingular elliptic curves over \(\mathbb {F}_p\) with \(p \equiv 7 \bmod 8\), highlighting two elliptic curves on the surface together with their three distinguished 2-torsion points and the corresponding 2-isogenies.

Full size image

The first observation, namely that \(\rho \) is a square, generalizes to all N satisfying \(\gcd (p-1, N) = 2\), where we continue to work over \(\mathbb {F}_p\) with \(p \equiv 7 \bmod 8\). More precisely, consider a curve E on the surface, let us say in Tate normal form with \(P = (0,0)\) a point of order \(N \ge 4\). The cyclic N-isogeny \(\varphi : E \rightarrow E' = E / \langle P \rangle \) is the composition of a horizontal N/2-isogeny, i.e. to another curve on the surface, and either (i) a horizontal 2-isogeny or (ii) a vertical 2-isogeny. Then we claim that we are in case (i) if and only if \(\rho \) is a square. To see this, note that we are in case (i) if and only if there exists a point \(P' \in E(\mathbb {F}_p)\) such that the composition of \(\varphi \) with \(E' \rightarrow E' / \langle P' \rangle \) is cyclic of degree \(N^2\). If \(\rho \) is a square then the existence of such a point simply follows from our radical isogeny formulae (7). Conversely, if there exists such a point \(P'\) then we necessarily have \(P = \lambda \hat{\varphi }(P')\) for some \(\lambda \in (\mathbb {Z}/ N)^*\), and it follows from

that \(\rho \) is a square.

Unfortunately, it seems harder to generalize the second observation, but based on experiments we conjecture the following statement for \(N = 4\), in which case we can take \(\mu = (p+1)/8\):

Conjecture 2

Assume that \(N = 4\) and that (E, P) is in Tate normal form

as above. If the isogeny \(E \rightarrow E / \langle P \rangle \) is horizontal then \(\rho = -b\) is a square. Moreover, applying the iterative formula (10) corresponds to the repeated action of \([(2, (\sqrt{-p} - 1)/2)]^2\) if one consistently computes \(\alpha = \root 4 \of {\rho }\) as \(-\rho ^\mu \) resp. \(\rho ^\mu \), depending on whether \(p \equiv 7 \bmod 16\) resp. \(p \equiv 15 \bmod 16\).

Note that we have just come to argue why \(\rho = -b\) is indeed a square. Also, since \(P = (0,0) \in E(\mathbb {F}_p)\), we necessarily have that 2P equals \(P_2^+\), the unique point of order 2 whose halves are \(\mathbb {F}_p\)-rational. As a result, since the isogeny \(\varphi : E \rightarrow E' = E / \langle P \rangle \) is cyclic and horizontal, it necessarily corresponds to the action of \([(2, (\sqrt{-p} - 1)/2)]^2\). Therefore, the main open problem in proving the conjecture is the last claim. So far, we did not succeed in giving a proof, nor did we manage to generalize its statement to larger values of N.

The core operation in CSIDH [7] is computing a composition of many horizontal isogenies of small prime degree \(\ell _i\) for \(i = 1, \ldots , n\), where the \(\ell _i\) are typically consecutive small primes starting from 2 or 3. The exact composition that needs to be computed can be specified as an exponent vector \([e_1, \ldots , e_n]\), where each \(e_i \in [-B_i, B_i]\) indicates how many horizontal isogenies of degree \(\ell _i\) have to be computed. In practice often \(B_i = B\) for all i, where B is some fixed small value such that \((2B+1)^n > 2^{2\lambda }\), with \(\lambda \) the (classical) security parameter. In the previous section, we considered this problem for a single \(\ell _i\) and showed that generating a new \(\ell _i\)-torsion point in every step is expensive.

In CSIDH this problem is (partly) remedied by chaining isogenies of distinct degrees, i.e. computing a horizontal isogeny of degree \(N = \prod _{i=1}^{n} \ell _i^{\delta _i}\) where \(\delta _i = 1\) if \(|e_i| > 1\) and zero otherwise. Without loss of generality we will assume that all \(\delta _i = 1\). Instead of generating an \(\ell _i\)-torsion point in every step, one first generates a point Q of order (possibly dividing) N and then pushes Q through the isogeny chain. Denote with \(Q_k = \varphi _k(Q)\) with \(\varphi _k\) the isogeny of degree \(N_k = \prod _{i=1}^{k} \ell _i\), then if Q had order N at the start, \(Q_k\) will have order \(M_k = N/N_k\). To generate a point of order \(\ell _{k+1}\) it therefore suffices to compute \([M_k/\ell _{k+1}]Q_k\), which is much cheaper than a full scalar multiplication, certainly for larger k. Note that in practice the original point P does not necessarily have order N, so this procedure might skip a few \(\ell _i\). This method therefore amortizes the cost of one full scalar multiplication (to generate the initial Q) over the different primes \(\ell _i\), and only requires a multiplication by \([M_k / \ell _{k+1}]\) in step k. Table 4 shows that pushing a point through an isogeny is a rather cheap operation, and the main costs are still the generation of the initial Q’s and the scalar multiplications by \([M_k / \ell _{k+1}]\). Table 4 also shows that excluding the computation of an N-torsion point, computing a radical isogeny of degree N is slower than a simple application of Vélu’s formulae.

For the above approach, it is clear that the number of initial Q’s that need to be generated is (at least) \(\max _i B_i\), so it typically does not make sense to sample the exponent vectors from a very skew box, i.e. to take \(B_1 \gg B_n\), even though computing an isogeny of degree \(\ell _1\) is much cheaper than computing an isogeny of degree \(\ell _n\). However, using radical isogenies it does make sense to really skew the box since for every prime \(\ell _i\) one only needs to generate one Q.

Implementation

To illustrate this approach, we implemented a variant of CSIDH that also uses radical isogenies to compute the class group action. Our implementation uses Magma v.2.25-2 [5] and is available at https://github.com/KULeuven-COSIC/Radical-Isogenies and builds upon the code from [1]. Concretely, for 128 bits of classical security, consider the field \(\mathbb {F}_p\), with p the CSURF-512 prime from [6], i.e.

In the implementation of [1], the authors used \(B_i = 5\) for all i, however using radical isogenies we propose the skew box

These vectors represent the action of classes of ideals of the form

on elements from the set of public keys \(S_p^-= \{A\in \mathbb {F}_p \mid y^2=x^3+Ax^2-x \text { is a supersingular elliptic curve} \}\). The set \(S_p^-\) is in 1-to-1 correspondence with \(\mathbb {F}_p\)-isomorphism classes of supersingular elliptic curves, which allows for a slightly easier key validation than using Montgomery curves. The set I contains approximately \(2^{256}\) integer vectors, and just as in [7], we heuristically assume that these vectors represent the elements in the class group quasi-uniformly.

The first step in computing the class group action is finding a 4-torsion point P, such that if we compute the isogeny \(\varphi :E\rightarrow E/\langle 2P\rangle \), it holds that \(\varphi (P)\) has halves defined over \(\mathbb {F}_p\). In accordance with Conjecture 2 and the discussion following it, this implies that the isogeny with kernel \(\langle P\rangle \) will then correspond to the action of \(\left( 2,\frac{\sqrt{-p}-1}{2} \right) ^2\). In order to iteratively compute this horizontal \(4^{\lfloor e_1/2\rfloor }\)-isogeny, we first swap to the Tate normal form by translating P to (0, 0). After iterating the 4-isogeny formula \(\lfloor e_1/2\rfloor \) times, we perform a vertical isogeny to a Montgomery representation of an elliptic curve on the floor. If \(e_1\) is odd, we do a single horizontal 2-isogeny on the Montgomery curve, as explained in [1].

The rest of the computation is done on Montgomery curves on the floor for two reasons. The first is that arithmetic on Montgomery curves is slightly more efficient than arithmetic on curves represented by elements of \(S_p^-\). The second reason is that, in order to compute 3-, 5-, 7-, 11- and 13-isogenies, we will need to swap between elliptic curves in Tate normal form and Montgomery curves. Computing the Montgomery representation of an elliptic curve is essentially finding a two-torsion point, which in practice means finding a solution to a cubic equation. If a cubic equation has three solutions, the explicit formulae to compute any single one of them require going through a quadratic field extension, even if all solutions are defined over the ground field.^{Footnote 2} An elliptic curve on the floor however, only has one nontrivial two-torsion point. In this case, the cubic equation has exactly one solution over \(\mathbb {F}_p\), and the formula to find it does not require field extensions.

We then compute a horizontal \(3^{e_2}\)-isogeny as follows. We first sample a 9-torsion point and swap to the Tate normal form by translating this point to (0, 0). Next, we calculate a \(9^{\lfloor e_2/2\rfloor -1}\)-isogeny iteratively. We perform one last 9-isogeny using Vélu’s formulae on the Tate normal form with kernel generator (0, 0), before swapping back to the Montgomery form of this curve. The reason for this choice is that one more iteration of the formulae would be more expensive, since we already know the final 9-torsion point and hence can simply use Vélu’s formulae. If \(e_2\) is odd, we will compute this final 3-isogeny together with the \(\ell \)-isogenies for \(\ell \ge 17\).

The \(\ell ^{e_i}\)-isogenies for \(\ell ^{e_i}=5^{e_3},7^{e_4},11^{e_5},13^{e_6}\) are then iteratively computed in a similar manner. We first compute an \(\ell \)-torsion point on a Montgomery curve to swap to the Tate normal form. Next, we iterate the formulae for \(\ell \)-isogenies \(e_i-1\) times, and the final \(\ell \)-isogeny is computed using Vélu’s formulae, at which point we go back to the Montgomery representation of the curve. The only noteworthy exception is that if \(|e_i|=1\), we use the original computation of the CSIDH class group action. The reason for this is that swapping to a Tate normal form requires sampling an \(\ell \)-torsion point, which means it is more efficient to perform this action together with the \(\ell \)-isogenies for \(\ell \ge 17\).

The rest of the \(\ell \)-isogenies for \(\ell \ge 17\) are performed as in [7], where optimizations such as those of [1] can be applied. At the end, we perform one final vertical isogeny to the surface to obtain a public key in \(S_p^-\).

We set the bound to swap to the new formulae of [1] at \(\ell > 113\), since this is the threshold where they start outperforming the formulae of [21] in Magma. The box I from which the exponent vectors are sampled was obtained heuristically over a large sample and is near optimal. Over a sample size of 100 000 class group actions each, our variant of CSIDH results in a speed-up of 19% over the one from [1]. We do note that this comparison is with respect to the CSIDH-512 parameter version, since the Magma code from [1] based on the CSURF-512 parameters did not seem to work. Since the CSIDH-512 parameters do not allow horizontal 2-isogenies, a small part of our speed-up can be ascribed to the work of [6].

Starting from a curve E with an N-torsion point P we have proved the existence of explicit formulae for the isogenous curve \(E' = E / \langle P \rangle \) and the coordinates of a point \(P'\) on \(E'\) of order N, such that the composition of \(E \rightarrow E' = E/ \langle P \rangle \) with \(E' \rightarrow E'/ \langle P' \rangle \) is cyclic of degree \(N^2\). This property implies that the formulae can be used repeatedly to compute chains of N-isogenies without generating N-torsion points in each step of the chain. Furthermore, the formulae, which we have described explicitly for \(N \le 13\), only involve basic arithmetic operations, except for the extraction of an Nth root. We have implemented these formulae and used them in two main applications: computing a chain consisting solely of N-isogenies, where we obtained a speed-up ranging from a factor 29 for \(N = 7\) to a factor 5 for \(N = 13\), and an improved implementation of CSIDH which is \(19\%\) faster than the state of the art implementation.

Open Problems. The following problems remain open and are interesting future work:

• Prove Conjecture 1, stating that our formulae have good reduction wherever there is no obvious obstruction.

• Devise a more efficient method for explicitly finding the radical isogeny formulae to avoid our current approach of factoring N-division polynomials as in Sect. 4, which is a major bottleneck.

• Optimize our formulae, e.g. is it indeed true that the P-distinguished case yields the most compact expressions? Using the relations \(\alpha ^N = \rho (b,c)\) and \(F_N(b,c) = 0\), using different instances of \(\rho \), or using different parametrizations of \(X_1(N)\) as in Table 2 or [26], can we rewrite our formulae such that they become more efficient?

• Prove Conjecture 2 on radical isogenies of degree \(N=4\) between supersingular elliptic curves over \(\mathbb {F}_p\) with \(p \equiv 7 \bmod 8\), and generalize it to larger even values of N.

• Measure the impact of our work on constant-time implementations of CSIDH and on the quantum circuits discussed in [3].