On the Parallelization of Square-Root Vélu’s Formulas

1. Introduction

Since the foundational work of De Feo, Jao and Plût [1,2], which resulted in the SIKE protocol [3], isogeny-based cryptography has received considerable interest in the development of post-quantum cryptography. Today, the security of SIKE has been broken by a wave of attacks which started with the work of Castryck and Decru [4] and was followed up by the works of Maino, Martindale, Panny, Pope, and Wesolowski [5] and of Robert [6], ultimately demonstrating the existence of a polynomial-time attack on SIKE with any starting curve. Nevertheless, there are still many isogeny-based protocols that remain unbroken, including the CSIDH [7] key exchange protocol and signature schemes like SeaSign [8], CSI-Fish [9] and SQISign [10], the last of which is currently under consideration in the NIST (National Institute of Standards and Technology) process for Standardization of Additional Digital Signature Schemes [11].

While SIKE only required the computation of isogenies of degrees 2 and 3, there has been a tendency in some of the newer isogeny-based protocols to move towards higher and higher prime degrees. This has brought increasing attention to the task of optimizing the evaluation of large-degree isogenies. For instance, the original CSIDH proposal used isogenies of degrees up to 587, and Chávez, Chi, Jacques and Rodríguez have suggested that a level 5 dummy-free implementation would need degrees as high as 2239 [12] after revising the security. More recently, the level 5 version of SQISign that was submitted to the NIST standarization process [13] brings the ceiling up with a colossal degree-318233 isogeny computation. Even the aforementioned attacks on SIKE require the evaluation of large prime-degree isogenies, with the techniques from [5] requiring the evaluation of an isogeny of degree 321193 over a large extension field for a direct key recovery attack. Therefore, optimizing the computation of large prime-degree isogenies is of great interest for constructive purposes as well as for cryptanalysis. Interestingly, a new generation has also risen of isogeny-based protocols which utilize the techniques of the attack constructively, such as FESTA [14], QFESTA [15], IS-CUBE [16] and constructions for VDFs [17] and VRFs [18], of which the last two also involve large-degree isogenies.

Traditionally, the main method for computing isogenies of prime degree ℓ when there exists a point of order ℓ in the curve has been Vélu’s formulas, initially introduced in [19] and expanded upon in [20] and [21]. In view of the ubiquitousness of large-degree isogenies, the more recent improvement from Bernstein, De Feo, Leroux, and Smith [22], which reduced the time complexity from $\mathcal{O}\left(\ell \right)$ to $\tilde{\mathcal{O}}\left(\sqrt{\ell }\right)$, has been a keystone for the viability of many protocols. In addition to the straightforward savings that it provides, this new algorithm also exhibits a substantial degree of parallelizability that has remained largely unexplored.

• Our Contributions. In this work we detail and showcase the parallelizability of the square-root Vélu formulas, focusing on the problem of computing a single large-degree isogeny from a fixed kernel generator and pushing a number of points through said isogeny. We analyze and provide a parallelism-friendly reformulation of the square-root Vélu formulas, and provide a proof-of-concept (PoC) implementation to illustrate the implications of our work. More precisely,

• We propose a new indexing system for the points in the kernel of the isogeny that is similar to the one presented in [22], but which can be naturally split into subsets that are assigned to each core in a multi-core implementation.
• We provide a parallelized expected cost function based on our square-root Vélu variant, as well as a revised expected cost of the sequential Vélu formulas that were presented in [23]. Assuming that the assigned task is to perform an isogeny construction and to push two points through it, the idealized speedup factors for large degrees is up to 1.79, 2.96 and 4.58 faster than the expected sequential cost function using two, four and eight cores, respectively.
• We give a PoC implementation utilising C and OpenMP to compute our isogeny formulas with large odd prime degrees. Our implementation assumes for concreteness that isogenies are evaluated over a prime field such that there exist rational points of the desired degree in the curve, and achieves speedup factors for large degrees of up to $1.74$, $2.54$ and $3.44$ for 2, 4 and 8 cores, respectively, even for medium-sized isogenies.

• Assumptions and generalizations. For the purposes of our implementation, we have searched for a 1792-bit prime such that $p+1$ contains evenly-spaced prime factors ranging from 19 to 321193 (the degree of the isogeny in the attack of [5]). Isogenies are performed over the base prime field, and two points are pushed through each isogeny. These choices were made for concreteness and closely mirror the scenario in CSIDH. On the other hand, they do not adhere exactly to the case of SQISign (which works in a quadratic field extension) nor that of the attack in [5] (which performs isogenies in a much larger field extension). Nevertheless, the techniques that we describe save on a fixed number of field multiplications, which are the dominating part of the total cost. Therefore we expect that similar speedups in percentage would be obtained in those other cases, with the caveat that the savings could even be slightly larger as the parallelization overhead becomes more negligible with respect to the arithmetic cost. As for the theoretical savings, our costs are expressed in terms of the total number of field multiplications and remain valid for any choice of base field.
• Related Work. Different forms of parallelism have been explored at various layers of isogeny-based protocols. For instance, the use of vectorization through Intel’s Advanced Vector Extensions (AVX-512) has been exploited for the finite field arithmetic layer, both in the context of CSIDH [24] and of the now obsolete SIKE protocol [25,26], obtaining speedup factors in the order of $1.5$. Additionally, it has also been proposed to use AVX-512 to batch multiple evaluations of the protocols, leading to increases in throughput of up to $3.6$ for CSIDH [24] and $4.6$ in SIKE [26]. When latency is the focus, vectorization can also be used to push multiple points through an isogeny, leading to the concept of parallel-friendly evaluation strategies which favour pushing points over point multiplications [26,27,28]. These parallel strategies can accelerate the evaluation of large composite-degree isogenies, but do not apply for isogenies of a large prime degree.

This work focuses exclusively on the multi-core optimization for the isolated evaluation of a large prime-degree isogeny, which has received much less attention. To the best of our knowledge, there are only simple parallelization strategies that have been proposed [29] for the linear-complexity Vélu formulas, not exploiting the improved complexity methods from [22]. We also do not consider vectorization: if it was to be used, we expect that the best way to exploit it would be at the field arithmetic layer. This would lead to improvements similar to those from [24,25,26] which are well documented and are largely parallel to our multi-core improvements.

• Outline. The remainder of the manuscript is organized as follows. In Section 2 we give the basic background and description of the original square-root Vélu formulas with their three main building blocks: KPS, xISOG and xEVAL, and derive their expected cost function. In Section 3 we introduce our new framework, explaining the need for a new indexing system in Section 3.1 and then detailing the new algorithms in Section 3.2 and deriving their expected cost function in Section 3.3. Finally, in Section 4, we present our experimental results, comparing the performance to the theoretical savings.

2. Background

We denote by ${\mathbb{F}}_q$ the finite field with $q=p^n$ elements, and p a prime number. We focus on supersingular Montgomery curves, E, given by

$$E:y^2=x^3+A{x}^2+x,A\in {\mathbb{F}}_q\setminus \{\pm 2\}.$$

Remark 1.

In general, a Montgomery curve is defined by the equation $E_{A,B}:B{y}^2=x^3+A{x}^2+x$ such that $B\ne 0$ and $A\ne \pm 2$, but all choices of B are equivalent up to isomorphism. For the sake of simplicity, we write $E_A$ instead of $E_{A,1}$ and omit the constant B.

Viewing the curve as a projective surface, the set $E\left({\mathbb{F}}_q\right)$ of ${\mathbb{F}}_q$-rational points in E forms a group where the point at infinity ∞ acts as the group identity. An order-d point P on E is a point on the curve such that d is the smallest positive integer satisfying

$$\left[d\right]P\sum _{i=1}^{d}P=\infty .$$

We write $E\left[d\right]$ to refer to the d-torsion subgroup $\{P\in E\left(\overline{{\mathbb{F}}_q}\right)\mid \left[d\right]P=\infty \}$ of E.

• Isogenies. An isogeny $\varphi :E\to E^{\prime }$ is a surjective morphism between curves with a finite kernel such that $\varphi \left({\infty }_E\right)={\infty }_{E^{\prime }}$. Such a map is always also a group homomorphism. Two curves E and $E^{\prime }$ are isogenous over ${\mathbb{F}}_q$ if there exists such $\varphi$ connecting them, or equivalently if $\#E\left({\mathbb{F}}_q\right)=\#E^{\prime }\left({\mathbb{F}}_q\right)$. The kernel of $\varphi$ is $\{P\in E\left({\mathbb{F}}_q\right)\mid \varphi \left(P\right)=\infty \}$, denoted by $ker\varphi$. When restricting only to separable maps, the degree of an isogeny is equal to the size of its kernel, and an isogeny is uniquely determined by its kernel up to an isomorphism. An isogeny $\varphi :E\to E^{\prime }$ of degree ℓ is referred to as a ℓ-isogeny, and it has a unique dual (up to isomorphism) $\widehat{\varphi }:E^{\prime }\to E$ of the same degree such that $\varphi \circ \widehat{\varphi }$ and $\widehat{\varphi }\circ \varphi$ are equivalent to the multiplication-by-ℓ map on E and $E^{\prime }$, respectively.

2.1. Computation of ℓ-Isogenies

Let $E_A$ be a Montgomery curve defined over ${\mathbb{F}}_q$ and $P\in E_A\left({\mathbb{F}}_q\right)$ a point of order ℓ. If $\varphi :E_A\to E_{A^{\prime }}$ is an isogeny with ker$\varphi =\langle P\rangle$, then it is possible to find polynomials $g\left(x\right),h\left(x\right)\in {\mathbb{F}}_q\left[x\right]$ such that

$$\varphi (x,y)=\left(\frac{g\left(x\right)}{h\left(x\right)},y{\left(\frac{g\left(x\right)}{h\left(x\right)}\right)}^{\prime }\right).$$

The polynomials $g\left(x\right)$ and $h\left(x\right)$ are related by a formula stated by Elkies [30], so the main task for computing an ℓ-isogeny is that of obtaining $h\left(x\right)$ from P. In the particular case of Montgomery curves, from the formulas of [31] the coefficient of the Montgomery curve $E_{A^{\prime }}$ can be obtained by,

$$A^{\prime }=\frac{2(1+d)}{1-d},d={\left(\frac{A-2}{A+2}\right)}^{\ell }{\left(\frac{h_S\left(1\right)}{h_S(-1)}\right)}^8,$$

(1)

and the x-coordinate ${\varphi }_x\left(x\left(Q\right)\right)$ of the point $\varphi \left(Q\right)$, by

$${\varphi }_x\left(x\left(Q\right)\right)={\left(\frac{x{\left(Q\right)}^{\ell }{h}_S(1/x\left(Q\right))}{h_S\left(x\left(Q\right)\right)}\right)}^2,$$

(2)

where ${\displaystyle h_S\left(X\right)=\prod _{s\in S}(X-x\left(\left[s\right]P\right))}$ and $S=\{1,3,\dots ,\ell -2\}$.

In 2020, Bernstein, De Feo, Leroux and Smith [22] proposed an important improvement in the computation of isogenies. Their key idea exploits the fact that for isogeny evaluations and constructions, it is not required to compute the polynomial h itself, but only its evaluation at a given set of points. We briefly describe their method for evaluating (1) and (2), starting by pointing out that, although the map $Q\mapsto x\left(Q\right)$ is not holomorphic, an algebraic relation exists between the values. More precisely,

Lemma 1

([22] Lemma 4.3). Let $E/{\mathbb{F}}_q$ be an elliptic curve, where q is a prime power. There exist biquadratic polynomials $F_0$, $F_1$, and $F_2$ in ${\mathbb{F}}_q[X_1,X_2]$ such that

$$\left(X-x(P+Q)\right)\left(X-x(P-Q)\right)=X^2+{\displaystyle \frac{F_1(x\left(P\right),x\left(Q\right))}{F_0(x\left(P\right),x\left(Q\right))}}+{\displaystyle \frac{F_2(x\left(P\right),x\left(Q\right))}{F_0(x\left(P\right),x\left(Q\right))}},$$

for all P and Q in E such that $0\notin \{P,Q,P+Q,P-Q\}$.

This allows us to rearrange $h_S\left(X\right)$ in a baby-step giant-step style:

$$h\left(X\right)=\left(\prod _{i\in I}\prod _{j\in J}(X-x\left([i+j]P\right))(X-x\left([i-j]P\right))\right)\left(\prod _{k\in K}(X-x\left(\left[k\right]P\right))\right),$$

(3)

where the index system $\{I,J\}$ is chosen such that $S=(I+J)\cup (I-J)\cup K$, with $K=S\setminus (I+J)\cup (I-J)$ and satisfying the following definition:

Definition 1.

Let S, I and J be finite sets of integers. The $(I,J)$ is an index system of S if:

1. 

the maps $I\times J\to \mathbb{Z}$ defined by $(i,j)⟼i+j$ and $(i,j)⟼i-j$ are both injective and have disjoint images.

2. 

$I+J$ and $I-J$ are both contained in S.

From here, they show that the factor of $h\left(X\right)$ that is related to I and J can be computed efficiently as a resultant. More precisely, letting $I\pm J$ be the union of the sets $I+J$ and $I-J$, the following result is proven.

Lemma 2

([22] Lemma 4.9). Let $E/{\mathbb{F}}_q$ be an elliptic curve, where q is a prime power. Let P be an order-n element of $E\left({\mathbb{F}}_q\right)$. Let $(I,J)$ be an index system such that $I,J,I+J,$ and $I-J$ do not contain any elements of $n\mathbb{Z}$. Then

$$h_{I\pm J}\left(X\right)=\frac{1}{{\Delta }_{I,J}}Re{s}_Z(h_I\left(Z\right),E_J(X,Z))$$

where

$${\displaystyle E_J(X,Z):=\prod _{i\in J}\left(F_0(Z,x\left(\left[j\right]P\right))X^2+F_1(Z,x\left(\left[j\right]P\right))X+F_2(Z,x\left(\left[j\right]P\right))\right)}$$

and ${\Delta }_{I,J}:=Re{s}_Z(h_I\left(Z\right);D_J\left(Z\right))$ where ${\displaystyle D_J\left(Z\right):=\prod _{i\in J}{F}_0(Z,x\left(\left[j\right]P\right))}$.

Remark 2.

Note that ${\Delta }_{I,J}$ in Lemma 2 does not depend on the variable X. The formulas (2) and (1) for computing the codomain curve and point evaluation respectively, include quotients of the polynomial h, so in practice, ${\Delta }_{I,J}$ in Lemma 2 does not need to be computed.

Notice that $h_{I\pm J}$ is precisely the first parenthesis that appears in (3), which can be computed quadratically faster as a resultant, while the second parenthesis (corresponding to $h_K\left(X\right)$) is a leftover product that is still computed in linear time. Based on this, Adj, Chi-Domínguez, and Rodríguez-Henríquez [23] introduced three explicit algorithms for evaluating ℓ-isogenies: KPS (which computes the necessary multiples of the kernel generator), xISOG (which obtains the coefficient of the image curve) and xEVAL (which obtains the image of a point). These algorithms are reproduced here as Algorithm 1, Algorithm 2 and Algorithm 3.

Algorithm 1  KPS

Inputs:  An elliptic curve $E_A/{\mathbb{F}}_q$; $P\in E_A\left({\mathbb{F}}_q\right)$ of order an odd prime ℓ   Output:  $\mathcal{I}$, $\mathcal{J}$, $\mathcal{K}$ such that $(I,J)$ is an index system for S  1: $b\leftarrow \lfloor \sqrt{(\ell -1)}/2\rfloor$  2: $b^{\prime }\leftarrow \lfloor (\ell -1)/4b\rfloor$  3: $J\leftarrow \{2j+1:0\le j<b\}$  4: $I\leftarrow \{2b(2i+1):0\le i<b^{\prime }\}$  5: $K\leftarrow S\setminus (I\pm J)$  6: $\mathcal{J}\leftarrow \left\{x\right(\left[j\right]P):j\in J\}$  7: $\mathcal{I}\leftarrow \left\{x\right(\left[i\right]P):i\in I\}$  8: $\mathcal{K}\leftarrow \left\{x\right(\left[k\right]P):k\in K\}$ v 9: return $\mathcal{I}$, $\mathcal{J}$ and $\mathcal{K}$

Algorithm 2  xISOG

Inputs:  An elliptic curve $E_A/{\mathbb{F}}_q:y^2=x^3+A{x}^2+x$; $P\in E_A\left({\mathbb{F}}_q\right)$ of order an odd prime ℓ; and $\mathcal{I}$,$\mathcal{J}$,$\mathcal{K}$ the output from KPS(P)  Output:  $A^{\prime }\in {\mathbb{F}}_q$ such that $E_A/{\mathbb{F}}_q:y^2=x^3+A^{\prime }{x}^2+x$ is the image curve of a separable isogeny with kernel $\langle P\rangle$  1: $h_I\leftarrow {\prod }_{x_i\in \mathcal{I}}(Z-x_i)\in {\mathbb{F}}_q\left[Z\right]$   2: $E_{0,J}\leftarrow {\prod }_{x_j\in \mathcal{J}}(F_0(Z,x_j)+F_1(Z,x_j)+F_2(Z,x_j))\in {\mathbb{F}}_q\left[Z\right]$   3: $E_{1,J}\leftarrow {\prod }_{x_j\in \mathcal{J}}(F_0(Z,x_j)-F_1(Z,x_j)+F_2(Z,x_j))\in {\mathbb{F}}_q\left[Z\right]$   4: $R_0\leftarrow Re{s}_Z(h_I,E_{0,J})\in {\mathbb{F}}_q$   5: $R_1\leftarrow Re{s}_Z(h_I,E_{1,J})\in {\mathbb{F}}_q$   6: $M_0\leftarrow {\prod }_{x_k\in \mathcal{K}}(1-x_k)\in {\mathbb{F}}_q$   7: $M_1\leftarrow {\prod }_{x_k\in \mathcal{K}}(-1-x_k)\in {\mathbb{F}}_q$   8: $d\leftarrow {\left(\frac{A-2}{A+2}\right)}^{\ell }{\left(\frac{M_0{R}_0}{M_1{R}_1}\right)}^8$ 9: return $\frac{2(1+d)}{1-d}$

Algorithm 3  xEVAL

Inputs:  An elliptic curve $E_A/{\mathbb{F}}_q:y^2=x^3+A{x}^2+x$; $P\in E_A\left({\mathbb{F}}_q\right)$ of order an odd prime ℓ; $\alpha =x\left(Q\right)=[Q_x:Q_z]\ne 0$ of a point $Q\in E_A\left({\mathbb{F}}_q\right)\setminus \langle P\rangle$ and $\mathcal{I}$,$\mathcal{J}$,$\mathcal{K}$ the output from KPS(P)  Output:  The image of $x\left(\varphi \right(Q\left)\right)$ by w, where $\varphi$ is a separable isogeny of kernel $\langle P\rangle$.  1: $h_I\leftarrow {\prod }_{x_i\in \mathcal{I}}(Z-x_i)\in {\mathbb{F}}_q\left[Z\right]$   2: $E_{0,J}\leftarrow {\prod }_{x_j\in \mathcal{J}}(F_0(Z,x_j)/{\alpha }^2+F_1(Z,x_j)/\alpha +F_2(Z,x_j))\in {\mathbb{F}}_q\left[Z\right]$  3: $E_{1,J}\leftarrow {\prod }_{x_j\in \mathcal{J}}(F_0(Z,x_j){\alpha }^2+F_1(Z,x_j)\alpha +F_2(Z,x_j))\in {\mathbb{F}}_q\left[Z\right]$   4: $R_0\leftarrow Re{s}_Z(h_I,E_{0,J})\in {\mathbb{F}}_q$   5: $R_1\leftarrow Re{s}_Z(h_I,E_{1,J})\in {\mathbb{F}}_q$   6: $M_0\leftarrow {\prod }_{x_k\in \mathcal{K}}(1/\alpha -x_k)\in {\mathbb{F}}_q$   7: $M_1\leftarrow {\prod }_{x_k\in \mathcal{K}}(\alpha -x_k)\in {\mathbb{F}}_q$ 8: return ${\left(M_0{R}_0{Q}_x\right)}^2/{\left(M_1{R}_1{Q}_z\right)}^2$

2.2. Computing the Resultants

The main computational task that is required for the square-root Vélu algorithm is the evaluation of the resultants $Re{s}_Z(h_I,E_{i,J})$, which can be achieved in constant-time via a remainder tree approach as described in [23]. This task accounts for approximately 85% of the total computation time, and depends on the size of the sets I, J, and K. For the following discussion, we set $\#J=b$, $\#I=b^{\prime }$ and $\#K=\frac{\ell -1}{2}-2b{b}^{\prime }$. As a remainder, if $f\left(Z\right)\in {\mathbb{F}}_q\left[Z\right]$ splits into linear factors as $f\left(Z\right)=a{\prod }_{0\le i<b}(Z-x_i)$ and $g\left(Z\right)\in {\mathbb{F}}_q\left[Z\right]$, then the resultant is given by

$$Re{s}_Z(f\left(Z\right),g\left(Z\right))=a^b\prod _{0\le i<b}g\left(x_i\right).$$

(4)

For our purposes, the polynomials f and g are provided as a product of $b^{\prime }$ linear polynomials and b quadratic polynomials, respectively. Since the $x_i$ values are readily available, one could compute (4) directly by evaluating g multiple times, but this would lead to a complexity of $\mathcal{O}\left(b{b}^{\prime }\right)=\mathcal{O}\left(\ell \right)$.

Instead, we begin by obtaining the polynomials f and g as a list of coefficients by multiplying together all the terms following a product tree approach, as shown in Figure 1.

After the product tree for f has been built, we compute a corresponding reciprocal tree. Unlike the product tree, the reciprocal tree is computed from the root down. If a node in the product tree contains a polynomial F of degree m, then the corresponding node in the reciprocal tree contains $(\mathfrak{F},c)$, where $\mathfrak{F}$ is a polynomial of degree m and c is a constant such that $re{v}_m\left(F\right)·\mathfrak{F}=cmod{x}^m$. Here, $re{v}_m(·)$ denotes the polynomial with the list of coefficients in reverse order: if ${\displaystyle F=\sum _{i=0}^m{F}_i{x}^i}$, then ${\displaystyle re{v}_m\left(F\right)=\sum _{i=0}^m{F}_{n-i}{x}^i}$. The reciprocal tree is depicted in Figure 2, and is used as an auxiliary tool to construct one last tree called the residue tree.

The residue tree is also built from the root down, and as depicted in Figure 3, each of its nodes contains $cgmodF$, where F is the polynomial in the corresponding node of the product tree of f and c the constant from the reciprocal tree. The leaves of the remainder tree contain a multiple of $g\left(Z\right)mod(Z-x_i)=g\left(x_i\right)$, so a multiple of the resultant is readily obtained by multiplying all the leaves together. The multiple cancels out when taking the ratios at the last step of Algorithm 3 and Algorithm 2. Remarkably, this method allows us to compute (a multiple of) all of the $g\left(x_i\right)$ with a complexity of just $\tilde{\mathcal{O}}(b+b^{\prime })=\tilde{\mathcal{O}}\left(\sqrt{\ell }\right)$. A detailed description for the cost of the product tree, reciprocal tree, and remainder tree is presented in the Appendix A.

2.3. Cost model for the Sequential Square-Root Vélu

An expected cost function for KPS, xISOG, and xEVAL has been presented in [23], but it takes into account some redundant computations and also fails to consider the cost of the remainder tree. We now present a more accurate and optimal cost function.

As in [23], we begin by noting that a better performance is achieved when $\#K$ is as small as possible and $\#I,\#J$ are similar in size. Since we need to have $2(\#I)(\#J)+\#K=(\ell -1)/2$, we set the size of $\#J$ to be $b:=\lceil \sqrt{(\ell -1)}/2\rceil$ and then perform division with remainder to obtain $\#I=\lfloor (\ell -1)/\left(4b\right)\rfloor$ and $\#K\in [0,2b]$. For simplicity, from now on we ignore rounding errors and assume the average value for $\#K$, so that $\#I\approx \#J\approx \#K\approx b\approx \sqrt{(\ell -1)}/2$.

Next, we provide explicit formulas for the number of field multiplications in each of the three procedures as a function of b. This cost model neglects the cost of field additions and subtractions, as well as the savings that a field squaring may have over a regular field multiplication. We consider this a fair compromise, since the total number of squarings is small and the benefit of having a single metric far outweighs the importance of the small error. The cost functions are expressed in terms of the cost of the various tree procedures, which are detailed in the Appendix A.

• Cost function for KPS. Notice that the computation of $h_I$ is performed in both Algorithm 3 and Algorithm 2, so instead we delegate this product tree and the corresponding reciprocal tree to KPS and perform it only once. The total cost of KPS is then reflected by obtaining b multiples of an elliptic curve point (at a cost of $6b$ field multiplications each) for each of $\mathcal{I}$, $\mathcal{J}$ and $\mathcal{K}$, and the cost of a product tree and reciprocal tree for b linear polynomials. We obtain

  $$\begin{array}{cc}\hfill Cos{t}_{\mathtt{KPS}}\left(b\right)& =18b+ProductTre{e}_1\left(b\right)+ReciprocalTree\left(b\right)=4b^{{log}_2\left(3\right)}+2{log}_2\left(b\right)+16b-2.\hfill \end{array}$$

• Cost function for xISOG. The degree-2 factors in $E_{0,J}$ and $E_{1,J}$ can be obtained at the cost of $5b$ field multiplications, and the complete polynomials are then obtained from two product trees of quadratic leaves. We then perform 2 residue trees to compute the two resultants of Algorithm 2, and another $2\ast (b-1)$ multiplications to compute $M_0$ and $M_1$. Algorithm 2 also requires an ℓ-th power exponentiation at an average cost of $1.5log\ell \approx 3({log}_2\left(b\right)+1)$ multiplications, which must be performed separately on the numerator and denominator when working in projective coordinates. Finally, there are another 10 multiplications for the final two steps. We obtain

  $$\begin{array}{cc}\hfill Cos{t}_{\mathtt{xISOG}}\left(b\right)=& 5b+2\times ProductTre{e}_2\left(b\right)+2\times ResidueTree\left(b\right)\hfill \\ \hfill & +2(b-1)+2\times 3({log}_2\left(b\right)+1)+10,\hfill \end{array}$$

  $$\begin{array}{cc}\hfill Cos{t}_{\mathtt{xISOG}}\left(b\right)=& 18b^{{log}_2\left(3\right)}+6{log}_2\left(b\right)-5b+12.\hfill \end{array}$$

• Cost function for xEVAL. In this case, computing the factors for $E_{0,J}$ requires $10b$ field multiplications. However, only one product tree is required since $E_{1,J}$ can be obtained by inverting the coefficient list of $E_{0,J}$, as noted in [23]. We then need two residue trees to compute the resultants in Algorithm 3, $2(2b-1)$ multiplications to compute $M_0$ and $M_1$, and 6 more multiplications for the final step. The total is

  $$\begin{array}{cc}\hfill Cos{t}_{\mathtt{xEVAL}}\left(b\right)=& 10b+ProductTre{e}_2\left(b\right)+2\times ResidueTree\left(b\right)+2(2b-1)+6\hfill \end{array}$$

  $$\begin{array}{cc}\hfill Cos{t}_{\mathtt{xEVAL}}\left(b\right)=& 15b^{{log}_2\left(3\right)}+5b+2.\hfill \end{array}$$

• Total cost function. The total cost function, assuming that two points need to be pushed through the isogeny as is the case in CSIDH and SQISign, is

  $$\begin{array}{cc}\hfill Cos{t}_{\mathtt{V}é\mathtt{lu}}\left(b\right)=& Cos{t}_{\mathtt{KPS}}\left(b\right)+Cos{t}_{\mathtt{xISOG}}\left(b\right)+2\times Cos{t}_{\mathtt{xEVAL}}\left(b\right)\hfill \\ \hfill =& 52b^{{log}_2\left(3\right)}+8{log}_2\left(b\right)+21b+14.\hfill \\ \hfill \end{array}$$

  (5)

3. Parallelizing Square-Root Vélu Formulas

In this section, we present our main results which focuses on the problem of parallelizing Vélu’s formulas for isogeny computations.

An immediate observation is that the two resultants that appear in each of Algorithm 2 and Algorithm 3, which represent the heaviest part of the computation, are completely independent of each other. Therefore, there is a simple way to parallelize the two algorithms in a 2-core implementation by assigning one resultant to each core. We show that these algorithms exhibit a good degree of parallelizability even beyond two cores, by assigning more than one core per resultant. However, doing so requires a new indexing system that allows for a clear separation of the computation between cores, which we introduce in the first subsection.

Additionally, other parts of the algorithms are also apt for parallelization, such as the computation of $M_0$ and $M_1$ where a partial product can be assigned to each core, or the computation of product trees where a subtree can be assigned to each core after a small bit of sequential work. We show also that the computation of point multiples in KPS can be parallelized with our new index system.

After detailing the new index system, provide explicit new algorithms for the square-root Vélu formulas which can be computed by assuming the availability of n-cores, where $n>1$ is a power of two.

3.1. Construction of a New Index System

The main observation that allows us to parallelize the computation of resultants beyond two cores is that they exhibit a multiplicative property. More precisely, if I can be written as the disjoint union

$$I=I_0\cup I_1\cup \dots \cup I_{n/2-1},$$

then

$$Re{s}_Z(h_I,E_{i,J})=Re{s}_Z(h_{I_0},E_{i,J})\times Re{s}_Z(h_{I_1},E_{i,J})\times \dots \times Re{s}_Z(h_{I_{n/2-1}},E_{i,J}),$$

so we can assign one subset $I_t$ to each of the cores, have them compute one subresultant each, and then multiply them all together. Doing so will require us to modify the sizes of I and J: since each resultant computation should have balanced sizes, we now require $\#J\approx \#I/(n/2)$. Accordingly, we now need $\#J\approx \sqrt{(\ell -1)/2n}$ and $\#I\approx \sqrt{(\ell -1)n/4}$. We will design such an indexing system according to the following lemma.

Lemma 3.

Let $n>1$ be a power of two, ℓ be an odd positive integer and consider the set $S=\{1,3,...,\ell \}$. If $b=\lfloor \sqrt{(\ell -1)/2n}\rfloor \ge 1$ and $b^{\prime }=\lfloor (\ell -1)/2nb\rfloor$, then $(I,J)$ is an index system for S, where

$$\begin{array}{cc}\hfill J:=& \{2j+1:0\le j<b\},\hfill \end{array}$$

$$\begin{array}{cc}\hfill I=& {\displaystyle \bigcup _{t=0}^{n/2-1}{I}_t,\mathrm{and}}\hfill \end{array}$$

$$\begin{array}{cc}\hfill I_t:=& \{2b(2t+1)+2bni:0\le i<b^{\prime }\}\mathrm{for}0\le t<n/2.\hfill \end{array}$$

Moreover, $I_t\cap I_{\tilde{t}}=\varnothing$ if and only if $t\ne \tilde{t}.$

Proof of Lemma .

Suppose that

$$(2b(2t+1)+2bn{i}_1))-(2b(2\tilde{t}+1)+2bn{i}_2))=(2j_1+1)-(2j_2+1),$$

where $0\le j_i,j_2<b$, $0\le i_1,i_2<b^{\prime }$ and $j_1\ne j_2$ and $i_1\ne i_2$ (t and $\tilde{t}$ can be equal). Then,

$$4b(t-\tilde{t})+2bn(i_1-i_2)=2(j_1-j_2),$$

so

$$2b(t-\tilde{t}+n(i_1-i_2)/2)=j_1-j_2$$

but this is impossible since $0<|t-\tilde{t}|<n$ and thus $0<|j_1-j_2|<b<2b|(t-\tilde{t}+n(i_1-i_2)/2)|$. Therefore, the map $I\times J\to \mathbb{Z}$ defined by $(i,j)⟼i+j$ is injective. Similarly, the map $(i,j)⟼i-j$ is injective and therefore, both have disjoint images.

On the other hand, the sets $I+J$ and $I-J$ are both contained in S. It is clear that the elements in these sets are odd integers. Since for all $0\le t<n/2$, $0\le j<b$ and $0\le i<b^{\prime }$ we have the following:

$$\left(2b\right(2t+1)+2bni)+(2j+1)\le 2b\left(2\right(n/2-1)+1+ni)+(2b+1)=2nbi+2nb+1$$

$$(2b(2t+1)+2bni)+(2j+1)\le 2nb\left(\frac{(\ell -1)}{2nb}-1\right)+2nb+1=\ell .$$

Similarly, $1\le \left(2b\right(2t+1)+2bni)-(2j+1)$ for all $0\le t<n/2$, $0\le j<b$ and $0\le i<b^{\prime }$ since the minimum of I is $\left(2b\right(2t+1)+2bni)=2b$ when $t=i=0$ and the maximum of J is $2j+1=2b-1$ when $j=b-1$.    □

Note that we have chosen to partition I into only $n/2$ subsets because, when n cores are available, we assign $n/2$ cores to each of the two resultants. However, when computing point multiples during KPS we do have all n cores available for working on the $n/2$ subsets. Therefore, it will be convenient to further divide each $I_t$ into two subsets to obtain a total of n half-sized subsets. Note that some subsets will have fewer elements than others because $\#I$ need not be a multiple of n. In practice, this is problematic because it is crucial for each core to perform exactly analogous tasks to guarantee performance. In order to achieve a balanced computation, we will add additional elements to the subsets that are lacking, ensuring that each subset contains the same number of elements, while ignoring redundant multiples in future steps. This idea is inspired by [22] and the parallelization technique proposed by Costello in [32] to compute the set ${\left[i\right]P}_{1\le i\le d}$. We provide an example to illustrate this approach.

Example 1.

Let $\ell =89$ and $n=4$. From the Lemma 3 we get

$$\begin{array}{cc}J:=\{2j+1:0\le j<3\},\hfill & \hfill {\displaystyle I=\bigcup _{t=0}^1{I}_t,}\end{array}$$

$$\begin{array}{cc}{I}_t:=\{2b((2t+1)+ni):0\le i<3\},\hfill & \hfill for all 0\le t<2.\end{array}$$

So $I_0:=\{2b,10b,18b\}$ and $I_1:=\{6b,14b,22b\}$. Instead of this, we can consider the following sets:

$$\begin{array}{cc}I{I}_t:=\{2b((2t+1)+2ni):0\le i<2\}\hfill & \hfill forall0\le t<4.\end{array}$$

So $I{I}_0:=\{2b,18b\}$, $I{I}_1:=\{6b,22b\}$, $I{I}_2:=\{10b,26b\}$ and $I{I}_3:=\{14b,30b\}$.

We will now demonstrate that the partitioning of the set I, without taking into account the additional elements added to ensure equal-sized subsets, still forms an index system.

Proposition 1.

Let us assume the same conditions of the Lemma 3. Let $c=\lceil b^{\prime }/2\rceil >2$. If

$$\begin{array}{cc}{\displaystyle I^{\prime }:=\bigcup _{t=0}^{n-1}I{I}_t\setminus \bigcup _{t=n/2}^{n-1}\left\{w_t\right\}}\hfill & \hfill with{w}_t:=2b(2t+1+2n(c-1))\left(b^{\prime }mod2\right),\\ \hfill & \hfill andI{I}_t:=\{2b(2t+1)+bn\left(2i\right):0\le i<c\}\\ \hfill & \hfill forall0\le t<n,\end{array}$$

$$\begin{array}{cc}I{I}_t:=\{2b(2t+1)+bn\left(2i\right):0\le i<c\},\hfill & \hfill forall0\le t<n,\end{array}$$

$$\begin{array}{cc}{\displaystyle I^{\prime }:=\bigcup _{t=0}^{n-1}I{I}_t\setminus \bigcup _{t=n/2}^{n-1}\left\{w_t\right\},}\hfill & \hfill with{w}_t:=2b(2t+1+2n(c-1))\left(b^{\prime }mod2\right),\end{array}$$

then $I^{\prime }=I.$ In particular, ($I^{\prime },J$) is an index system for S.

Proof of Proposition 1.

If $b^{\prime }$ is an even integer, it is clear that $I{I}_t\subset I$, for all $0\le t<n$ and thus ${\displaystyle \bigcup _{t=n/2}^{n-1}\left\{w_t\right\}=\left\{0\right\}}$. Note that $0\le 2i\le 2c=b^{\prime }$ and considering that ${\displaystyle \bigcup _{t=0}^{n-1}I{I}_t}$ and I have the same cardinality, we conclude that they are equal. If $b^{\prime }$ is an odd integer, all $0\le t<n/2$ we have that $I{I}_t\subset I$. When $n/2\le t<n$, dropping the last element $w_t:=2b(2t+1+2n(c-1))$ in each $I{I}_{n/2},...,I{I}_{n-1}$, we have $I{I}_t\subset I$. Note that the cardinality of ${\displaystyle \bigcup _{t=n/2}^{n-1}\left\{w_t\right\}}$ is $n/2$ and the cardinality of ${\displaystyle \bigcup _{t=0}^{n-1}I{I}_t}$ is $nc=n((b^{\prime }+1)/2)=n{b}^{\prime }+n/2=\#I+n/2$.    □

Similarly, we consider a partition of $K=S\setminus (I\pm J)=\{2b{b}^{\prime }n+1,...,\ell -2,\ell \}$ into n sets $K_0,..,K_{n-1}$ with $d=\lceil \#K/n\rceil$ elements,

$$K_t=\{2b{b}^{\prime }n+2(t+1)+2nk:0\le k<d\},0\le t<n$$

However, when $\#K$ is not divisible by n, we will add one additional element to each set $K_{dmod\left(n\right)},...,K_{n-1}$, ensuring that each subset contains the same number of elements.

The additional elements that were added to balance the size of the sets are not taken into account when calculating the corresponding products for each set.

3.2. Parallelized Algorithms

We now describe the construction of parallelized versions of algorithms KPS, xISOG, and xEVAL assuming that n cores are available, which we refer to as KPS-Parallel, xISOG-Parallel, and xEVAL-Parallel.

The KPS-Parallel algorithm, presented as Algorithm 4, is based on the new index system from Proposition 1: instead of computing the sets $\mathcal{I}$ and $\mathcal{K}$ as in the original algorithm, KPS-Parallel computes the sets $\mathcal{I}{\mathcal{I}}_t$ and ${\mathcal{K}}_t$, as well as the polynomial product tree for $h_{I{I}_t}$, by assigning one core to each value of $0\le t<n$. Note that the computation of the resultants requires a reciprocal tree for each of the $h_{I_t}$, so the polynomials $\{h_{I{I}_t},0\le t<n\}$ are multiplied pair-wise to obtain the polynomials $\{h_{I_t},0\le t<n/2\}$. We then compute the residue trees, which are built from the root down: the root is computed for each tree using $n/2$ cores, and then the two subtrees of each tree are computed in parallel, using two cores per tree, for a total of n cores. Note that although the computations for the root of the product and reciprocal trees, corresponding to lines 11 and 12 of 4, should only be ran with $n/2$ cores, we instead have all cores running with half of them repeating the work over dummy arrays ${\mathcal{I}}_t$ for $t>n/2$. These computations are of course redundant, but ensure that the workload is balanced.

Algorithm 4 KPS-Parallel

For the algorithms xISOG-Parallel and xEVAL-Parallel, the set J is partitioned into n subsets $J_t$, and each core computes the polynomial product tree corresponding to one subset. The roots of the n trees are then multiplied together sequentially to obtain the full polynomial $E_{i,J}$. Next, each core computes one of the subresultants $Re{s}_Z(h_{I_t},E_{i,J})$ using the reciprocal trees from KPS-Parallel, and the subresultants are multiplied sequentially at the end to obtain the two principal resultants. As for the computations related to $\mathcal{K}$, it is also split into subsets ${\mathcal{K}}_t$ so that each core can compute a subproduct of $M_i$, and the subproducts are multiplied together at the end.

Algorithm 5 xISOG-Parallel

In the next subsection, we present a detailed analysis that estimates the cost of computing the image curve of an isogeny and the evaluation of a point using the procedures KPS-Parallel+ xISOG-Parallel+ xEVAL-Parallel.

Algorithm 6 xEVAL-Parallel

3.3. Cost Analysis of the Parallel Square-Root Vélu

In this section, we will provide an estimation of the cost that is expected when performing the KPS-Parallel, xISOG-Parallel, and xEVAL-Parallel procedures, which is analogous to the one presented in subSection 2.3 for the sequential algorithms. Recall that for our new index system from subSection 3.1 we have $\#J=\lfloor \sqrt{(\ell -1)/\left(2n\right)}\rfloor$, where ℓ is the degree of the isogeny and n represents the number of cores available, while each of the $I_t$ subsets has a size of $b^{\prime }:=\lfloor (\ell -1)/\left(2n(\#J)\right)\rfloor$ for a total size of $\#I=(n/2)\lfloor (\ell -1)/\left(2n\right(\#J\left)\right)\rfloor$. It follows that

$$\begin{array}{cc}\hfill \#K:=& \frac{\ell -1}{2}-2(\#I)(\#J),\hfill \end{array}$$

$$\begin{array}{cc}\hfill \#K\le & \frac{\ell -1}{2}-2\left(\sqrt{\frac{\ell -1}{2n}}-1\right)\left(\frac{n}{2}\left(\sqrt{\frac{\ell -1}{2n}}-1\right)\right),\hfill \end{array}$$

$$\begin{array}{cc}\hfill \#K=& \sqrt{2n(\ell -1)}-n.\hfill \end{array}$$

As before, we will ignore rounding errors and assume that $\#K$ takes the middle value of this range (neglecting the second term), so that $\#I\approx nb/2$, $\#J\approx b$, and $\#K\approx nb$ for $b=\sqrt{(\ell -1)/\left(2n\right)}$.

Proposition 2.

The cost of computing the parallel AlgorithmKPS-Parallelwith n cores is

$$\begin{array}{c}\hfill Cos{t}_{\mathtt{KPS}-\mathtt{Parallel}}(b,n)=\frac{8}{3}{b}^{{log}_2\left(3\right)}+4{log}_2\left(b\right)+15b+12n-18,\end{array}$$

field multiplications.

Proof of Proposition 2.

In order to find all the point multiples in K, we first compute the first n multiples sequentially and then core t start from the t-th multiple and take steps of size n. This leads to wall-clock time of $n+\#K/n-1$ point operations, which involve 6 field multiplications each. A similar approach is taken for computing all the multiples in the $\mathcal{I}{\mathcal{I}}_t$ sets, while the ones in $\mathcal{J}$ are computed sequentially as this is the smallest of the sets. A product tree is constructed for each of the $h_{{II}_t}$ in parallel, and line 11 involves multiplying two polynomials of degree b/2. As for the reciprocal trees, we incur the cost of the root node and then just half of the remainder cost since there is one core working on each subtree, so its cost is (Reciprocal(b) + ReciprocalTree(b))/2, where Reciprocal(b) = b^log2(3) + b + 2log₂(b) − 2 is the cost of the root node alone. The total cost is then

$$\begin{array}{cc}\hfill Cos{t}_{\mathtt{KPS}-\mathtt{Parallel}}(b,n)=& 6\left(n+b-1\right)+6\left(n+b/2-1\right)+6b\hfill \\ \hfill & +ProductTre{e}_1(b/2)+PolyMul(b/2)\hfill \\ \hfill & +\left(Reciprocal\right(b)+ReciprocalTree(b\left)\right)/2,\hfill \end{array}$$

and the proposition follows after considering the cost functions in the Appendix A.

   □

We now focus on Algorithm 5:xISOG-Parallel which computes the coefficient of the image curve, and on Algorithm 6:xEVAL-Parallel which computes the isogeny evaluation at a specific point.

Proposition 3.

The cost of computing the parallel AlgorithmxISOG-Parallelwith n cores is

$$\begin{array}{cc}\hfill Cost& { }_{\mathtt{xISOG}-\mathtt{Parallel}}(b,n)=\hfill \\ \hfill & 6b^{{log}_2\left(3\right)}+3(n^2-n+2){\left(\frac{b}{n}\right)}^{{log}_2\left(3\right)}+3{log}_2\left(b\sqrt{n}\right)-b-\frac{b}{n}+\frac{5}{2}n+\frac{11}{2},\hfill \end{array}$$

field multiplications.

Proof of Proposition 3.

We begin by considering the cost of computing the polynomials $E_{0,J}$ and $E_{1,J}$ in Algorithm 5, where $\mathcal{J}$comes from KPS-Parallel. We partition $\mathcal{J}$ into the subsets ${\mathcal{J}}_0,\dots ,{\mathcal{J}}_{n-1}$ of size $(b/n)$ each, and compute one sub-product tree $E_{0,J_t}$ per subset $J_t$ concurrently. The cost of computing the factors of each subset is given by $5(b/n)$ field multiplications, then a product tree is computed for $b/n$ quadratic polynomials on each core, and the tree roots are multiplied together sequentially to obtain the complete polynomial $E_{0,J}$. This last step involves multiplying together n polynomials of degree $2b/n$ each, which we denote as $LinearProduc{t}_{2b/n}\left(n\right)$. An identical procedure is used to compute $E_{1,J}$ as well.

Using the residue trees $rtre{e}_{h_{I_0}},...,rtre{e}_{h_{I_{n/2-1}}}$ from KPS-Parallel, each core then compute a subresultant using a residue tree of size $b/n$, and the subresultants are multiplied sequentially at a cost of $n/2-1$ multiplications per resultant.

For each of the $M_i$, each core computes a subproduct of size b at a cost of $b-1$ multiplications, and then combining the subproducts takes another $n-1$ multiplications.

Finally, we use two cores to compute the ℓ-th power exponentiation in Algorithm  for the numerator and denominator concurrently at a cost of about $1.5log\ell \approx 3/2+3{log}_2\left(b\right)+3log\left(n\right)/2$, and the last few products take 10 more multiplications.

The total cost is then

$$\begin{array}{cc}\hfill Cos{t}_{\mathtt{xISOG}-\mathtt{Parallel}}(b,n)=& 5b/n+2\times ProducTre{e}_2(b/n)+2\times LinearProduc{t}_n(2b/n)\hfill \\ \hfill & +ResidueTree\left(b\right)+n/2-1\hfill \\ \hfill & +2\times (b+n-2)\hfill \\ \hfill & +3/2+3{log}_2\left(b\right)+3log\left(n\right)/2+10,\hfill \end{array}$$

and the proposition follows after considering the cost functions in the Appendix A.    □

Proposition 4.

The cost of computing the parallel AlgorithmxEVAL-Parallelwith n cores is

$$\begin{array}{ccc}\hfill Cost& { }_{\mathtt{xEVAL}-\mathtt{Parallel}}(b,n)=\hfill & \hfill 6b^{{log}_2\left(3\right)}+\frac{3}{2}(n^2-n+2){\left(\frac{b}{n}\right)}^{{log}_2\left(3\right)}+b+\frac{7b}{n}+\frac{5}{2}n.\end{array}$$

field multiplications.

Proof. 

The proof of the current proposition follows a similar structure to the previous proposition, with the main distinction lying in the polynomials $E_{0,J}$$E_{1,J}$$M_{0,t}$ and $M_{1,t}$. As before, the factors of $E_{0,J}$ can be computed with a cost of $10(b/n)$, and only one product tree is needed since $E_{1,J}$ is obtained at no additional cost. For each of $M_0$ and $M_1$, each core is still working with a subset of size b, but now there is a cost of b scalar multiplications for computing the factors in a subset, $b-1$ for the subproduct of each subset, and $n-1$ for combining the subproducts. There is no exponentiation in this case, and the final steps require only 6 additional multiplications. The total cost is then

$$\begin{array}{cc}\hfill Cos{t}_{\mathtt{xEVAL}-\mathtt{Parallel}}(b,n)=& 10b/n+ProducTre{e}_2(b/n)+LinearProduc{t}_n(2b/n)\hfill \\ \hfill & +ResidueTree\left(b\right)+n/2-1\hfill \\ \hfill & +2(2b+n-2)\hfill \\ \hfill & +6,\hfill \end{array}$$

and the proposition follows after considering the cost functions in the Appendix.    □

The next theorem summarizes our cost analyzis, and follows immediately from the previous propositions.

Theorem 1.

The expected total cost of the algorithmsKPS-Parallel,xISOG-Paralleland twicexEVAL-Parallelexpressed in field multiplications, assuming that two points need to be pushed through the isogeny, is given by

$$\begin{array}{cc}\hfill Cos{t}_{\mathtt{V}é\mathtt{lu}-\mathtt{Parallel}}(b,n)=& Cos{t}_{\mathtt{KPS}-\mathtt{Parallel}}(b,n)+Cos{t}_{\mathtt{xISOG}-\mathtt{Parallel}}(b,n)\hfill \\ \hfill & +2\times Cos{t}_{\mathtt{xEVAL}-\mathtt{Parallel}}(b,n)\hfill \\ \hfill =& \frac{62}{3}{b}^{{log}_2\left(3\right)}+6(n^2-n+2){\left(\frac{b}{n}\right)}^{{log}_2\left(3\right)}+7{log}_2\left(b\right)\hfill \\ \hfill & +\frac{13b}{n}+16b+{\displaystyle \frac{3}{2}}{log}_2\left(n\right)+\frac{39}{2}n-{\displaystyle \frac{25}{2}}.\hfill \end{array}$$

(6)

Remark 3.

While it may be tempting to compare (6) directly to (5), a direct comparison would be incorrect since the definitions of b (the size of the set J) are different in each case. The fair comparison would be in terms of ℓ, where (5) has $b=\sqrt{(\ell -1)}/2$, whereas (6) has $b=\sqrt{(\ell -1)/\left(2n\right)}$.

4. Experimental Results

We implemented our parallel version of the square-root Vélu algorithm in the C programming language, leveraging OpenMP for parallelization. This implementation has been integrated into a fork of the SQALE’d CSIDH library [12], where we can measure the savings due to parallelism for isogenies of each prime degree $\ell |p+1$. In order to visualize the progression of the performance with the isogeny degree, we have incorporated a new 1792-bit prime

$$\begin{array}{cc}\hfill p_{test}:=0\mathrm{x}& 4\mathrm{FA}4\mathrm{E}8\mathrm{C}57\mathrm{C}4\mathrm{EFF}02\mathrm{D}0650\mathrm{FE}3\mathrm{AFC}0413536\mathrm{E}72253101645\mathrm{B}1387\mathrm{DA}2\mathrm{DD}519\mathrm{C}17\mathrm{FBEC}\dots \hfill \\ \hfill & 69843\mathrm{C}04\mathrm{DEAEA}2\mathrm{DB}59\mathrm{CDDDA}7876\mathrm{B}514\mathrm{C}101\mathrm{A}1\mathrm{DF}0\mathrm{D}96778\mathrm{BD}72\mathrm{A}3\mathrm{C}51844\mathrm{BB}0196\dots \hfill \\ \hfill & \mathrm{F}73\mathrm{F}1\mathrm{DDBFEC}980\mathrm{A}4\mathrm{BB}3\mathrm{B}200\mathrm{A}4\mathrm{E}618\mathrm{C}54621\mathrm{ADB}35\mathrm{B}5\mathrm{E}4\mathrm{B}0545\mathrm{F}5\mathrm{BE}025\mathrm{D}63\dots \hfill \\ \hfill & \mathrm{BC}914\mathrm{AB}11\mathrm{A}882\mathrm{AD}78\mathrm{B}6203\mathrm{C}57\mathrm{A}31031\mathrm{B}98\mathrm{B}6\mathrm{C}104\mathrm{DC}99\mathrm{AC}9\mathrm{A}4532\mathrm{DEC}0\mathrm{C}293\dots \hfill \\ \hfill & 0\mathrm{F}8\mathrm{AE}51\mathrm{B}008\mathrm{E}4\mathrm{BA}6\mathrm{D}26\mathrm{E}56\mathrm{C}736\mathrm{D}3\mathrm{C}067\mathrm{C}8\mathrm{F}2\mathrm{DFDF}7\mathrm{F}8206\mathrm{B}444\mathrm{A}42\mathrm{D}39\mathrm{E}0\mathrm{F}4\mathrm{D}82\mathrm{FF}\dots \hfill \\ \hfill & 3\mathrm{FD}0\mathrm{EB}1\mathrm{DF}44\mathrm{B}31\mathrm{DDCDE}876\mathrm{E}658489\mathrm{E}1\mathrm{CA}359\mathrm{DAF}5868\mathrm{A}6\mathrm{C}22\mathrm{E}8455\mathrm{B}4\mathrm{A}4\mathrm{F}7679\mathrm{F}6\dots \hfill \\ \hfill & 2\mathrm{B}0\mathrm{C}30\mathrm{D}8883\mathrm{D}2\mathrm{B}79931\mathrm{C}19\mathrm{E}4737\mathrm{C}3\mathrm{CC}33056461\mathrm{E}9\mathrm{C}96\mathrm{A}175\mathrm{D}94\mathrm{B}594\mathrm{B}2\mathrm{A}9\mathrm{EAB}4\mathrm{B}6\mathrm{B}6303,\hfill \end{array}$$

where $p_{test}+1=4{\prod }_{i=0}^{107}{\ell }_i$ for odd primes ${\ell }_i$ that are roughly evenly spaced between ${\ell }_0=19$ and ${\ell }_{107}=321193$ (the degree of the isogeny in the attack of [5]). This prime was chosen so that isogenies of each degree ${\ell }_i$ can be evaluated using only ${\mathbb{F}}_{p_{test}}$ arithmetic (for which we also provide an assembly code implementation). As previously mentioned, this model most closely resembles the scenario in CSIDH as opposed to applications where isogenies are performed over quadratic fields or larger extensions. However, we make this choice for concreteness while still expecting our results to extrapolate well to other scenarios since the main savings come from reducing the number of field multiplications by a factor that is agnostic to the field itself. Our library can be found at https://github.com/TheParallelSquarerootVeluTeam/Parallel-Squareroot-Velu. The benchmark results presented in this section are from experiments conducted using an Intel(R) Xeon(R) Gold 6230R processor equipped with 26 physical cores, operating at $2.10\mathrm{GHz}$. Turbo boost and Hyper-threading technologies were disabled during these experiments.

Our experimental results are illustrated in Figure 4. We measured the total computational time for executing KPS-Parallel, xISOG-Parallel, and two iterations of xEVAL-Parallel for each of the prime degrees ${\ell }_i$ using 2, 4 and 8 cores, and compared to the computational time of the original sequential implementation of [12]. In Figure 4 we compare the observed speedup factors against the theoretical speedup using (6) and (5). We observe there is a more important difference for the eight-core implementation. Moreover, we notice that in situations where the degree of the isogeny is not sufficiently large, the overhead associated with synchronization and the computational steps involved tends to offset any potential performance gains achieved by employing eight cores. However, the speedup factor can be seen to achieve levels close to the asymptotic value starting from $\ell \approx 10000$.

To showcase our results in a more practical setting, we also implemented and benchmarked our parallelized square-root Vélu algorithms in the context of SQALE’d CSIDH-9216. In its dummy-free variant, SQALE’d CSIDH-9216 performs the operations corresponding to KPS + xISOG + 2xEVAL exactly once for each of the 333 smallest odd primes (excluding 263), over a 9216-bit prime field. The experimental timings for each of the three routines, summed over all the degrees, are shown in Table 1. Our two-core implementation achieves speedup factors of 1.23, 1.70, and 1.92 for KPS-Parallel, xISOG-Parallel, and xEVAL-Parallel, respectively, while for the four-core implementation we obtained acceleration factors of 1.42, 2.27, and 3.02. While these accelerations span the entirety of the Vélu-related computations in the protocol, it is important to stress that they cannot be translated easily into an acceleration for the protocol as a whole, since the protocol includes other non-negligible computations (most notably elliptic curve point scalar multiplications).

5. Discussion

Since KPS-Parallel performs some of its computational tasks sequentially, while others only benefit from half of the total number of cores, it is not surprising that its speedup factor is the lowest out of the three square-root Vélu routines. Nonetheless, as illustrated by Table 1, this is not too important as its total cost is comparatively small. Fortunately, the largest accelerations come from the xEVAL-Parallel, which apart from being the most expensive of the three is often required to be evaluated multiple times for different points.

Both our experimental and theoretical results show that parallelizing with at least two cores yields strong improvements, which are important considering that they can be combined with other improvements such as vectorized field arithmetic. On the downside, it is also noticeable that the speedup factors in our implementations do not scale well when transitioning to four- and eight-core implementations. From eight cores onward, even the maximum speedup derived from the theoretic estimates suggests that the utilization that can be achieved may not be attractive for practical applications. The intuition behind these diminishing returns on the parallelization effectiveness is easy to explain, as the size of the resultants in algorithm 5 and 6 only decreases as $1/\sqrt{n}$.