Practical Usage of Radical Isogenies for CSIDH

Recently, a radical isogeny was proposed to boost commutative supersingular isogeny Diffie–Hellman (CSIDH) implementation. Radical isogenies reduce the generation of a kernel of a small prime order when implementing CSIDH. However, when the size of the base field increases, field exponentiation, a core component of computing radical isogenies, becomes more computationally intensive. As the size of the field inevitably grows to resist a quantum attack, so it is necessary to discuss the practical utilization of the radical CSIDH. This paper presents an optimized implementation of radical isogenies and analyzes its ideal use in CSIDH-based cryptography with a review of quantum analysis. We tailored the formula for transforming Montgomery curves into the Tate normal form and further optimized the radical 2-isogeny formula and projective versions of the radical 5- and 7-isogenies. Except for CSIDH-512, using only the radical 2-isogeny for all parameters improves performance by 6% to 10%.


I. INTRODUCTION
Isogeny-based cryptography was first proposed by Couveignes in [15]. Couveignes presented a non-interactive key exchange using ordinary elliptic curves defined over F q , whose endomorphism ring is equivalent to a given order O in an imaginary quadratic field. A Diffie-Hellman-like key exchange protocol can be constructed from the commutativity of Cl(O). This work was later rediscovered independently by Rostovtsev and Stolbunov [13], which is now called the CRS scheme. However, the quantum-subexponential attack exists for the scheme [14], and the scheme is inefficient for practical use.
The isogeny-based cryptography regained attention after the introduction of supersingular isogeny Diffie-Hellman (SIDH) by Jao and De Feo [12]. As SIDH uses supersingular elliptic curves, the endomorphism ring is non-commutative, so it resists the attack proposed in [14]. The security of SIDH is based on the difficulty of finding an isogeny between two given isogenous elliptic curves over a finite field, known to be quantum-exponential. The supersingular The associate editor coordinating the review of this manuscript and approving it for publication was S. K. Hafizul Islam .
isogeny key encapsulation (SIKE), a key encapsulation mechanism based on SIDH, was selected as an alternative candidate for NIST PQC standardization round three. However, due to a polynomial-time key recovery attack by Castryck and Decru, SIDH-based cryptosystems are no longer safe [11]. Although various masking methods are presented in [8], [10], and [9], masked variants of SIDH are not yet attractive in terms of performance and key size. Thus, commutative SIDH (CSIDH), described later, could be a more attractive choice.
The CRS scheme was revisited by De Feo, Kieffer, and Smith in [17] and independently by Castryck et al. in [16]. The advantage of the CRS scheme is that it offers efficient and safe public key validation, making it suitable for constructing a noninteractive key exchange [17]. In [17], they modernized the CRS construction by offering a more efficient method to compute the group action and select algorithm parameters. The CRS scheme was further improved by Castryck et al. in [16] by proposing CSIDH, which solves the parameter selection problem of the CRS schemes using supersingular elliptic curves defined over F p . As SIDH-based cryptosystems become inefficient, CSIDH has attracted more researcher interest because various cryptographic primitives can be constructed [24], [25]. The average performance of one group action of CSIDH is around tens of milliseconds, which is faster than other CRS-based protocols.
The advantage of CSIDH-based cryptography is that its key size is smaller than that of any other PQC primitives. However, unlike other PQC primitives, which use simple matrix-vector multiplication as building blocks, isogeny-based cryptography uses complicated elliptic curve arithmetic over a finite field larger than 500 bits. The disadvantage of isogeny-based cryptography is that it is much slower than other PQC primitives. Hence, numerous studies have been proposed to optimize the performance of isogeny-based cryptography. One line of work is to optimize isogeny computation, which can be performed using another form of elliptic curve or by optimizing the isogeny formula. In [27], [28], and [26], hybrid methods employing the birational equivalence between Montgomery and twisted Edwards curves have been proposed for faster implementation. To optimize isogeny computation, Bernstein et al. recently proposed a new method of computing an ℓ-isogeny, reducing the computational cost fromÕ(ℓ) toÕ( √ ℓ) field operations [23]. Another line of work is to tweak the current schemes for faster implementation. In [29], Costello proposed a new type of SIDH called B-SIDH. In this scheme, Alice computes isogenies from a (p + 1)-torsion supersingular curve subgroup, while Bob computes on the (p − 1)-torsion subgroup of the quadratic twist of the curve. In addition, B-SIDH can be viewed as a tweak to SIDH, allowing faster computation on Alice's side with a more reduction-friendly prime field.
For CSIDH, CSURF was proposed in [22], exploiting the horizontal 2-isogenies using the supersingular elliptic curves defined on the surface. Further, CSURF uses supersingular elliptic curves with the endomorphism ring Z[(1 + √ −p)/2] for p ≡ 7 mod 8. They demonstrated that these elliptic curves could be identified with tweaked Montgomery curves (Montgomery − curves), which have elliptic curve arithmetic and isogeny formulae similar to Montgomery curves (Montgomery + curves). Over this prime field, the prime number 2 splits in Q( √ −p), allowing for the use of horizontal 2-isogenies. As a 2-isogeny merely consists of a single exponentiation over F p , adjusting the private key exponent can lead to better performance, and the desired security level can be tailored more precisely. The CSURF method is slower than CSIDH, as the elliptic curve arithmetic and isogeny formula using projective coordinates are slower on Montgomery − curves than on Montgomery + curves.
However, the idea of exploiting the 2-isogeny has extended to the introduction of the radical isogeny in [21]. The CSIDH-based algorithms require isogeny computations of various degrees, and for this operation, a point on an elliptic curve of a specific order must be created to generate a kernel of an isogeny. A random point Q is selected in F p to generate a kernel of a given order, which costs approximately 1.5 log p field multiplications, and is multiplied by some cofactor k, which costs approximately 11 log p field multiplications in CSIDH-based settings. If P = [k]Q equals the identity, another random point is selected to repeat the process. Hence, generating a kernel is a painstaking process, especially for small torsion points where the failure probability is 1/ℓ [18], [21].
Hence, in [21], a novel approach called radical isogeny is introduced that computes chains of n-isogenies. This approach requires sampling at most one n-torsion point. Similar to CSURF, the maximum value of the private key exponent corresponding to primes using radical isogeny can be enlarged, and the maximum value of the private key exponent corresponding to primes not using radical isogeny can be reduced to minimize the number of kernel point generations.

A. OUR CONTRIBUTIONS
This work analyzes the optimal usage of radical isogenies for implementing CSIDH. The following list details the main contributions of this work.
• In this paper, we optimize the radical isogeny formulae in affine and projective versions proposed in earlier studies [21], [30]. We can implement it more efficiently in C by rationalizing the denominator and tailoring the conversion between various curves. In addition, we analyze the radical 3-and 4-isogeny formula in [7] from an implementation perspective. Through these studies, we present the optimized C implementation results of CSIDH with the N -isogeny (N ∈ {2, 3, 5, 7}).
• We review the quantum complexity of CSIDH and derive CSIDH parameters that satisfy NIST security Level 1 according to the power of the quantum adversary. For the first time, we provide the C implementation result of CSIDH with the sliding window method, improving the cost of field exponentiation, a core component of computing radical isogenies. Through several experiments, we conclude that using only the radical 2-isogeny is better with a larger prime field. Except for CSIDH-512, using only the radical 2-isogeny for all parameters improves performance by 6% to 10%. The results of the implementation are presented in Section IV.

B. ORGANIZATION
This paper is organized as follows. Section II introduces the required background. Next, Section III briefly details the radical isogeny and presents the optimization results for degrees of 2, 3, 4, 5, and 7 for implementation. The implementation results are presented in Section IV, and we draw conclusions in Section V.

II. PRELIMINARY
This section introduces two types of Montgomery elliptic curves. Then, CSIDH and the idea of radical isogeny are presented.

A. MONTGOMERY CURVE AND TWEAKED MONTGOMERY CURVE
We let K be a field with characteristics not equal to 2 or 3. The Montgomery curves over K are defined by the following equation: where b(a 2 − 4) ̸ = 0. Throughout the paper, an elliptic curve in the above form is called the Montgomery + curve. When b = 1, we express it as M a . The tweaked Montgomery curves over K are denoted by where b(a 2 + 4) ̸ = 0. Throughout the paper, an elliptic curve in the above form is called the Montgomery − curve. When b = 1, we express it as M − a . It is well known that point arithmetic on M a can be efficiently performed using only the x-coordinates. We let P = (x p , y p ) and Q = (x q , y q ) be points on M a such that x p ̸ = x q , and P − Q = (x p−q , y p−q ). Then, the x coordinates of their sum P+Q, denoted as x p+q , and the doubling of [2]P, denoted as x [2]P , can be computed as follows: We can induce a similar formula for a Montgomery − curve, M − a [22]. We let P = (x p , y p ) and Q = (x q , y q ) be points on M − a such that x p ̸ = x q , and P − Q = (x p−q , y p−q . Then, the x coordinates of their sum P + Q, denoted as x p+q , and the doubling of [2]P, denoted as x [2]P , can be computed as follows: As defined in the above equations, the elliptic curve arithmetic formula on M − a is similar to the case of M a , except for some sign flips in the numerator. However, these sign flips cause changes in the computational costs when using projective coordinates and projective curve coefficients for implementation. In addition, as the isogeny formula is induced using the differential addition formula, the elliptic curve arithmetic and isogeny on M − a are slower than on M a .

B. CSIDH PROTOCOL AND SECURITY 1) CSIDH PROTOCOL
The CSIDH is an isogeny-based Diffie-Hellman-like key exchange protocol proposed by Castryck et al. [16] and uses commutative group action on supersingular elliptic curves defined over a finite field F p . We let O be an imaginary quadratic order and Eℓℓ p (O) denote the set of elliptic curves defined over F p with the endomorphism ring O.
It is well known that the class group Cl(O) acts freely and transitively on Eℓℓ p (O). This group action is represented by where ℓ i values are small, distinct odd primes. We let E be a supersingular elliptic curve over F p such that End p (E) = Z[π], where End p (E) is the endomorphism ring of E over F p and π = √ −p. Note that End p (E) is a commutative subring of the quaternion maximal order End(E). Then, the trace of Frobenius is zero; hence, #E(F p ) = p + 1. As Suppose Alice and Bob want to exchange a secret key. Alice chooses a vector (e 1 , · · · , e n ) ∈ Z n , where e i ∈ [−m, m] for a positive integer m. The vector represents an isogeny associated with the group action by the ideal class

2) QUANTUM SECURITY OF CSIDH
In [5], the quantum security of CSIDH was thoroughly investigated. They revealed that the quantum security of CSIDH depends on the size of the prime field, not on the size of the private key exponent. Hence, to achieve a 128-bit quantum security level, the authors recommended using a prime field of at least 4096 bits. In [5], a 4096-bit prime is presented using 417 small primes. Using all 417 primes for a group action degrades the performance and exceeds the target classical security level.
The meet-in-the-middle type of attack is the best-known classical attack; thus, based on the complexity of this attack, the number of primes to be used varies according to the maximum value of the private key exponent. For example, for a constant-time CSIDH using the method in [19], if the maximum value of the private key exponent is 5, then we can use the 64 smallest primes. If the maximum value of the private key exponent is 1, then we can use the 139 smallest primes. The group action of CSIDH-4096 using the method in [19] takes approximately 23 gigacycles. For details on the quantum analysis, please refer to [5].

C. RADICAL ISOGENIES
Castryck et al. proposed an efficient method to compute small-degree isogenies in [21]. Computing an ℓ-isogeny from an elliptic curve E(F p ) consists of two steps in CSIDH. First, a point P over F p of order ℓ is generated. Second, an isogenous curve E(F p )/⟨P⟩ is generated.
To generate a kernel of a given order, a random point Q is selected in F p , which costs approximately 1.5 log p field multiplications, and is multiplied by the cofactor k = #E(F p )/ℓ, which costs approximately 11 log p field multiplications.
If P = [k]Q equals the identity, another random point is selected to repeat the process. Hence, generating a kernel is a painstaking process, especially for small torsion points where the failure probability is 1/ℓ [18], [21].
Thus, when computing ℓ i -isogenies for 1 ≤ i ≤ n, it is more efficient to sample a n i=1 ℓ i -torsion point and push it through the isogeny to create a chain of isogenies of degrees ℓ 1 , . . . , ℓ n than to generate ℓ i -torsion points for each ℓ i -isogeny. Nevertheless, the probability of failure is higher when creating a small torsion point; therefore, more random points are selected than are needed.
In [21], they proposed an innovative approach to construct a formula to compute chains of n-isogenies for small n. For an elliptic curve E, we let φ : E → E ′ be an n-isogeny, where ker(φ) = ⟨P⟩ for a n-torsion point P on E. The aim is to express the n-torsion point P ′ on E ′ in terms of the coefficients of E and the coordinates of P. Then, by composing an isogeny E ′ → E ′ /⟨P⟩ with φ, we obtain an isogeny of degree n 2 . More explicitly, they applied the fact that an elliptic curve E over a field K with a K -rational point P of order n ≥ 4 can be represented by the Tate normal form: Then, using Vélu's formula, we can compute the isogenous curve E ′ = E/⟨P⟩. The n-torsion point P ′ on E ′ can be expressed in terms of the coefficients of E and coordinates of P through the corresponding dual isogeny E ′ → E. Then, the composition E → E ′ → E ′ /⟨P ′ ⟩ is an isogeny of order n 2 . This method allows for computing chains of n-isogenies of arbitrary length and requires only one n-torsion point for the first step.
In the next section, we specifically state the formula for radical isogeny of degrees 2, 3, 4, 5, and 7, which we use to implement CRADS n . For general formula details, please refer to [21].

III. INTEGRATION AND OPTIMIZATION OF RADICAL ISOGENIES
This section presents the optimization techniques for implementing radical isogenies. To exploit radical isogenies for applications in CSIDH-based algorithms, we chose radical isogenies of degrees 2, 3, 4, 5, and 7 for the following reasons. Other than radical 2-or 4-isogenies, to compute an n e -isogeny, we must have one n-torsion point to start the process. Hence, using m different degrees of radical isogeny requires processing m torsion point generations over a finite field, which is costly. Although this can be minimized using the Elligator method in [20], the advantage of the radical isogeny is that it can minimize the number of randomly generated points with a certain order. However, the radical isogeny formula itself is costly because it requires n-th root computation (exponentiation in this setting). Additionally, as the radical isogeny formula becomes more complicated as the degree increases, we infer that 7 is the upper bound for CSIDH and implementation in C.
In [30], it was noted that using projective curve coefficients for computing radical isogenies is more efficient because it can reduce inversions during the computation of a chain of isogenies. As this applies to radical isogeny of degrees 4, 5, and 7, we apply the optimized version of the following formulas. Moreover, we tailor the transformation between forms of elliptic curves for further optimization. In [7], Onuki and Moriya proposed new representations of the radical isogeny with degrees 3 and 4. Including these formulae, we analyze the radical isogeny formulae comprehensively from the perspective of implementation.
The notation M and E refer to field multiplication and exponentiation, respectively, and we assume 1M ≈ 1S. We consider the field inversion and n-th root computation to be field exponentiation.
Remark 1: Recently, further optimization of the radical isogeny formulae was proposed by Castryck et al. [6]. According to the paper, computing the radical N -isogeny was optimized or newly proposed for N ∈ {2, 3, . . . , 17} ∪ {19}. However, we do not discuss the formulae of [6] because they have not resulted in noticeable improvement, at least in this paper.
A. RADICAL 2 e -ISOGENY 1) RADICAL 2 E -ISOGENY USING THE MONTGOMERY − CURVE For the radical 2-isogeny, we briefly define the formula for supersingular elliptic curves E defined over a finite field F p with p ≡ 7 mod 8, which is the main field used in CSIDH-based algorithms. Over this prime field, curves (E) can be divided into two groups: those located on the floor with the endomorphism ring Z[ √ −p] and a unique F p -rational point of order two or those located on the surface with the endomorphism ring Z[(1 + √ −p)/2] and three distinguished F p -rational points of order two. These three points of order two are categorized as follows: • P − : whose halves have x-coordinates not defined over F p ; • P + 1 : whose halves are not defined over F p , but their xcoordinates are; and • P + 2 : whose halves are defined over F p . As denoted in Lemma 9 in [21], using points P + 1 or P + 2 allows us to compute the chain of 2-isogenies. Additionally, as stated in Proposition 4 of [22], supersingular elliptic curves with the endomorphism ring Z[(1+ √ −p)/2] are F p -isomorphic to the curve M − a . Hence, we optimized the 2 e -isogeny formula in [22]. In [22], an algorithm that computes a chain of 2isogenies is presented by composing the 2-isogeny formula on Montgomery + curves and transformations between a Montgomery − curve and Montgomery + curve.
Step 4 in Algorithm 1 can be rewritten as follows: VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
return sign(e) · a 10: end if Compared to the direct implementation of Step 4, the above equation saves one inversion. The cost for computing 2 eisogeny is 5M+3E+ (e − 1)·(2M+1E), where E refers to a field exponentiation. Over a finite field F p , where p ≡ 3 mod 4, for a ∈ F p , the square root of a is computed as a (p+1)/4 , the inverse of a is computed as a (p−2) , and the square root inverse is computed as a (p+1)(p−2)/4 mod p−1 . As exponentiation dominates the performance of the 2-isogeny, we used the sliding window method.

2) RADICAL 2 E -ISOGENY USING THE TATE NORMAL FORM
When computing the 2 e -isogeny, it is sometimes better to work with a chain of 4-isogenies and compute the 4 e/2isogeny. This approach is applied to implement SIDH-based algorithms, and we describe the corresponding process as stated in [21]. To use a chain of 4-isogenies, we transformed the Montgomery − curve into the Tate normal form: where P = (0, 0) is a 4-torsion point on E b and b ∈ F p . We let r be the x-coordinate of the 4-torsion point on the Montgomery − curve. Then, r is expressed as follows: In addition, b, expressed in terms of r, is where γ = 2r. Applying Vélu's formula results in where E ′ = E/⟨P⟩. To compute consecutive 4-isogenies, we transformed a 4-torsion point P ′ on E ′ to (0, 0). Then, E ′ is isomorphic to the following elliptic curve: If p ≡ 7 mod 16, then α = −b µ . If p ≡ 15 mod 16, then α = b µ , where 4µ ≡ 1 mod (p − 1)/2. After computing a chain of 4-isogenies, we transformed E ′ back to the corresponding Montgomery − curve, The computational cost for transforming a Montgomery − curve to a Tate normal form (i.e., computing b) is 9M+3E.
The cost for computing one 4-isogeny in affine curve coefficients is 3M+2E. Using projective curve coefficients for computing radical isogenies is more efficient [30], as it can reduce inversions during the computation of a chain of isogenies. We let α = A/C for A, C ∈ F p . Then, Now, in the next round of computing the 4-isogeny, we must calculate 4 In projective coordinates, this is equivalent to ( 4 √ X : XZ 4 : Z 2 ), which saves one inversion.

3) RADICAL 2 E -ISOGENY USING THE MONTGOMERY + CURVE
Onuki and Moriya proposed an optimized representation of the radical 4-isogeny in [7]. If M a is a Montgomery curve with coefficient a ∈ F p , and β is a fourth root of 4(a + 2), then the 4-isogeny M a → M a ′ can be computed by To compute the 4-isogeny chain, we calculated the intermediate value corresponding to 4(a ′ + 2) and computed a ′ at the end of the radical 4 e/2 -isogeny. Thus, the cost of computing the radical 4-isogeny once is 3M+2E. The computational cost of transforming the modified Montgomery coefficient of the form 4(a + 2) into the Montgomery coefficient is only 1M + 1A through precomputed constants. If e is an odd number, the last 2-isogeny must be computed at the end of the 4-isogeny chains. This 2-isogeny M a ′ → M a ′′ is as follows: where β is a square root of 4(a ′ + 2). This 2-isogeny can be computed by 1M + 2E. There is an advantage of not transforming the curve form; therefore, we apply this method in the implementation. VOLUME 11, 2023

B. RADICAL 3-ISOGENY 1) RADICAL 3-ISOGENY USING THE WEIERSTRASS CURVE
For a given 3-torsion point Q on M a , we can transform M a into an isomorphic curve of the form: for some a 1 , a 3 ∈ F p , where Q on M a is mapped to a point P = (0, 0) on E. We let r be the x-coordinate of Q. Then, a 1 and a 3 , expressed in terms of r, are as follows: Applying Vélu's formula to E by letting ⟨P⟩ be the kernel results in the 3-isogenous curve E ′ = E/⟨P⟩. A translation that maps the 3-torsion point Q ′ on E ′ to (0, 0) is required to construct a formula for computing the chain of 3-isogenies. Hence, the final curve, obtained by translating Q ′ to (0, 0), is of the form: where a ′ 1 = −6α + a 1 and a ′ 3 = 3a 1 α 2 − a 2 1 α + 9a 3 , for α = 3 √ −a 3 . After computing the chain of 3-isogenies, we transform E ′ back into the corresponding Montgomery curve. The formula for transforming a Weierstrass curve to a Montgomery curve is presented in Magma code in [21]. Like the Weierstrass coefficient a 2 = a 4 = a 6 = 0 in the case of a radical 3-isogeny, we specifically optimized for these circumstances. The computational cost for transforming the Montgomery curve to a (close) Tate curve (i.e., computing a 1 and a 3 ) is 4M+2E. The cost for computing one 3-isogeny is 2M+1E, and the cost for transforming a Weierstrass curve back into the Montgomery form is 16M+4E.

2) RADICAL 3-ISOGENY USING THE MONTGOMERY CURVE
Onuki and Moriya proposed an optimized representation of the radical 3-isogeny in [7]. We let M a be a Montgomery curve with coefficient a ∈ F p and let α be a cube root of t(t 2 − 1), where t is the x-coordinate of the 3-torsion point of M a . Then, the 3-isogeny M a → M a ′ and t ′ , the x-coordinate of the 3-torsion point of M a ′ , can be computed by 3 . The computational cost of a 3-torsion point on the image curve is 4M+1E, slightly more expensive than the Weierstrass version. However, recovering the Montgomery coefficient is much cheaper at 4M+1E. Thus, for the Weierstrass version to be better than this method, the radical 3-isogeny must be performed at least 2.5 k + 8 times, with k ≈ 1.5 log 2 p. Therefore, we implement this radical 3-isogeny.

C. RADICAL 5-ISOGENY
The computation of the radical 5-isogeny follows a process similar to that of the radical 3-isogeny. For a given 5-torsion point Q on M a , we transformed M a into an isomorphic curve of the following form: Applying Vélu's formula to E by letting ⟨P⟩ be the kernel results in a 5-isogenous curve E ′ = E/⟨P⟩. Again, to construct a formula for computing a chain of 5-isogenies, a translation that maps the 5-torsion point Q ′ on E ′ to (0, 0) is necessary. Hence, the final curve, obtained by translating Q ′ to (0, 0), is of the following form: After computing a chain of 5-isogenies, we transformed E ′ back into a corresponding Montgomery curve. The computational cost for transforming a Montgomery curve into a Tate curve (i.e., computing b) is 8M+1E. The cost for computing one 5-isogeny is 5M+2E, and the cost for transforming a Weierstrass curve back into the Montgomery form is 18M+4E. Similar to the 4-isogeny, using the projective curve coefficient as in [30] can save one exponentiation. If α = X /Z , for X , Z ∈ F p , then b ′ can be expressed Because (X ′ : Z ′ ) = (X ′ : Z ′5 ), only one fifth-root computation is required, which can save an inversion.
where N ′ can be expressed in terms of α for α = 7 √ N 5 − N 4 . As N ′ is too large, we do not explicitly state it in this paper. After computing the chain of 7-isogenies, we transformed E ′ back into the corresponding Montgomery curve. The computational cost for transforming a Montgomery curve into a Tate curve (i.e., computing N ) is 10M+1E. The cost for computing one 7-isogeny is 10M+2E, and the cost for transforming a Weierstrass curve back into a Montgomery form is 20M+4E. Table 4 lists the computational cost of a radical isogeny of various degrees when the affine curve coefficient is used. In Table 4, Mont_to_E refers to the transformation from the Montgomery curve to the Tate normal form for odd-degree isogenies and refers to the transformation from the Montgomery + curves to the Montgomery − curves for 2-and 4-isogenies. The ℓ-isogeny refers to the computation of the one ℓ-isogeny, where ℓ ∈ {2, 3, 4, 5, 7}, and E_to_Mont refers to the transformation from the Weierstrass or Tate curve to the Montgomery curve for odd-degree isogenies and refers to the transformation from the Montgomery − curves to the Montgomery + curves for 2-and 4-isogenies.  Table 3 compares the computational cost of radical isogenies in [30] and in this work. In Table 3, the 2-and 3isogenies refer to the computational cost of a radical isogeny using the affine curve coefficients. Hence, affine curve coefficients are used to implement 2-and 3-isogenies, whereas projective curve coefficients are used to implement 4-, 5-, and 7-isogenies.
The computational cost of E_to_Mont combines the transformation from the projective curve coefficients to the affine curve coefficients and the transformation between curves. Last, extra_2_isog refers to the additional computation of the 2-isogeny when the 4-isogeny formula is used to compute the 2 e -isogeny for an odd integer e. Remark 2: Table 3 excludes the multiplication by a small constant as a multiplication count in [30] because multiplication by a constant in radical isogenies can be substituted with addition.
As denoted in Table 4, computing a 4-isogeny once saves only 1M compared to computing a 2-isogeny twice. Moreover, the transformation from the Montgomery − curve to a certain form of an elliptic curve is more costly in the 4-isogeny case than in the 2-isogeny case; therefore, using the radical 4-isogeny in affine coordinates does not reduce computational cost. As denoted in Table 3, computing a 4-isogeny once saves 1E, which can offset the increased computation necessary to change the curve compared to a 2-isogeny.

IV. IMPLEMENTATION RESULTS
This section discusses parameter selection for CSIDH against quantum attacks and presents the implementation results of   CSIDH using various radical isogenies. From this section onward, CSIDH that uses a radical isogeny up to prime degree n is denoted as CRADS n . We measured the performance of CRADS n and CSIDH with various parameters. For CSIDH-k, an approximately k-bit prime was used in the implementation. These primes are in the form p ≡ 7 mod 8 to use the radical 2-isogeny.
To estimate the implementation results, we executed CSIDH and CRADS n with the maximum exponent private key. In addition, to implement large odd-degree isogenies for CSIDH and CRADS n , we used the square root Vélu formula in [23]. In this paper, the square root Vélu formula was applied for isogenies of degrees greater than 100. The strategy for computing a group action in our implementation is summarized in Fig. 1.
Other optimization methods, such as the action strategy, were not applied in the results presented in Tables 9 and 11 to understand the pure influence and availability of radical isogenies. We measured only the group action without validation and averaged over 10 000 rounds. All cycle counts . . , ℓ n are the first n odd prime and ℓ lf refers to last factor of p + 1. k is the number of used odd primes and m is derived from [5] for classical security level 1.
were obtained on a one-core Intel(R) Xeon(R) Gold 6230R CPU at 2.10 GHz, running Ubuntu 22.04 LTS. For the compilation, we used GCC version 11.3.0 with the optimization level -O3.

A. PARAMETER SELECTION FOR CSIDH
As mentioned in Section II-B2, a quantum-subexponential attack occurs against CSIDH. Earlier studies are presented to estimate the quantum complexity of CSIDH [3], [4], [5]. In this paper, we chose parameters with the potential to satisfy NIST security Level 1. With the DW cost, which is the product of the quantum circuit depth and width, we estimated the quantum security of each parameter and compared it with AES-128. We used the c-sieve estimator provided by [5], including the CSIDH oracle cost. We refer to the quantum security of AES algorithms according to MAXDEPTH from [1] and [2]. Table 6 reveals that decreasing the limit of the quantum depth results in an increase in the DW cost of solving AES because the depth limitation does not guarantee enough Grover iterations, resulting in higher costs. Ironically, AES-128 and CSIDH-512/1024 have similar security levels for a stronger quantum adversary. However, in MAXDEPTH 2 40 (a more realistic assumption), CSIDH-3072 is similar to AES-128, indicating that the quantum security level must be comprehensively reviewed by analyzing the development status of quantum computers and various factors in real time.  Thus, we considered p 512 to p 4096 for candidates for quantum security Level 1 and experimented.

B. PERFORMANCE OF CSIDH AND CRADS n
In this implementation, all primes are of the form p = 8( n i=1 ℓ i ) · ℓ lf − 1. We set the basic parameters to satisfy classical security Level 1, as listed in Table 8. In this field, we chose a supersingular Montgomery + curve, M + 0 : y 2 = x 3 +x, as a base curve for CSIDH and a Montgomery − curve, M − 0 : y 2 = x 3 − x, as a base curve for CRADS n . Field exponentiation is the core operation of radical isogeny; thus, we used the sliding window method for effective exponentiation. The window sizes of this implementation are 5 and 7, obtained from the results in Table 5.
As a result of Section III, we used the radical 2-/4-and 3-isogenies of Onuki's method with some adjustments and used the projective radical 5-and 7-isogenies in this paper. The performance of the group action and comparison results of CSIDH and CRADS n are presented in Table 9. The parameter settings are based on Table 8, and radical isogenies are used only for a fixed m.
As indicated in Table 9, even considering that we do not adjust the interval of radical isogenies, it seems inefficient to use odd-degree radical isogenies. Moreover, as the size of the prime field increases, the inefficiency of odd-degree radical isogenies also increases because odd radical isogenies still require the point sampling process and have various operations that require considerable field exponentiation. We executed further experiments with modified intervals to determine the optimal exponents.  Table 11).
In Table 10, (e 2 , e 3 , e 5 , e 7 ) in the RAD row indicate that we ran CSIDH/CRADS n using radical 2 e 2 -, 3 e 3 −, 5 e 5 −, 7 e 7isogenies. Further, (m : k) in the ODD row means that each k odd-degree isogeny is iterated at most m times, respectively. For example, a group action of CRADS 7 with prime p 512 is computed over the following interval: As listed in Table 11, CRADS 2 using the window method leads to a 6% to 10% performance improvement in all prime fields. In p 512 , where the cost of exponentiation is relatively small, CRADS n (n ∈ {3, 5, 7}) outperforms CRADS 2 . However, as the size of the prime field increases, using only a radical 2-isogeny makes CSIDH more efficient.

V. CONCLUSION
There is a risk that a quantum-subexponential attack can change CSIDH parameters at any time. Thus, studying various optimization techniques is valuable. This paper analyzed the optimal use of radical isogenies to implement CSIDH over various prime fields. In this regard, we further optimized the radical isogeny formulae in [21] and [30] and compared them with the results of previous studies.
However, according to the results, odd-degree radical isogenies appear impractical for implementing large CSIDH due to the inefficiency of exponentiation in a large finite field. Nevertheless, the experiments demonstrate that a radical 2-isogeny is still valuable, leading to a 6% to 10% improvement, although other optimizing methods are not considered.
Finally, as the size of the finite field must be increased to resist quantum attacks, more research and optimization are required to use radical isogenies effectively. In particular, a more effective approach to utilizing radical isogeny will be achieved by combining it with previously studied techniques, such as those discussed in [31], to derive the optimal group action strategies. We expect our study to form the basis for extension to more practical use of a radical isogeny in the future.