# Accelerating the Delfs-Galbraith Algorithm with Fast Subfield Root Detection

Maria Corte-Real Santos, Craig Costello, Jia Shi

► To cite this version:

Maria Corte-Real Santos, Craig Costello, Jia Shi. Accelerating the Delfs-Galbraith Algorithm with Fast Subfield Root Detection. CRYPTO 2022, Aug 2022, Santa Barbara, United States. pp.285-314, 10.1007/978-3-031-15982-4_10. hal-05446980

HAL Id: hal-05446980

https://hal.science/hal-05446980v1

Submitted on 7 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Accelerating the Delfs–Galbraith algorithm
with fast subfield root detection

Maria Corte-Real Santos<fnref n="1" />*, Craig Costello<fnref n="2" />, and Jia Shi<fnref n="3" />†

**Abstract.** We give a new algorithm for finding an isogeny from a given supersingular elliptic curve $E/\mathbb{F}_{p^2}$ to a subfield elliptic curve $E'/\mathbb{F}_p$, which is the bottleneck step of the Delfs–Galbraith algorithm for the general supersingular isogeny problem. Our core ingredient is a novel method of rapidly determining whether a polynomial $f \in L[X]$ has any roots in a subfield $K \subset L$, while avoiding expensive root-finding algorithms. In the special case when $f = \Phi_{\ell,p}(X, j) \in \mathbb{F}_{p^2}[X]$, i.e., when $f$ is the $\ell$-th modular polynomial evaluated at a supersingular $j$-invariant, this provides a means of efficiently determining whether there is an $\ell$-isogeny connecting the corresponding elliptic curve to a subfield curve. Together with the traditional Delfs–Galbraith walk, inspecting many $\ell$-isogenous neighbours in this way allows us to search through a larger proportion of the supersingular set per unit of time. Though the asymptotic $\bar{O}(p^{1/2})$ complexity of our improved algorithm remains unchanged from that of the original Delfs–Galbraith algorithm, our theoretical analysis and practical implementation both show a significant reduction in the runtime of the subfield search. This sheds new light on the concrete hardness of the general supersingular isogeny problem (i.e. the foundational problem underlying isogeny-based cryptography), and has immediate implications on the bit-security of schemes like B-SIDH and SQISign for which Delfs–Galbraith is the best known classical attack.

**Keywords:** Isogeny-based cryptography, supersingular isogeny problem, Delfs–Galbraith algorithm.

# 1 Introduction

In its most general form, the *supersingular isogeny problem* asks to find an isogeny

$$ \phi: E_1 \to E_2 $$

<footnote n="1">¹ University College London, London, UK
maria.santos.20@ucl.ac.uk</footnote>

<footnote n="2">² Microsoft Research, Redmond, USA
craigco@microsoft.com</footnote>

<footnote n="3">³ University of Waterloo, Waterloo, Canada
j96shi@uwaterloo.ca</footnote>

<footnote>*Supported by EPSRC grant EP/S022503/1.</footnote>

<footnote>†Part of this work was done while Jia was an intern at Microsoft Research.</footnote>
between two given supersingular curves, $E_1/\bar{\mathbb{F}}_p$ and $E_2/\bar{\mathbb{F}}_p$. We emphasize that this is the general problem, where we do not assume knowledge of the degree of the isogeny, or any torsion point information. The best known classical attack against the supersingular isogeny problem is the Delfs-Galbraith algorithm [13], which, for two curves $E_1$ and $E_2$ defined over $\mathbb{F}_{p^2}$, has two steps. The first step computes random walks in the $\ell$-isogeny graph (for some choice of $\ell$) to find isogenies $\phi_1: E_1 \to E'_1$ and $\phi_2: E_2 \to E'_2$, such that $E'_1/\mathbb{F}_p$ and $E'_2/\mathbb{F}_p$ are *subfield curves*. There are around $\lfloor p/12 \rfloor$ supersingular curves up to isomorphism and $O(p^{1/2})$ of them are subfield curves, therefore this step runs in $\tilde{O}(p^{1/2})$ bit operations. The second step searches for a *subfield isogeny* $\phi': E'_1 \to E'_2$ that connects $\phi_1$ and $\phi_2$, and it requires $\tilde{O}(p^{1/4})$ bit operations [13]. It follows that the entire algorithm runs in $\tilde{O}(p^{1/2})$ operations on average, with the cost dominated by the first step, i.e., the search for paths to subfield curves.

**Solver.** To our knowledge, a precise complexity analysis of the Delfs-Galbraith algorithm has not been conducted. We fill this gap by presenting an optimised implementation of the Delfs-Galbraith algorithm, called **Solver**, and conducting experiments over many thousands of instances of the subfield search problem to determine its concrete complexity. Though **Solver** finds the full path, we focus on the optimisation and complexity of the bottleneck step: finding subfield curves. These optimisations include:

- *Choice of $\ell$.* In their high-level description of the algorithm, Delfs and Galbraith do not specify which $\ell$-isogeny graph to walk in. Framing the problem of taking a step in the $\ell$-isogeny graph as computing the roots of a polynomial of degree $\ell$, in **Solver** we chose the simplest and most efficient choice: $\ell = 2$.

- *Fast square root finding in $\mathbb{F}_{p^2}$.* We use the techniques presented in [25, §5.3] to construct an optimised algorithm for finding square roots in $\mathbb{F}_{p^2}$, which only requires two $\mathbb{F}_p$ exponentiations and a few $\mathbb{F}_p$ multiplications and additions.

- *Random walks in the 2-isogeny graph.* We implement a depth-first search to find subfield nodes in the 2-isogeny graph and give a precise complexity analysis on the number of $\mathbb{F}_p$ operations required.

**SuperSolver.** The main contribution of this paper is a new state-of-the-art algorithm for solving the general supersingular isogeny problem, called **SuperSolver**. This is a variant of the Delfs-Galbraith algorithm that exploits a combination of our new *subfield root detection* algorithm and the use of modular polynomials. We show that we can efficiently determine whether a polynomial $f \in L[X]$ has a root in a subfield $K \subset L$, without finding any roots explicitly. Though this algorithm works for general fields and polynomials (and may be of use in other contexts), we apply it to the case where $f = \Phi_{\ell,p}(X, j) \in \mathbb{F}_{p^2}[X]$, i.e., where $f$ is the $\ell$-th modular polynomial evaluated at a supersingular $j$-invariant. This provides a means of quickly determining whether there is an $\ell$-isogeny connecting
the corresponding elliptic curve to a subfield curve: we develop this Neighbour-InFp subroutine in Section 4, and use it as the core of our SuperSolver algorithm in Section 5.

In Section 7, we conduct extensive experiments using both our Solver and SuperSolver libraries, all of which show that SuperSolver performs much faster than Solver. In Table 1, we give a taste of the types of improvements we see in searching for subfield nodes over supersingular sets of various sizes, taking a number of primes from the isogeny-based literature. These primes were specifically chosen because the Delfs-Galbraith algorithm for the general supersingular isogeny problem is the best known classical attack against the cryptosystems they target.

Our Solver and SuperSolver algorithms are written in Sage [30] and Python and can be found at

https://github.com/microsoft/SuperSolver.

**Cryptographic implications.** This paper has implications on the classical bit-security of any supersingular isogeny-based scheme for which the Delfs-Galbraith algorithm is the best known attack; this includes the key exchange scheme B-SIDH [11], the signature scheme in [16, §4], and the signature scheme SQISign [15]. For any proposed instantiation of such schemes, our SuperSolver suite allows the analysis in Section 7 to be conducted on input of any prime $p$, and determines a precise estimate on the number of operations required (on average) to solve the corresponding supersingular isogeny problem. This is especially accurate when the cardinality of the class group is known, which has recently been shown to be feasible for primes up to 512 bits [5]. On the other hand, we point out that the improvements in this paper have no direct impact on the classical security of SIDH [14] and SIKE [19]. Though the Delfs-Galbraith algorithm can be used to attack any supersingular isogeny-based cryptosystem, there are much faster *claw-finding* algorithms (see [14,1]) for solving the special instances of isogeny problems that arise in those schemes.

**Roadmap.** We give the preliminaries in Section 2. In Section 3, we present our optimised instantiation of the traditional Delfs-Galbraith algorithm, called Solver. In Section 4, we construct an efficient algorithm to detect whether a polynomial has a root in a subfield. We use this algorithm to build SuperSolver in Section 5. In Section 6, we present a worked example to highlight the differences between both algorithms, and in Section 7 we present a number of implementation results that illustrate the concrete improvements offered by SuperSolver.

## 2 Preliminaries

In this section we briefly set notation and give the requisite background for this paper. Readers familiar with the paragraph headings below are welcome to skip to the final two paragraphs.
<table><caption>**Table 1.** The number of nodes inspected per $10^8$ field multiplications and for primes targeting schemes where Delfs-Galbraith is the best known classical attack. The Solver column corresponds to optimised Delfs-Galbraith walks in $\mathcal{X}(\mathbb{F}_p, 2)$ – see Section 3. The SuperSolver columns correspond to enabling our fast subfield root detection algorithm with the three fastest sets of $\ell$'s (top to bottom) – see Section 5. We also give the approximate number of $\mathbb{F}_p$-multiplications per node inspected at each step, as computed during the precomputation phase that predicts which sets will perform fastest; for Solver we have used the empty set {}.</caption><thead><tr><th rowspan="2">Prime p</th><th colspan="2">Solver</th><th colspan="2">SuperSolver</th></tr><tr><th>Nodes inspected</th><th>F<sub>p</sub>-mults. per node</th><th>Fastest Sets</th><th>Nodes inspected</th><th>F<sub>p</sub>-mults. per node</th></tr></thead><tbody><tr><td rowspan="3">B-SIDH-p247 [11]</td><td rowspan="3">248,913</td><td rowspan="3">402</td><td>{3,5,7,9,11}</td><td>1,716,751</td><td>55.6</td></tr><tr><td>{3,5,7,9,11,13}</td><td>1,727,601</td><td>56.0</td></tr><tr><td>{3,5,7,11,13}</td><td>1,731,625</td><td>57.8</td></tr><tr><td rowspan="3">TwinSmooth-p250 [12]</td><td rowspan="3">233,507</td><td rowspan="3">427</td><td>{3,5,7,9,11}</td><td>1,680,337</td><td>59.1</td></tr><tr><td>{3,5,7,9,11,13}</td><td>1,699,825</td><td>59.1</td></tr><tr><td>{3,5,7,11,13}</td><td>1,697,769</td><td>59.8</td></tr><tr><td rowspan="3">SQISign-p256 [15]</td><td rowspan="3">248,915</td><td rowspan="3">403</td><td>{3,5,7,9,11}</td><td>1,716,751</td><td>55.7</td></tr><tr><td>{3,5,7,9,11,13}</td><td>1,727,601</td><td>55.9</td></tr><tr><td>{3,5,7,11,13}</td><td>1,731,625</td><td>57.8</td></tr><tr><td rowspan="3">TwinSmooth-p384 [12]</td><td rowspan="3">163,331</td><td rowspan="3">610</td><td>{3,5,7,9,11,13}</td><td>1,529,025</td><td>63.1</td></tr><tr><td>{3,5,7,9,11}</td><td>1,464,709</td><td>65.1</td></tr><tr><td>{3,5,7,9,11,13,17}</td><td>1,487,919</td><td>65.4</td></tr><tr><td rowspan="3">TwinSmooth-p512 [12]</td><td rowspan="3">127,511</td><td rowspan="3">784</td><td>{3,5,7,9,11,13}</td><td>1,397,761</td><td>69.1</td></tr><tr><td>{3,5,7,9,11,13,17}</td><td>1,391,645</td><td>70.0</td></tr><tr><td>{3,5,7,9,11,13,19}</td><td>1,351,509</td><td>72.1</td></tr></tbody></table>

**Modular polynomials.** We will use $\Phi_{\ell}(X, Y) \in \mathbb{Z}[X, Y]$ to denote the classical modular polynomial (see [29]) that parameterises pairs of elliptic curves with cyclic $\ell$-isogeny in terms of their $j$-invariants: $\Phi_{\ell}(j_1, j_2) = 0$ if and only if $j_1$ and $j_2$ are the $j$-invariants of $\ell$-isogenous elliptic curves. Readers unfamiliar with modular polynomials are encouraged to look at Sutherland's database$^4$, which contains $\Phi_{\ell}(X, Y)$ for all $\ell \le 300$ and for all primes $\ell \le 1000$. The polynomial $\Phi_{\ell}$ is symmetric in $X$ and $Y$, i.e., $\Phi_{\ell}(X, Y) = \Phi_{\ell}(Y, X)$, and if $\ell = \prod_{i=1}^n \ell_i^{e_i}$ is

<footnote>$^4$See [28], a database computed using techniques from various joint works of his [29,6].</footnote>
ℓ's prime decomposition, the degree of $\Phi_{\ell}(X, Y)$ in both X and Y is

$$N_{\ell} := \deg(\Phi_{\ell}(X, Y)) = \prod_{i=1}^{n} (\ell_i + 1) \ell_i^{e_i - 1}. \qquad (1)$$

The difficulty in computing $\Phi_{\ell}(X, Y)$ is in the size, rather than the number, of its coefficients. As discussed in [29], storing $\Phi_{\ell}(X, Y)$ requires $O(\ell^3 \log \ell)$ bits, which corresponds to several gigabytes for $\ell \approx 1000$ and many terabytes for $\ell \approx 10^4$. Fortunately, for our purposes, the modular polynomials already contained in Sutherland's database are more than sufficient. Moreover, we will be using them in the context of cryptanalysing instances of the supersingular isogeny problem over a fixed finite field $\mathbb{F}_{p^2}$, meaning we can reduce all of the large coefficients modulo $p$ as a precomputation. Indeed, even before the target $j$-invariants are known, $\Phi_{\ell}(X, Y) \in \mathbb{Z}[X, Y]$ will be preprocessed into

$$\Phi_{\ell,p}(X, Y) \in \mathbb{F}_p[X, Y],$$

where we note the additional subscript, defined by reducing all coefficients of $\Phi_{\ell}(X, Y)$ modulo $p$. By the symmetry of $\Phi_{\ell}(X, Y)$, this means we must store around $N_{\ell}^2/2$ coefficients in $\mathbb{F}_p$, requiring only $O(\ell^2 \log p)$ bits.

**Supersingular isogeny graphs.** Following [13, §1], let $p > 3$ be a prime and let $S_{p^2}$ denote the set of all supersingular $j$-invariants in $\mathbb{F}_{p^2}$. The number of such $j$-invariants is $\#S_{p^2} = \lfloor p/12 \rfloor + b$, where $b \in \{0, 1, 2\}$ is determined by the value of $p$ mod 12 [27, Theorem V.4.1(c)]. For any positive integer $\ell$ with $p \nmid \ell$, we use $\mathcal{X}(\mathbb{F}_p, \ell)$ to denote the supersingular isogeny graph whose nodes correspond to the $j$-invariants in $S_{p^2}$ and whose edges are $\ell$-isogenies defined over $\mathbb{F}_p$. When $\ell$ is prime, these graphs are fully connected [23], and (with the possible exception of a few nodes) are $(\ell + 1)$-regular expander graphs that satisfy the Ramanujan property [24]. Crucial to both the Delfs-Galbraith algorithm and this paper is the subset $S_p$ of supersingular $j$-invariants defined over $\mathbb{F}_p$. The size of this set is $\#S_p = \tilde{O}(p^{1/2})$ [13, Equation 1], and since $\#S_{p^2} = O(p)$, the expected number of randomly chosen elements in $S_{p^2}$ we would have to take before finding one in $S_p$ is in $\tilde{O}(p^{1/2})$.

**The Delfs-Galbraith algorithm.** The Delfs-Galbraith paper largely focusses on the problem of finding an isogeny $\phi': E_1' \to E_2'$ between two supersingular curves, $E_1'/\mathbb{F}_p$ and $E_2'/\mathbb{F}_p$, whose $j$-invariants are in $S_p$. One of their main results is an algorithm [13, Algorithm 1] that computes such a $\phi'$ in $\tilde{O}(p^{1/4})$ bit operations. At the end of their paper [13, Section 4], they show how this can be used as a subroutine to give an algorithm for the general supersingular isogeny problem, which asks to find an isogeny

$$\phi: E_1 \to E_2$$

between two supersingular curves, $E_1/\mathbb{F}_{p^2}$ and $E_2/\mathbb{F}_{p^2}$, whose $j$-invariants are in $S_{p^2}$. The idea is to perform simple non-backtracking random walks in $\mathcal{X}(\mathbb{F}_p, \ell)$
until hitting an elliptic curve with a $j$-invariant defined over $\mathbb{F}_p$. Finding a walk from $E_1/\mathbb{F}_{p^2}$ to $E'_1/\mathbb{F}_p$ yields an isogeny $\psi_1: E_1 \to E'_1$, and finding a walk from $E_2/\mathbb{F}_{p^2}$ to $E'_2/\mathbb{F}_p$ yields an isogeny $\psi_2: E_2 \to E'_2$. A full isogeny $\phi: E_1 \to E_2$ is then found as the composition $\phi = (\hat{\psi}_2 \circ \phi' \circ \psi_1)$, where $\hat{\psi}_2: E'_2 \to E_2$ is the dual of $\psi_2$, and $\phi': E'_1 \to E'_2$ is the subfield isogeny above that can be computed in $\tilde{O}(p^{1/4})$ bit operations. The bottleneck in the Delfs-Galbraith algorithm is finding the paths from the curves with $j \in S_{p^2} \setminus S_p$ to the curves with $j \in S_p$. From the above discussion, the number of $j$-invariants in $S_{p^2}$ we expect to search over before finding one in $S_p$ is $\tilde{O}(p^{1/2})$. Following [13, Section 4], the steps taken in $\mathcal{X}(\tilde{\mathbb{F}}_p, \ell)$ are non-backtracking, meaning that one stores the current $j$-invariant, $j_c$, and the previous $j$-invariant $j_p$. To take the next step, one then chooses one of the $N_\ell - 1$ roots (see Equation 1) of

$$ \Phi_{\ell}(X, j_c)/(X - j_p) $$

at random. Since $\ell$ and $N_\ell$ are fixed and small, it follows that the asymptotic complexity of the search for subfield $j$-invariants is $\tilde{O}(p^{1/2})$. Before presenting our improved search for subfield $j$-invariants, in Section 3 we present an optimised version of this algorithm, and subsequently replace the $\tilde{O}$ above with a precise, concrete complexity.

**Factoring polynomials in finite fields.** Let $f(X) \in \mathbb{F}_q[X]$ be a monic polynomial of degree $\ell$ with $q = p^k$ for a prime $p$, and for the purposes of this paper, assume that $p$ is very large (i.e., cryptographically sized) and $\ell$ is relatively small (i.e., $\ell < 100$). The literature contains a number of methods for finding the irreducible factors of $f$ in $\mathbb{F}_q[x]$, and we briefly mention the most applicable and well-known algorithms for our scenario. Berlekamp's algorithm [4] factors $f$ using an expected number of $O(\ell^3 + \ell^2 \log \ell \log q)$ operations in $\mathbb{F}_q$ [26, Theorem 20.12]. This appears to be superior to the Cantor-Zassenhaus algorithm [8], which uses an expected number of $O(\ell^3 \log q)$ operations in $\mathbb{F}_q$ [26, Theorem 20.9], however one can take advantage of certain time-memory trade-offs to implement Cantor-Zassenhaus so that it requires $O(\ell^3 + \ell^2 \log q)$ operations in $\mathbb{F}_q$ [26, Exercise 20.13]. Note that both of these big-$O$ complexities hide a number of subtleties, that $\mathbb{F}_q$-inversions are included as $\mathbb{F}_q$ operations, and moreover that both of these algorithms are probabilistic. Their deterministic variants have worse complexities [26, §20.6].

**Polynomial GCD.** Euclid's integer GCD algorithm is easily adapted to compute polynomial GCD's [26, §17.3]. Computing the GCD of two polynomials $g, h \in \mathbb{F}_q[x]$ requires $O(\deg(g) \cdot \deg(h))$ operations in $\mathbb{F}_q$. Again, here each $\mathbb{F}_q$ inversion is counted as an $\mathbb{F}_q$ operation. In order to make our algorithms run as fast as possible, one of the necessary subroutines we derive in Section 4 is an inversion-free polynomial GCD algorithm, for which we state a tight upper bound on the concrete complexity.
**Measuring complexity.** Throughout this paper we will avoid stating asymptotic (i.e., big-*O*-style) complexities in favour of stating concrete ones. One of our goals in Section 3 is to replace the $\tilde{O}(p^{1/2})$ complexity of the original Delfs-Galbraith algorithm with a closed formula that can be used to give precise estimates on the classical security of the relevant cryptographic instantiations. We will use the metric of $\mathbb{F}_p$ multiplications as convention, noting that it is relatively straightforward to convert this into a more fine-grained metric (e.g. bit operations, machine operations, cycle counts, gate counts, circuit depth, etc.) depending on the context and on the implementation of the $\mathbb{F}_p$ arithmetic. For simplicity, we will count $\mathbb{F}_p$ squarings as multiplications and ignore additions. We justify this by noting that, roughly speaking, the ratio of multiplications to additions in all of the algorithms in this work are similar, and the complexity of $\mathbb{F}_p$ additions have a minimal impact on any of the aforementioned metrics.

**Subfield search complexity determines concrete bit security.** Both the Solver implementation detailed in Section 3 and the SuperSolver implementation detailed in Section 5 solve all instances of the general supersingular isogeny problem. On input of any prime $p$ and any two supersingular $j$-invariants in $S_{p^2}$, both implementations will always terminate with an isogeny that solves the corresponding problem. We emphasise that henceforth our sole focus is on the $\tilde{O}(p^{1/2})$ subfield search phase of the Delfs-Galbraith algorithm. Finding a path between subfield nodes requires $\tilde{O}(p^{1/4})$ operations, which is negligible in both the asymptotic sense and in the sense of obtaining cryptographic security estimates. To see this, suppose the asymptotic $\tilde{O}(p^{1/2})$ complexity of the first phase is replaced by a concrete complexity of $c_p \cdot p^{1/2}$, and the asymptotic $\tilde{O}(p^{1/4})$ complexity of the second phase is replaced by a concrete complexity of $d_p \cdot p^{1/4}$, where $c_p$ and $d_p$ are polynomials in $\log p$. The total complexity of the Delfs-Galbraith algorithm is then

$$c_p \cdot p^{1/2} + d_p \cdot p^{1/4}.$$

For primes of cryptographic size, small changes in $c_p$ have an immediate influence on the total runtime of the algorithm, while much larger changes in $d_p$ will not play a part in the bit security of the problem. For $p > 2^{200}$, a factor 2 change in $c_p$ changes the bit security of the problem by 1, while $d_p$ would have to change by a factor of at least $2^{50}$ to have the same impact on the bit security.

### 3 Solver: optimised Delfs-Galbraith subfield searching in $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$

Recall from the previous section that the non-backtracking walks in $\mathcal{X}(\bar{\mathbb{F}}_p, \ell)$ store the current $j$-invariant, $j_c$, and the previous $j$-invariant $j_p$, and then take a step in $\mathcal{X}(\bar{\mathbb{F}}_p, \ell)$ by choosing one of the $N_\ell - 1$ roots of $\Phi_\ell(X, j_c)/(X - j_p)$. In determining the asymptotic $\tilde{O}(p^{1/2})$ complexity of these walks, Delfs and Galbraith did not need to analyse the cost of a single step. However, to set the
stage for our improved search in Section 5, we must optimise this process and determine its concrete cost. The first parameter that must be specified is $\ell$, i.e., the isogeny graph to walk around in. Considering both Equation (1) and the complexity of the factorisation algorithms in Section 2, we chose $\ell = 2$ to obtain most efficient and simplest choice where we are able to take advantage of fast explicit methods for computing square roots in $\mathbb{F}_{p^2}$.

**Scott's fast square roots in $\mathbb{F}_{p^2}$**. Optimal computation of square roots in extension fields of large characteristic requires careful attention to detail. A 2013 paper by Adj and Rodríguez-Henríquez [2] cost the process of computing square roots in $\mathbb{F}_{p^2}$ at two $\mathbb{F}_p$ residuosity tests, two $\mathbb{F}_p$ square roots, and one $\mathbb{F}_p$ inversion, for a total of five exponentiations in $\mathbb{F}_p$. In [25, §5.3], Scott shows that these operations can be combined in a clever way to significantly reduce this cost. The inputs into the Tonelli-Shanks $\mathbb{F}_p$ square root algorithm [22, Algorithm 3.34] can be tweaked in such a way that the two residuosity tests are absorbed into the two square roots. Moreover, he shows that most of the inversion cost can also be absorbed by application of Hamburg's combined 'square-root-and-inversion' trick [17]. This reduces the bulk of cost of an $\mathbb{F}_{p^2}$ square root from five $\mathbb{F}_p$ exponentiations to just two. In addition, there are a handful of $\mathbb{F}_p$ multiplications and additions that either update the Tonelli-Shanks outputs depending on the residuosity outcomes or collect and combine the results according to the "complex" formula in [25, §5.3]. We use this to construct a general square root algorithm in our implementation that is highly optimised with respect to the number of $\mathbb{F}_p$ operations it incurs<fnref n="5" />.

**Taking a step in $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$**. After stepping from $j_p \in \mathbb{F}_{p^2}$ to $j_c \in \mathbb{F}_{p^2}$, a non-backtracking walk in $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$ will step to one of two new nodes: $j_0$ and $j_1$. These are computed by solving the quadratic equation that arises from the modular polynomial $\Phi_\ell(X, Y)$ with $\ell = 2$:

$$ \Phi_2(X, Y) = -X^2Y^2 + X^3 + Y^3 + 1488 \cdot (X^2Y + Y^2X) - 162000 \cdot (X^2 + Y^2) + 40773375 \cdot XY + 8748000000 \cdot (X + Y) - 1574640000000000. $$

The three neighbours of $j_c$ in $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$ are $j_p, j_0$, and $j_1$, meaning that $\Phi_2(X, j_c)$ factorises as

$$ \Phi_2(X, j_c) = (X - j_p)(X - j_0)(X - j_1). $$

This yields a quadratic equation, whose solutions are $j_0, j_1$, defined by $X^2 + \alpha X + \beta = 0$, where

$$ \alpha = -j_c^2 + 1488 \cdot j_c + j_p - 162000, $$

$$ \beta = j_p^2 - j_c^2 j_p + 1488 \cdot (j_c^2 + j_c j_p) + 40773375 \cdot j_c - 162000 \cdot j_p + 8748000000. $$

<footnote n="5">⁵Note that the fixed exponentiations that take place in the calls to Tonelli-Shanks could be further optimised for a specific $p$ by tailoring a larger window or a different addition chain, but the impact (for our purposes and comparisons) of this improvement would be minor.</footnote>
Computing these coefficients costs a small, constant number of $\mathbb{F}_p$ operations, so the process of computing both $j_0$ and $j_1$ from $j_p$ and $j_c$ boils down to solving the quadratic equation, which essentially requires one $\mathbb{F}_{p^2}$ square root. Since this square root incurs two $\mathbb{F}_p$ exponentiations and a few additional $\mathbb{F}_p$ operations, it follows that the cost of computing each new $j \in S_{p^2}$ during the walks in $\mathcal{X}(\mathbb{F}_p, 2)$ is (on average) approximately one $\mathbb{F}_p$ exponentiation.

**The depth first search in $\mathcal{X}(\mathbb{F}_p, 2)$**. Repeating the process described above allows us to perform the search for subfield nodes using a depth first search in a binary tree with $d$ levels as follows. We write $j_{m,n}$ for the $n$-th node at level $m$, where $0 \le m \le d$ and $0 \le n \le 2^m - 1$. The first three levels are depicted in Figure 1. We initialise the root node $j_{0,0}$ as the target $j \in S_{p^2}$, and set $j_{1,0}$ and

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Levels 0, 1, and 2 of the binary tree in the depth first search of $\mathcal{X}(\mathbb{F}_p, 2)$.</figcaption></figure>

$j_{1,1}$ as two of its three neighbours<fnref n="6" /> in $\mathcal{X}(\mathbb{F}_p, 2)$. The depth first search starts by setting $j_c = j_{1,0}$ and $j_p = j_{0,0}$. We then solve the quadratic equation above to obtain $j_{2,0}$ and $j_{2,1}$, and repeat this procedure with $j_c = j_{i+1,0}$ and $j_p = j_{i,0}$ for $1 \le i \le d-1$ until the leftmost leaf $j_{d,0}$ is computed and the **path** stack is fully initialised as

$$ \mathbf{path} = [j_{0,0}, j_{1,0}, \dots, j_{d-1,0}, j_{d,0}]. $$

To avoid any waste, we also maintain a stack of the other solution to the quadratic equations that were computed along the way, which we call **sibling nodes**

$$ \mathbf{siblings} = [j_{1,1}, \dots, j_{d-1,1}, j_{d,1}]. $$

The algorithm then proceeds back up the levels by popping **path** until its last element is the root of a subtree that has not been checked in its entirety. At this point **siblings** is popped and pushed into **path**. When the last element of **path** is the root of a subtree that has not been exhausted, we initialise the process of solving quadratic equations, pushing one of the two solutions into

<footnote n="6">⁶Initially we do not have a $j_p$, so all three neighbours can be computed using generic root finding; our code does this during the setup phase.</footnote>
path and the other into siblings until path contains $d + 1$ elements. Each time the quadratic equation solver is called, the two roots (i.e., $j$-invariants) are immediately checked; if either of them lie in $\mathbb{F}_p$, it is added to path and the process is terminated. Otherwise, the process is repeated recursively until path = $[j_{0,0}]$, in which case the $2^{d+1} - 1$ nodes in the tree have been exhausted without finding a solution. To guarantee that a solution is found, one could increase $d$ and start again, but our code proceeds by simply storing the first (leftmost) leaf and its parent in separate memory so that the process can restart here and avoid recomputing any prior $j$'s. As Delfs and Galbraith point out, setting the depth $d = \frac{1}{2}\log_2 p$ should be enough. Since the number of nodes in the tree is $2^d$, increasing $d$ by $\epsilon$ makes the failure probability diminish by $1/2^\epsilon$. Setting $\epsilon = 10$ was sufficient in all of our experiments. Finally, as pointed out by Delfs and Galbraith in [13, §4], this process parallelises perfectly. For $P$ processors, one can simply compute a binary tree of depth $\lceil \log_2 P \rceil$ during setup and distribute $P$ of the leaf nodes as individual starting points.

**The concrete complexity of Delfs-Galbraith.** Table 2 reports on experiments conducted using Solver, the optimised instantiation of the traditional Delfs-Galbraith walk. For each bitlength between 21 and 40, we solved 10,000 instances of the subfield search. In each case we chose 100 random primes and, for each prime, 100 pseudo-random $j$-invariants in $S_{p^2}$. The numbers in each column report the averages (as base-2 logarithms) of these search complexities.

<table><caption>**Table 2.** The concrete cost of the subfield search phase of the Delfs-Galbraith over small fields of various bitlengths. Further explanation in text.</caption><thead><tr><th>Bitlengths of primes <i>p</i></th><th>21</th><th>22</th><th>23</th><th>24</th><th>25</th><th>26</th><th>27</th><th>28</th><th>29</th><th>30</th></tr></thead><tbody><tr><td>Av. number of nodes visited</td><td>8.8</td><td>9.4</td><td>10.0</td><td>10.3</td><td>10.9</td><td>11.4</td><td>11.9</td><td>12.3</td><td>13.1</td><td>13.5</td></tr><tr><td>Av. number of <i>F</i><sub><i>p</i></sub>-multiplications</td><td>14.5</td><td>15.0</td><td>15.7</td><td>16.0</td><td>16.7</td><td>17.2</td><td>17.8</td><td>18.2</td><td>19.0</td><td>19.5</td></tr></tbody></table><table><thead><tr><th>Bitlengths of primes <i>p</i></th><th>31</th><th>32</th><th>33</th><th>34</th><th>35</th><th>36</th><th>37</th><th>38</th><th>39</th><th>40</th></tr></thead><tbody><tr><td>Av. number of nodes visited</td><td>13.5</td><td>14.2</td><td>14.7</td><td>15.3</td><td>15.8</td><td>16.3</td><td>17.1</td><td>17.3</td><td>17.6</td><td>18.1</td></tr><tr><td>Av. number of <i>F</i><sub><i>p</i></sub>-multiplications</td><td>19.5</td><td>20.5</td><td>20.8</td><td>21.3</td><td>21.9</td><td>22.4</td><td>23.2</td><td>23.6</td><td>24.1</td><td>24.6</td></tr></tbody></table>

In all cases the number of $\mathbb{F}_p$ multiplications is found to be

$$ \#(\mathbb{F}_p \text{ mul.s.}) = c \cdot \sqrt{p} \cdot \log_2 p, $$

with $0.75 \le c \le 1.05$. In Section 7, we shed more light on the concrete complexity of both Solver and SuperSolver.

**Remark 1 (Vélu's formulas).** There is no traditional elliptic curve arithmetic found in either Solver or SuperSolver. All of the steps taken within $\mathcal{X}(\mathbb{F}_p, 2)$
and the rapid inspections conducted in $\mathcal{X}(\bar{\mathbb{F}}_p, \ell)$ use the modular polynomials. We point out there may be specific instances of $p$ where one could perform walks faster than repeatedly solving the $\Phi_{2,p}(X, j)$ quadratic by, say, employing Vélu's formulas [31] with the optimal strategies of De Feo-Jao-Plût [14]. For example, with a prime $p = 2^e 3^f - 1$, the price of computing a $2^e$-isogeny (i.e., walking through $e$ nodes in $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$) in this way may be cheaper than the price of computing $e$ square roots in $\mathbb{F}_{p^2}$ (note that the latter reveals 2 nodes each time). However, we argue that these scenarios are likely to only exist for special instances of the supersingular isogeny problem that are geared towards cryptosystems like SIDH [14] and SIKE [19]. As discussed in Section 1, here there are claw-finding algorithms that are much faster than the Delfs-Galbraith algorithm (though the number of $\mathbb{F}_p$ operations required to compute an $\ell^e$-isogeny still grows with $p$, and therefore our fast subfield root detection would also be useful in that context). In the case of both general primes and the types of primes in Table 1, it is highly unlikely that using Vélu's formulas [31] will be competitive with the binary tree depth-first search in $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$; computing general $(\prod \ell_i^{e_i})$-isogenies from kernel elements is much more expensive than $\ell^e$-isogenies when $\ell \in \{2, 3\}$, and one travels through fewer nodes in $S_{p^2}$ per $(\prod \ell_i^{e_i})$-isogeny when the $\ell_i$ grow larger.

**Remark 2 (Radical isogenies).** Another alternative to solving the quadratic equation that arises from $\Phi_2(X, j_c)/(X - j_p)$ is to instead take steps in $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$ using formulas for radical isogenies [10]. For example, given a supersingular Montgomery curve parameterised as $E_A : y^2 = x^3 + Ax^2 + x$ or as $E_\alpha = x(x - \alpha)(x - 1/\alpha)$, one can compute non-backtracking chains of 2-isogenies as either $A \to A' \to A''...$, or as $\alpha \to \alpha' \to \alpha''...$, rather than computing the chain of $j$-invariants $j \to j' \to j''...$, as we do. Computing the next value in all of these chains requires one square root (which dominates the cost for primes of cryptographic size) and a small handful of additional field operations, the number of which depends on the choice of chain. In the case of computing the chains $A \to A' \to A''...$ or $\alpha \to \alpha' \to \alpha''...$, the number of additional operations are fewer (see [9] and [7]) than those which we incur using the modular polynomial, however we have not opted to exploit this minor speedup for the following reasons. Indeed, it is not true in general that $j(E_A) \in \mathbb{F}_p$ implies $A \in \mathbb{F}_p$ or that $j(E_\alpha) \in \mathbb{F}_p$ implies $\alpha \in \mathbb{F}_p$. Since $j(E_A) = 256(A^2 - 3)^3/(A^2 - 4)$, in general there are six values of $A$ corresponding to a given $j$. Similarly, since $j(E_\alpha) = 256(\alpha^4 - \alpha^2 + 1)^3/(\alpha^4(\alpha^2 - 1)^2)$, in general there are twelve values of $\alpha$ corresponding to a given $j$. For large primes it is typically the case that most (or all) of the $A$'s and $\alpha$'s corresponding to a given $j \in S_p$ are not defined over $\mathbb{F}_p$. Thus, if radical isogenies were used to compute chains of $\alpha$'s or $A$'s in the context of Delfs-Galbraith, we would need to compute a value that determines whether the corresponding $j$ lies in $\mathbb{F}_p$. We note that this can be achieved without inverting the denominators in the expressions for $j(E_\alpha)$ or $j(E_A)$, i.e., $(a+b\cdot\beta)/(c+d\cdot\beta)$ is in $\mathbb{F}_p$ if and only if $ad = bc$ for $a, b, c, d \in \mathbb{F}_p$ and $\mathbb{F}_{p^2} = \mathbb{F}_p(\beta)$. Thus, the original Delfs-Galbraith walk in $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$ is likely to save a small, fixed number of multiplications per 2-isogeny by computing chains of $A$'s or $\alpha$'s instead
of $j$'s. However, when invoking our fast subfield root detection in the sections that follow, it is critical (for Algorithm 2) that the $j$-invariants of each node are computed explicitly, so that the higher $\ell$-degree modular polynomials can be used to probe for $\ell$-isogenous subfield neighbours. This subsequent computation of the $j$-invariant seems to require an additional field exponentiation (we could not see a way to merge the square roots and inversions into one exponentiation in these instances), which would kill the potential advantage of radical isogenies in the optimised SuperSolver algorithm.

**Remark 3 (Alternative modular functions).** There are several well-known modular functions other than the $j$-function – see [29]. A natural question in the context of this paper is whether any such functions can be used to make the search for subfield nodes in supersingular isogeny graphs more efficient. For example, the modular polynomials for Weber's $f$-function [28] are the same degree as those of the $j$-function, but have much smaller coefficients, many of which are zero. If these more compact modular polynomials could be used in the same way as those for the $j$-function, the practical gains would be significant. However, their applicability in the context of SuperSolver appears to be hampered by reasons similar to those discussed in Remark 2. Weber's $f$ is related to $j$ via $j = (f^{24} - 16)^3/f^{24}$, meaning there can be as many as 72 $f$'s corresponding to a single $j$-invariant, and it is not true in general that given $j \in \mathbb{F}_p$, the corresponding $f \in \mathbb{F}_p$. Although this makes the Weber polynomials unreliable replacements in the context of the SuperSolver algorithm, our search for alternative modular functions that would be compatible with SuperSolver was far from exhaustive, and it is likely that the $j$-function is not optimal across all of them. We leave any further investigations in this direction as future work.

# 4 Fast subfield root detection

In this section we derive a method for determining whether a polynomial $f(X) = a_nX^n + \dots + a_1X + a_0 \in \mathbb{F}_{q^d}[X]$ with $d \ge 2$ has a root lying in the subfield $\mathbb{F}_q$, where $q$ is a power of prime $p$. Though this can be achieved by factoring the polynomial, the methods described in Section 2 become too costly for our purposes; the number of $\mathbb{F}_q$ operations required depends on the size of $q$, which hampers their relative efficiency as $q$ grows large. Our aim in this section is to detail a much faster algorithm that detects whether a root lies in a subfield and show that the number of $\mathbb{F}_q$ operations required by our algorithm only depends on the degree of $f$ and the degree of the extension $d$.

As the algorithms in this section may be of independent interest, we leave them as general as possible before specialising back to the application at hand in Section 5. The results up to Proposition 1 are presented for general finite field extensions of the form $\mathbb{F}_{q^d}/\mathbb{F}_q$, but we will later specialise to the quadratic extensions of prime fields, i.e., where $q = p$ and $d = 2$. The inversion-free GCD in Algorithm 1 is derived for an arbitrary polynomial ring $K[x]$, but we will only need to use it in $\mathbb{F}_p[x]$.
In this section, for a polynomial in $\mathbb{F}_{q^d}[X]$, we will reduce the problem of detecting a root in $\mathbb{F}_q$ to computing the greatest common divisor of $d$ related polynomials $g_1, ..., g_d$. In the case where $d > 2$, we will need to compute the GCD of more than two polynomials. This can be done by recursively computing the GCD of two polynomials and using the following identity:

$$ \gcd(g_1, g_2, ..., g_d) = \gcd(g_1, \gcd(g_2, ..., g_d)). \qquad (2) $$

We aim to minimise the number of $\mathbb{F}_q$ multiplications needed to compute the GCD and so we construct these polynomials so that they are defined over $\mathbb{F}_q$. To achieve this, we will need two results. The first is a theorem by Lidl and Niederreiter [21, Theorem 2.24].

**Theorem 1.** Let $F$ be a finite extension of a finite field $K$, both considered as vector spaces over $K$. Then the linear transformations from $F$ into $K$ are exactly the mappings $L_\beta(\alpha)$, for $\beta \in F$, where $L_\beta(\alpha) = \text{Tr}_{F/K}(\beta\alpha)$ for all $\alpha \in F$. Furthermore, we have $L_\beta \neq L_\gamma$ whenever $\beta, \gamma$ are distinct elements of $F$.

The second result we will need is the following lemma.

**Lemma 1.** For $n \in \mathbb{N}$, let $f_1, ..., f_n \in \mathbb{F}_{q^d}[X]$ be polynomials and $A \in \text{GL}_n(\mathbb{F}_{q^d})$. Defining $(g_1, ..., g_n) := A \cdot (f_1, ..., f_n)$, we have

$$ \gcd(f_1, ..., f_n) = \gcd(g_1, ..., g_n). $$

*Proof.* If a polynomial $h \in \mathbb{F}_{q^d}[X]$ divides $f_1, ..., f_n$, then $h$ divides any linear combination of the $f_1, ..., f_n$. Therefore, $h$ divides $g_1, ..., g_n$. Since $A$ is invertible, by swapping the roles of $g_i$ and $f_i$ we see that the converse holds. $\square$

We are now ready to present the main result of this section.

**Proposition 1.** For some $d \ge 2$, let $\pi$ be the $q$-power Frobenius endomorphism in $\text{Gal}(\mathbb{F}_{q^d}/\mathbb{F}_q)$ and consider a polynomial $f(X) = a_nX^n + ... + a_1X + a_0 \in \mathbb{F}_{q^d}[X]$. Let $\beta$ be a primitive element of the extension $\mathbb{F}_{q^d}/\mathbb{F}_q$, in the sense that the field extension is generated by a single element $\beta$, i.e., $\mathbb{F}_q(\beta) = \mathbb{F}_{q^d}$. For $i = 1, ..., d$, define the following polynomials over $\mathbb{F}_{q^d}$:

$$ g_i := \sum_{j=0}^{d-1} \pi^j(\beta^{i-1}f). $$

Then $g_i(X) \in \mathbb{F}_q[X]$, and $\gcd(g_1, ..., g_d)$ divides $f$. In particular, if $\gcd(g_1, ..., g_d)$ is of degree 1, then $f$ has a root in $\mathbb{F}_q$. Furthermore, if $\gcd(g_1, ..., g_d) = 1$, then $f(X)$ does not have any roots in $\mathbb{F}_q$.

*Proof.* Using the notation in Theorem 1, we have

$$ \begin{aligned} g_i(X) &= [(\beta^{i-1}a_n + \cdots + \pi^{d-1}(\beta^{i-1}a_n))X^n + \cdots + (\beta^{i-1}a_0 + \cdots + \pi^{d-1}(\beta^{i-1}a_0))] \\ &= \sum_{m=0}^{n} L_{\beta^{i-1}}(a_m)X^m. \end{aligned} $$
By Theorem 1, for all $i = 1, ..., d$ and $m = 0, ... n$, we have $L_{\beta^{i-1}}(a_m) \in \mathbb{F}_q$, implying that $g_i(X) \in \mathbb{F}_q[X]$. Setting $(d \times d)$ matrix $A$ to be

$$A = \begin{bmatrix} 1 & 1 & \dots & 1 \\ \beta & \pi(\beta) & \dots & \pi^{d-1}(\beta) \\ \vdots & \vdots & \ddots & \vdots \\ \beta^{d-1} & \pi(\beta^{d-1}) & \dots & \pi^{d-1}(\beta^{d-1}) \end{bmatrix} = \begin{bmatrix} 1 & 1 & \dots & 1 \\ \beta & \beta^q & \dots & \beta^{q^{d-1}} \\ \vdots & \vdots & \ddots & \vdots \\ \beta^{d-1} & (\beta^{d-1})^q & \dots & (\beta^{d-1})^{q^{d-1}} \end{bmatrix},$$

we have $(g_1, ..., g_d) := A \cdot (f, \pi(f), ..., \pi^{d-1}(f))$. As for Vandermonde matrices [18, §6.2], we find $\det(A) = \prod_{0 \le i < j \le d-1} (\beta^{q^j} - \beta^{q^i})$, which is non-zero for $\beta$ a primitive element of the extension $\mathbb{F}_{q^d}/\mathbb{F}_q$ and so $A \in \text{GL}_d(\mathbb{F}_{q^d})$. By Lemma 1, we have

$$\text{gcd}(f, \pi(f), ..., \pi^{d-1}(f)) = \text{gcd}(g_1, ..., g_d),$$

therefore $\text{gcd}(g_1, ..., g_d) | f$. If $\text{gcd}(g_1, ..., g_d)$ is of degree 1, then $(X - r) | f$ for some $r \in \mathbb{F}_q$, and so $f$ has a root in $\mathbb{F}_q$.

We further note that $\text{gcd}(f, \pi(f), ..., \pi^{d-1}(f))$, and therefore $\text{gcd}(g_1, ..., g_d)$, is precisely the largest divisor of $f$ that is defined over $\mathbb{F}_q$. As a result, if $\text{gcd}(g_1, ..., g_d) = 1$, then $f(X)$ does not have any roots in $\mathbb{F}_q$. $\square$

**Applying Proposition 1 to detect subfield nodes.** The proof of Proposition 1 tells us that $\text{gcd}(g_1, ..., g_d)$ is precisely the largest divisor of $f \in \mathbb{F}_{q^d}[X]$ that is defined over $\mathbb{F}_q[X]$. In our target application of searching for subfield nodes in large supersingular isogeny graphs, i.e., when $d = 2$ and $q = p$, we will most commonly encounter $\text{gcd}(g_1, g_2) = 1$, which immediately rules out subfield neighbours in the $\ell$-isogeny graph. Non-trivial GCD's will, with overwhelmingly high probability, be of degree 1 and reveal a single subfield node; this is why our implementation of Algorithm 1 below terminates and returns true when the degree of the GCD is 1.

For large supersingular isogeny graphs, the only way for the degree of $\text{gcd}(g_1, g_2)$ to be larger than 1 is when a given $j$-invariant is $\ell$-isogenous to multiple subfield nodes, or when a given $j$-invariant is $\ell$-isogenous to conjugate $j$-invariants in $\mathbb{F}_p^2$.$^7$

In our scenario where $d = 2$, we see that $\pi(\beta) + \beta = 0$, meaning that $\pi^k(\beta) = (-1)^k \beta$. As a result, to detect a subfield root, we compute $\text{gcd}(g_1, \beta g_2)$ where $g_1 = f + \pi(f)$ and $g_2 = f - \pi(f)$. In this case we do not need to calculate any more powers of $\beta$ and we only need to do one GCD computation.

**Inversion-free polynomial GCD.** To complete the detection of roots in a subfield, we must compute the GCD of polynomials in polynomial ring $K[X]$,

<footnote>$^7$A real-world attack should check any non-trivial GCD, since either of these scenarios are a win for the cryptanalyst; the latter case reveals information about the secret endomorphism ring of the target isomorphism class (see [20, §5.3]), and the former case gives multiple solutions to the subfield search problem.</footnote>
where $K$ is a field. In Algorithm 1, we modify Euclid's polynomial-adapted algorithm [26, §17.3] to compute the GCD of two polynomials $g, h \in K[X]$ while avoiding inversions in $K$. We use $LC(f)$ to denote the leading coefficient of the polynomial $f$. Note that, for the purposes of incorporating it into our target application of subfield searching in the next section, the algorithm outputs the boolean true when the GCD has degree 1 in $K[X]$.

<caption>**Algorithm 1 InvFreeGCD(): Inversion-free GCD**</caption>

**Input:** Polynomials $g, h \in K[X]$, such that $\deg g \ge \deg h$
1: Initialize $r, s \leftarrow LC(h) \cdot g, LC(g) \cdot h$
2: **while** $\deg r \ge 1$ and $r \ne s$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;$r \leftarrow r - X^{\deg r - \deg s} \cdot s$
4: &nbsp;&nbsp;&nbsp;&nbsp;$r, s \leftarrow LC(s) \cdot r, LC(r) \cdot s$
5: &nbsp;&nbsp;&nbsp;&nbsp;**if** $\deg r \le \deg s$ **then**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$r, s \leftarrow s, r$
7: &nbsp;&nbsp;&nbsp;&nbsp;**return** $\neg(\deg r = 1$ and $r \ne s)$

**Proposition 2.** *Given input $g, h \in K[X]$ such that $\deg g \ge \deg h$, Algorithm 1 terminates using at most*

$$ \frac{1}{2}(\deg g + \deg h + 2)(\deg g + \deg h + 3) - 6 $$

multiplications in $K$.

*Proof.* Line 1 incurs at most $\deg g + \deg h + 2$ multiplications in $K$. Setting $r_0 := r, s_0 := s$, we define this to be loop 0. For $i \ge 1$, we denote by $r_i, s_i$ (where $\deg s_i \ge \deg r_i$) the polynomials in loop $i$ of Lines 2-6. Using this notation, we move to Line 7 when $\deg r_i \le 1$ or $r_i = s_i$. Now, in loop $i \ge 1$ we replace $r_i$ by $r_i - X^{\deg r_i - \deg s_i} s_i$, meaning $\deg r_{i-1} - \deg r_i \ge 1$, and compute $r_i \cdot LC(s_i)$ and $s_i \cdot LC(r_i)$. This requires $\deg r_i + \deg s_i + 2$ multiplications in $K$. In the worst case, we have $\deg r_{i-1} - \deg r_i = 1$ for $i \ge 1$, where the number of multiplications will decrease by exactly 1 after each loop. In the final loop we have $\deg r_i, \deg s_i = 1$, so we compute 4 multiplications in $K$. In summary, in the worst case we begin with $\deg g + \deg h + 2$ multiplications, decreasing by 1 until we get to 4. Therefore, the total number of multiplications is at most $\sum_{n=4}^{\deg g + \deg h + 2} n$, which is the bound above. $\square$

In summary, Proposition 1 shows that detecting subfield roots of $f \in \mathbb{F}_{q^d}[X]$ amounts to computing the GCD of $d$ related polynomials in $\mathbb{F}_q[X]$. We showed that computing this GCD is simpler when $d=2$. Proposition 2 gives an upper bound on the number of $\mathbb{F}_q$ multiplications required to compute such a GCD in $\mathbb{F}_q[X]$. In the next section we use these tools to build a faster algorithm for finding subfield nodes in supersingular isogeny graphs.
## 5 SuperSolver: optimised subfield searching with fast subfield root detection in $\mathcal{X}(\overline{\mathbb{F}}_p, \ell)$

SuperSolver is an algorithm which, given two $j$-invariants in $S_{p^2}$ corresponding to two supersingular curves $E_1/\overline{\mathbb{F}}_{p^2}$ and $E_2/\overline{\mathbb{F}}_{p^2}$, will, on average, solve the supersingular isogeny problem with lower concrete complexity than the traditional Delfs-Galbraith Solver algorithm described in Section 3. As in the Delfs-Galbraith algorithm, SuperSolver takes non-backtracking walks in $\mathcal{X}(\overline{\mathbb{F}}_p, 2)$ until they hit a $j$-invariant in $\overline{\mathbb{F}}_p$. However, at each step of the random walk, SuperSolver also inspects $\mathcal{X}(\overline{\mathbb{F}}_p, \ell)$, for carefully chosen $\ell > 2$, to efficiently detect whether $j$ has any $\ell$-isogenous neighbours in $\overline{\mathbb{F}}_p$. Traditionally, inspecting $\mathcal{X}(\overline{\mathbb{F}}_p, \ell)$ for a subfield neighbour requires fully factoring a degree-$N_\ell$ polynomial and determining whether any of the roots lie in $\overline{\mathbb{F}}_p$. Performing this for each $\ell$ would require $O(\ell^3 + 2\ell^2 \log p)$ operations in $\overline{\mathbb{F}}_{p^2}$ using the modified Cantor-Zassenhaus algorithm (see Section 2), which is prohibitively costly. Following the results from Section 4, however, SuperSolver conducts the inspection of $\mathcal{X}(\overline{\mathbb{F}}_p, \ell)$ with $O(\ell^2)$ multiplications in $\overline{\mathbb{F}}_p$. We make this count precise later in this section. Crucially, the number of $\overline{\mathbb{F}}_p$ operations is no longer dependent on the size of $p$, and this means that as $p$ grows large, the set of $\ell$'s that are optimal to use also grows, and the more profitable (relatively speaking) SuperSolver becomes. We reiterate that, although both Solver and SuperSolver return the full isogeny between $E_1/\overline{\mathbb{F}}_{p^2}$ and $E_2/\overline{\mathbb{F}}_{p^2}$, our discussion focusses on the bottleneck problem of finding an isogeny from $E_1/\overline{\mathbb{F}}_{p^2}$ (resp. $E_2/\overline{\mathbb{F}}_{p^2}$) to $E'_1/\overline{\mathbb{F}}_p$ (resp. $E_2/\overline{\mathbb{F}}_p$). If, at some node $j$, we detect an $\ell$-isogenous neighbour in $\overline{\mathbb{F}}_p$, SuperSolver will then factorise the degree-$N_\ell$ polynomial $\Phi_{\ell,p}(X, j)$ to determine the subfield $j$-invariant. We view this as a post-computation step, since we are only interested in the concrete complexity of the average step taken in the walk (which we assume does not find a subfield node). Note that the paths between $E_1/\overline{\mathbb{F}}_{p^2}$ and $E_2/\overline{\mathbb{F}}_{p^2}$ returned by both Solver and SuperSolver both look the same: in general, both start and finish with a chain of 2-isogenies that is connected in the middle by a chain of different prime-degree isogenies. The main difference, as the results in Section 7 illustrate, is that 2-isogeny chains at each end are much shorter. Recall that in the original Delfs-Galbraith algorithm, each step consists of finding the roots of a quadratic equation in $\overline{\mathbb{F}}_{p^2}[X]$, which reveals two neighbouring nodes in $\mathcal{X}(\overline{\mathbb{F}}_p, 2)$. In SuperSolver, after forming a list of carefully chosen $\ell > 2$, each step will also include the rapid inspection of $\mathcal{X}(\overline{\mathbb{F}}_p, \ell)$ for every $\ell$ in this list. Though the inspection of the neighbours in $\mathcal{X}(\overline{\mathbb{F}}_p, \ell)$ increases the total number of $\overline{\mathbb{F}}_p$ multiplications at each step, more nodes are checked. We first describe the process of taking a step in SuperSolver, and then move to describing how to choose the list of $\ell > 2$ in order to minimise the number of $\overline{\mathbb{F}}_p$ multiplications per node inspected.

**Remark 4 (Odd $\ell$ only).** With the exception of the leaf nodes in the last level of the binary tree, it is redundant to perform rapid node inspections in $\mathcal{X}(\overline{\mathbb{F}}_p, 2\ell)$ if rapid inspections in $\mathcal{X}(\overline{\mathbb{F}}_p, \ell)$ are also part of the routine, since the latter inspections will detect (or exclude) subfield nodes at the next level of the walk down the tree. We therefore find it optimal to only include odd $\ell_i$ in the lists
constructed at the end of this section. Note that there is no redundancy in including odd composite $l_i$'s in our lists, even if they have proper divisors that are also in the list.

**Rapid inspection of the $l$-isogenous neighbours.** Here we describe Algorithm 2: NeighbourInFp. On input of $l, j \in \mathbb{F}_{p^2}$ and $p$, it outputs true if $j$ is $l$-isogenous to a $j' \in \mathbb{F}_p$, and false otherwise. Recall from Equation (1) that the degree of $\Phi_{l,p}$ in $X$ and $Y$ is $N_l$. The first subroutine of NeighbourInFp is EvalModPoly(j, l, p): it evaluates $\Phi_{l,p}(X, Y)$ at $Y = j$ by computing $j^2, \dots, j^{N_l}$, and then multiplying these by the corresponding coefficients of $\Phi_{l,p}$, returning the coefficients $a_{N_l}, \dots, a_0$ of $X$ in $\Phi_{l,p}(X, j)$. Note that, since we typically have a list of multiple $l$, i.e., $l_1 < \dots < l_t$, the powers of $j$ (up to $N_{l_t}$) are computed once-and-for-all at every $j$, and recycled among the $l_i < l_t$. We follow Section 4 to detect whether $\Phi_{l,p}(X, j) \in \mathbb{F}_{p^2}[X]$ has a root in $\mathbb{F}_p$. Letting $\beta \in \mathbb{F}_{p^2}$ be such that $\mathbb{F}_{p^2} = \mathbb{F}_p(\beta)$, we first compute the related polynomials

$$g_1 := (1/2) \cdot [\Phi_{l,p}(X, j) + \pi(\Phi_{l,p}(X, j))] \quad \text{and}$$

$$g_2 := (-\beta/2) \cdot [\Phi_{l,p}(X, j) - \pi(\Phi_{l,p}(X, j))],$$

where $\pi \in \text{Gal}(\mathbb{F}_{p^2}/\mathbb{F}_p)$ is the Frobenius endomorphism. By Proposition 1, we have $g_1, g_2 \in \mathbb{F}_p[X]$ and

$$\deg(\gcd(g_1, g_2)) = 1 \implies \Phi_{l,p}(X, j) \text{ has a root in } \mathbb{F}_p.$$

We then complete the inspection of $\mathcal{X}(\mathbb{F}_p, l)$ by using Algorithm 1 to calculate $\gcd(g_1, g_2) \neq 1$, then (for large enough $p$) it is overwhelmingly likely that $\deg(\gcd(g_1, g_2)) = 1$, which is why our implementation uses the degree of the GCD as the criterion for terminating the subfield search. Another possibility is to terminate whenever $\gcd(g_1, g_2)$ is non-constant, and then to inspect the higher degree GCD according to the two possible scenarios discussed in Section 4.

Note that if we have a polynomial $f(X) = a_nX^n + a_{n-1}X^{n-1} + \dots + a_1X + a_0 \in \mathbb{F}_{p^2}[X]$ then

$$\frac{1}{2}[f + \pi(f)] = \text{Re}(a_n)X^n + \text{Re}(a_{n-1})X^{n-1} + \dots + \text{Re}(a_1)X + \text{Re}(a_0) \in \mathbb{F}_p[X],$$ $$\frac{-\beta}{2}[f - \pi(f)] = \text{Im}(a_n)X^n + \text{Im}(a_{n-1})X^{n-1} + \dots + \text{Im}(a_1)X + \text{Im}(a_0) \in \mathbb{F}_p[X],$$

where, for $a+b\beta \in \mathbb{F}_{p^2}$, $\text{Re}(a+b\beta) = a$ and $\text{Im}(a+b\beta) = b$, in analogy with the notation used for complex numbers. As a result, we can obtain $g$ and $h$ directly from $f = \Phi_{l,p}$ by computing

$$g_1 = X^{N_l} + \dots + \text{Re}(a_0), \quad \text{and} \quad g_2 = \text{Im}(a_{N_l-1})X^{N_l-1} + \dots + \text{Im}(a_0).$$

This avoids having to compute any $\mathbb{F}_{p^2}$ multiplications to calculate the related polynomials $g_1, g_2$.
**Algorithm 2 NeighbourInFp()**: Detect whether $j \in \mathbb{F}_{p^2}$ is $\ell$-isogenous to a $j' \in \mathbb{F}_p$

**Input:** $\ell, j, p$
1: $a_{N_\ell}, \dots, a_0 \leftarrow \text{EvalModPolyj}(\ell, j, p)$
2: $g_1 \leftarrow X^{N_\ell} + \dots + \text{Re}(a_0)$
3: $g_2 \leftarrow \text{Im}(a_{N_\ell-1})X^{N_\ell-1} + \dots + \text{Im}(a_0)$
4: **return** $\text{InvFreeGCD}(g_1, g_2)$

**Cost of Inspecting the $\ell$-isogeny Graph.** Evaluating $\Phi_{\ell,p}(X, Y)$ at $Y = j$ with EvalModPolyj requires at most $9N_\ell(N_\ell-1)$ multiplications in $\mathbb{F}_p$, noting that one $\mathbb{F}_{p^2}$ multiplication is equivalent to $3 \mathbb{F}_p$ multiplications. By Proposition 2, we compute $\text{InvFreeGCD}(g_1, g_2)$ with at most $(2N_\ell+1)(N_\ell+1)-6 \mathbb{F}_p$ multiplications. Therefore, for a fixed $\ell$, the cost of inspecting $\mathcal{X}(\mathbb{F}_p, \ell)$ is

$$ \begin{aligned} \text{cost}_{\ell} &= \frac{1}{N_{\ell}} [\#\mathbb{F}_{p} \text{ multiplications needed to inspect } \ell\text{-isogenous neighbours}] \\ &\le \frac{1}{N_{\ell}} [11N_{\ell}^2 - 6N_{\ell} - 5], \end{aligned} $$

which depends only on $\ell$. This means that, for each $\ell$, $\text{cost}_\ell$ can be computed once for all primes. In Table 3 we present the $\ell$ with the lowest cost, ordering them by increasing $\text{cost}_\ell$ from left to right.

<table><caption>**Table 3.** The cost of inspecting $\ell$-isogenous neighbours, $\text{cost}_\ell$, for $\ell$ ordered by increasing cost from left to right.</caption><tbody><tr><th scope="row">ℓ</th><td>3</td><td>5</td><td>7</td><td>9</td><td>11</td><td>13</td><td>17</td><td>19</td><td>15</td><td>23</td></tr><tr><th scope="row">N<sub>ℓ</sub></th><td>4</td><td>6</td><td>8</td><td>12</td><td>12</td><td>14</td><td>18</td><td>20</td><td>24</td><td>24</td></tr><tr><th scope="row">F<sub>p</sub> mul per node</th><td>16.3</td><td>24.5</td><td>32.6</td><td>48.8</td><td>48.8</td><td>56.8</td><td>72.8</td><td>80.9</td><td>96.9</td><td>96.9</td></tr></tbody><tbody><tr><th scope="row">ℓ</th><td>25</td><td>29</td><td>21</td><td>31</td><td>27</td><td>37</td><td>41</td><td>43</td><td>33</td><td>35</td></tr><tr><th scope="row">N<sub>ℓ</sub></th><td>30</td><td>30</td><td>32</td><td>32</td><td>36</td><td>38</td><td>42</td><td>44</td><td>48</td><td>48</td></tr><tr><th scope="row">F<sub>p</sub> mul per node</th><td>120.9</td><td>120.9</td><td>128.9</td><td>128.9</td><td>144.9</td><td>152.9</td><td>168.9</td><td>177.0</td><td>192.9</td><td>192.9</td></tr></tbody></table>

The important takeaway from Table 3 is that the number of $\mathbb{F}_p$ multiplications incurred by our algorithm does not grow with $p$. This count is fixed and depends only on $\ell$. Looking back at the root solving algorithms in Section 2, we see a stark difference in expected performance. Those algorithms have many constants hidden by the big-O, have a leading $\ell^3$ term (compared to our $\ell^2$ term), and, importantly, the number of field operations they incur grows as the field grows due to their implicit dependency on $\log p$. Moreover, as mentioned in Section 2, the complexities cited are for *probabilistic* root finding algorithms. Their deterministic variants have even worse complexities [26, §20.6].
**Choosing the $\ell_i$ to minimise the cost of a step.** We consider the cost of each step in SuperSolver, which we denote by the ratio

$$ \text{cost} = \frac{\text{total \# of } \mathbb{F}_p \text{ multiplications}}{\text{total \# of nodes revealed}}. \qquad (3) $$

The aim of this section is to describe how to construct a list of $\ell_i$ that minimises the cost. Recall from in Table 3 that the $\ell$'s that give the cheapest cost per node inspected are (from left to right)

$$ [3, 5, 7, 11, 13, 9, 17, 19 \dots]. \qquad (4) $$

We will use $L_b$ to denote each list of $\ell_i$ and $\text{cost}_{L_b}$ to denote the corresponding cost, where the bit representation of $b$ specifies the set of $\ell$'s from Equation (4); the least significant bit of $b$ determines if 3 is included, the second least significant bit of $b$ determines if 5 is included, and so on. For example, $L_0 = \{\}$, $L_2 = \{5\}$, and and $L_{11} = \{3, 5, 7\}$. Each step will always include revealing 2 neighbours in $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$, therefore for a node $j$ we have for each step:

total \# of $\mathbb{F}_p$ mulfs. $\ge \#\mathbb{F}_p$ mulfs. needed to find roots of $\Phi_{2,p}(X, j)$;

total \# of nodes revealed $\ge 2$.

Here, equality holds only when we take the list to be $L_0$, which corresponds to the original Delfs-Galbraith algorithm. Minimising the cost in Equation (3) is a non-trivial task. We first restrict the $L_b$ to only contain $\ell$ such that $\text{cost}_\ell < \text{cost}_{L_0}$, otherwise it would be more advantageous to take another step by moving to a neighbouring node in $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$. We emphasise that $\text{cost}_{L_0}$ grows with $p$, whereas $\text{cost}_\ell$ stays fixed. This signifies that the condition on $\ell$ becomes less restrictive as $p$ increases. Suppose that, imposing this condition we get $L_b \subseteq [\ell_1, \dots, \ell_n]$. We then exhaust all $b < 2^n$, corresponding to subsets of $[\ell_1, \dots, \ell_n]$, to determine the $L_b$ that minimise Equation (3). It is important to note that, as this optimisation depends only on the prime $p$, $L_b$ can be determined in the precomputation.

# 6 A worked example

We now use a worked example to illustrate how the Solver and SuperSolver programs solve the supersingular isogeny problem, and to highlight the differences between them. Our SuperSolver suite is written in Sage/Python and a boolean variable supersolver specifies whether Solver or SuperSolver is used. For a prime $p$, and two supersingular $j$-invariants $j_1$ and $j_2$ defined over $\mathbb{F}_{p^2} = \mathbb{F}_p(\beta)$, Solver runs by entering

Solver(p, j10, j11, j20, j21, false)

and SuperSolver runs by calling

Solver(p, j10, j11, j20, j21, true),

where $j10 = \text{Re}(j_1)$, $j11 = \text{Im}(j_1)$ and similar for $j20, j21$.
We picked

$$p = 2^{20} - 3,$$

the smallest of the primes from Table 4 (of Section 7), and generated two pseudo-random<fnref n="8" /> $j$-invariants in $S_{p^2} \setminus S_p$:

$$j_1 = 129007\beta + 818380 \quad \text{and} \quad j_2 = 97589\beta + 660383.$$

**Preprocessing.** The preprocessing phase of both programs starts by constructing the extension field $\mathbb{F}_{p^2} = \mathbb{F}_p(\beta)$, where $\beta^2$ is the first non-square in the sequence $-1, -2, 2, -3, 3, \dots$. It then computes a list of constants for the Tonelli-Shanks subroutine, most notably the exponent $(p - 2^e - 1)/2^{e+1}$, where $e$ is the maximum integer such that $2^e \mid (p-1)$. This exponent is Scott's 'progenitor' [25, p. 3], which essentially determines the complexity of $\mathbb{F}_p$ square roots, and therefore of $\mathbb{F}_{p^2}$ square roots. As a result, it determines the cost of taking a step in $\mathcal{X}(\mathbb{F}_p, 2)$ - see Section 3. The preprocessing phase then computes a set of integers $\ell \ge 3$ (according to the optimisations in Section 5 and the relevant heuristics in [13]), fetches the associated files (originally from Sutherland's database [28]) containing $\Phi_\ell(X, Y) \in \mathbb{Z}[X, Y]$ and reduces all of the coefficients to store a set of new, more compact files containing elements of $\mathbb{F}_p$ that define each of the $\Phi_{\ell,p}(X, Y) \in \mathbb{F}_p[X, Y]$. Note that this is done for both Solver and SuperSolver, since both of these programs use the original Delfs-Galbraith subfield path algorithm [13, Algorithm 1] after the searches for subfield nodes is complete. It is important to note, especially in the cryptanalytic context, that all of these preprocessing steps only depend on $p$ and can therefore be done without knowledge of $j_1$ and $j_2$.

**Solver.** The optimised walk in $\mathcal{X}(\mathbb{F}_p, 2)$ proceeds exactly as described in Section 3, i.e., using the depth first search through the binary trees rooted at $j_1$ and $j_2$, until both searches find the subfield nodes $j'_1 \in \mathbb{F}_p$ and $j'_2 \in \mathbb{F}_p$. In the case of our example, paths were found to $j'_1 = 760776$ and $j'_2 = 35387$, depicted in Figure 2 and Figure 3. They correspond to $\phi_1: E_1 \to E'_1$ and $\phi_2: E_2 \to E'_2$, where $j(E_1) = j_1$, $j(E'_1) = j'_1$, $j(E_2) = j_2$, and $j(E'_2) = j'_2$.

Solver then computes a connecting path between the subfield nodes following Delfs-Galbraith [13, Algorithm 1]. This is depicted in Figure 4. Solver simply reverses the steps in $\phi_2$ to obtain its dual, $\hat{\phi}_2$, and outputs the full path as $\phi: E_1 \to E_2$ as $\phi = \hat{\phi}_2 \circ \phi' \circ \phi_1$.

**SuperSolver.** With $p = 2^{20} - 3$, the preprocessing phase determined that SuperSolver is optimal with $L_3 = \{3, 5\}$ (see also Table 4 in the next section). Before departing the starting node $j_1 = 129007\beta + 818380$, SuperSolver performs the rapid inspection of its 3- and 5-isogenous neighbours as described in Section 5. It then takes steps in $\mathcal{X}(\mathbb{F}_p, 2)$ as in Section 3, but at each new node it performs

<footnote n="8">⁸We do this by taking long walks in $\mathcal{X}(\mathbb{F}_p, 3)$ away from a known subfield curve.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** A walk through $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$ for $p = 2^{20} - 3$ during Solver. The walk starts at $j_1 = 129007\beta + 818380 \in S_{p^2}$ and finds the subfield node $j'_1 = 760776 \in S_p$ after 21 steps.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. A walk through $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$ for $p = 2^{20} - 3$ during Solver. The walk starts at $j_2 = 97589\beta + 660383 \in S_{p^2}$ and finds the subfield node $j'_2 = 35387 \in S_p$ after 21 steps.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. A path connecting two subfield $j$-invariants by taking steps in $\mathcal{X}(\bar{\mathbb{F}}_p, \ell)$ with $\ell \in \{17, 29, 31, 37\}$. The walk starts at $j'_1 = 760776 \in S_p$ and connects to $j'_2 = 35387 \in S_p$ after 8 steps.</figcaption></figure>
the rapid inspection of the 3- and 5-isogenous neighbours. In our example, both walks found a subfield node after 2 steps in $\mathcal{X}(\mathbb{F}_p, 2)$. The walk from $j_1$ found a 3-isogenous neighbour and the walk from $j_2$ found a 5-isogenous neighbour. The final step that finds $\phi'$ is implemented in SuperSolver exactly as it was for Solver. The three isogenies $\phi_1, \phi_2$, and $\phi'$, comprising the full isogeny $\phi = \phi_2 \circ \phi' \circ \phi_1$, are depicted in Figure 5.

$$ \phi_1: j_1 \xrightarrow{2} 219247\beta + 863507 \xrightarrow{2} 489342\beta + 132142 \xrightarrow{3} j'_1 = 35387 $$

$$ \phi_2: j_2 \xrightarrow{2} 867493\beta + 220256 \xrightarrow{2} 252807\beta + 1011175 \xrightarrow{5} j'_2 = 292917 $$

$$ \phi': j'_1 \xrightarrow{17} 658300 \xrightarrow{29} 343840 \xrightarrow{31} 560315 \xrightarrow{17} j'_2 \xleftarrow{37} 439276 $$

<caption>Fig. 5. The three paths found comprising an isogeny from $E_1$ to $E_2$ as found by SuperSolver.</caption>

To illustrate the core idea in this paper, we focus on the isogeny $\phi_1$ depicted at the top of Figure 5 and walk through the steps of the NeighbourInFp algorithm. Evaluating the third modular polynomial at the intermediate $j$-invariants (Step 1 of Algorithm 2) yields

$$ \Phi_{3,p}(X, 219247\beta + 863507) = X^4 + (212814\beta + 479338)X^3 + (408250\beta + 920025)X^2 + (811739\beta + 93038)X + 942336\beta + 847782; $$

$$ \Phi_{3,p}(X, 489342\beta + 132142) = X^4 + (872004\beta + 13960)X^3 + (1031755\beta + 822066)X^2 + (969683\beta + 747785)X + 813010\beta + 255391. $$

Though the theory tells us that these two polynomials split over $\mathbb{F}_{p^2}[X]$, to the naked eye there is no way to distinguish which (if any) of these polynomials has a root in $\mathbb{F}_p$. In both cases, setting $g_1 = 1/2 \cdot (\Phi_{3,p} + \pi(\Phi_{3,p}))$ (Step 2 of Algorithm 2) and $g_2 = -\beta/2 \cdot (\Phi_{3,p} - \pi(\Phi_{3,p}))$ (Step 3 of Algorithm 2) respectively yields

$$ g_1 = X^4 + 479338X^3 + 920025X^2 + 93038X + 847782; $$

$$ g_2 = 425628X^3 + 816500X^2 + 574905X + 836099, $$

and

$$ g_1 = X^4 + 13960X^3 + 822066X^2 + 747785X + 255391; $$

$$ g_2 = 695435X^3 + 1014937X^2 + 890793X + 577447. $$
In the first case, Step 4 of Algorithm 2 outputs $\gcd(g_1, g_2) = 1$, meaning that $\Phi_{3,p}(X, 219247\beta + 863507)$ has no subfield roots. In the second case, we see $\gcd(g_1, g_2) = X + 1013186$, meaning that $-1013186 = 35387$ is a subfield root. In our example, we note that the total number of steps between $j_1$ and $j_2$ returned by SuperSolver is 10, which is much shorter than the 50 steps taken by Solver. Since the middle subfield path finding algorithm is the same in both routines, there is no guarantee that the total path will always be smaller for SuperSolver. It is worth pointing out, however, that the two outer paths from elements in $S_{p^2} \setminus S_p$ to $S_p$ (i.e., $\phi_1$ and $\phi_2$) returned by SuperSolver will never be longer than those returned by Solver. Indeed, Solver can be viewed as a special case of SuperSolver where the list of $\ell$'s is chosen to be $L_0$. Finally, we note that both Solver and SuperSolver always conclude by checking the correctness of the full path from $j_1$ to $j_2$.

# 7 Implementation results

In this section we present some experimental results highlighting the efficacy of SuperSolver. The experiments focus solely on the search for subfield nodes (i.e., the bottleneck step of Delfs-Galbraith) and come in two flavours: many $j$-invariants over small primes, and one $j$-invariant over a large, cryptographic prime.

**Small primes and many walks.** Table 4 and Table 5 report experiments that were run on the largest primes of the 30 bitlengths from 20 to 49. We started at 5000 pseudo-random<fnref n="9" /> supersingular $j$-invariants in $S_{p^2} \setminus S_p$ for the primes of bitlengths 20-24, at 1000 $j$'s for the primes of bitlengths 25-29, at 500 $j$'s for the primes of bitlengths 30-34, at 100 $j$'s for the primes of bitlengths 35-39, at 50 $j$'s for the primes of bitlengths 40-44, and at 10 $j$'s for the primes of bitlengths 45-49. For every $j$, we ran both Solver and SuperSolver (with the five sets of $\ell$'s that were predicted to perform best during preprocessing) until all walks hit a subfield $j$-invariant. Throughout, we will denote these fast sets of $\ell$'s by $L_b$, as in Section 5. In all cases we counted the exact number of $\mathbb{F}_p$ multiplications, squarings and additions required to find the subfield node. Following our metric in Section 2, Table 5 reports the average number of $\mathbb{F}_p$ multiplications by counting squarings as multiplications, and highlights in bold which of the five predicted sets of $\ell$'s performed best on average.

Table 4 reports the average number of nodes visited in each of the walks, along with $\#S_{p^2}/\#S_p$, the expected number of random elements in $S_{p^2}$ that would need to be sampled to find a subfield element in $S_p$. Here, the primes are small enough that $S_p$ can be computed precisely (see Section 2). For each prime, Table 4 highlights in bold the column that matches up with the least multiplications reported in Table 5. Note that, for SuperSolver, the number of

<footnote n="9">⁹Just as in Section 6, we used long walks in $\mathcal{X}(\mathbb{F}_p, 3)$ away from a known starting curve to achieve uniformity in $S_{p^2}$.</footnote>
nodes visited is the number of nodes that are actually walked onto in $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$, not the number of nodes inspected using our fast subfield detection algorithm. Thus, in general, the lowest average number of nodes visited does not correspond to the lowest average number of multiplications. Indeed, the walks with fewer $\ell$'s spend less compute time inspecting $\ell$-isogenous neighbours and therefore move onto new nodes faster, but do not cover as much of the supersingular set during the fast inspection.

<table><caption>**Table 4.** The average number of nodes visited in the search for subfield $j$-invariants in Solver and SuperSolver. Further explanation in text.</caption><thead><tr><th rowspan="2">Prime p</th><th rowspan="2">p mod 8</th><th rowspan="2">[#S<sub>p</sub><sup>2</sup>/#S<sub>p</sub>]</th><th rowspan="2">F<sub>p</sub>-mults. per step</th><th colspan="2">Solver</th><th colspan="5">SuperSolver</th></tr><tr><th>Fastest L<sub>j</sub>'s [L<sub>(i)</sub> . . . , L<sub>(v)</sub>]</th><th>L<sub>0</sub></th><th>L<sub>(i)</sub></th><th>L<sub>(ii)</sub></th><th>L<sub>(iii)</sub></th><th>L<sub>(iv)</sub></th><th>L<sub>(v)</sub></th></tr></thead><tbody><tr><td>2<sup>20</sup> - 3</td><td>5</td><td>530</td><td>54</td><td>[L<sub>3</sub>, L<sub>1</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>2</sub>]</td><td>812</td><td><b>127</b></td><td>257</td><td>76</td><td>107</td><td>193</td></tr><tr><td>2<sup>21</sup> - 9</td><td>7</td><td>156</td><td>53</td><td>[L<sub>3</sub>, L<sub>1</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>2</sub>]</td><td>459</td><td><b>86</b></td><td>218</td><td>53</td><td>87</td><td>111</td></tr><tr><td>2<sup>22</sup> - 3</td><td>5</td><td>584</td><td>60</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>1</sub>, L<sub>5</sub>, L<sub>6</sub>]</td><td>885</td><td>170</td><td>108</td><td><b>288</b></td><td>146</td><td>145</td></tr><tr><td>2<sup>23</sup> - 15</td><td>1</td><td>583</td><td>71</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>838</td><td><b>172</b></td><td>106</td><td>169</td><td>121</td><td>430</td></tr><tr><td>2<sup>24</sup> - 3</td><td>5</td><td>1277</td><td>64</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>1</sub>, L<sub>6</sub>]</td><td>1897</td><td><b>318</b></td><td>209</td><td>311</td><td>618</td><td>273</td></tr><tr><td>2<sup>25</sup> - 39</td><td>1</td><td>1231</td><td>71</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>1</sub>, L<sub>6</sub>]</td><td>1873</td><td><b>360</b></td><td>223</td><td>359</td><td>933</td><td>259</td></tr><tr><td>2<sup>26</sup> - 5</td><td>3</td><td>732</td><td>62</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>1</sub>, L<sub>5</sub>, L<sub>6</sub>]</td><td>1362</td><td>352</td><td><b>194</b></td><td>691</td><td>271</td><td>233</td></tr><tr><td>2<sup>27</sup> - 39</td><td>1</td><td>2348</td><td>73</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>3455</td><td>917</td><td><b>438</b></td><td>579</td><td>497</td><td>1766</td></tr><tr><td>2<sup>28</sup> - 57</td><td>7</td><td>2965</td><td>64</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>1</sub>, L<sub>5</sub>, L<sub>6</sub>]</td><td>9748</td><td>1788</td><td>1022</td><td><b>3065</b></td><td>1314</td><td>1306</td></tr><tr><td>2<sup>29</sup> - 3</td><td>5</td><td>2953</td><td>74</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>4384</td><td>1053</td><td><b>526</b></td><td>712</td><td>603</td><td>2161</td></tr><tr><td>2<sup>30</sup> - 35</td><td>5</td><td>3965</td><td>75</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>5555</td><td>1443</td><td>749</td><td><b>961</b></td><td>849</td><td>2825</td></tr><tr><td>2<sup>31</sup> - 1</td><td>7</td><td>9009</td><td>75</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>27103</td><td>4501</td><td><b>2602</b></td><td>3755</td><td>3136</td><td>8794</td></tr><tr><td>2<sup>32</sup> - 5</td><td>3</td><td>5142</td><td>75</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>10149</td><td>2520</td><td><b>1445</b></td><td>2108</td><td>1702</td><td>5335</td></tr><tr><td>2<sup>33</sup> - 9</td><td>7</td><td>6638</td><td>77</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>20387</td><td><b>3832</b></td><td>2342</td><td>3756</td><td>2676</td><td>10562</td></tr><tr><td>2<sup>34</sup> - 41</td><td>7</td><td>10526</td><td>78</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>32640</td><td>6443</td><td>3790</td><td>6094</td><td><b>4531</b></td><td>16320</td></tr><tr><td>2<sup>35</sup> - 31</td><td>1</td><td>117571</td><td>99</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>150101</td><td><b>14893</b></td><td>27873</td><td>23076</td><td>20921</td><td>9850</td></tr><tr><td>2<sup>36</sup> - 5</td><td>3</td><td>29040</td><td>83</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>63384</td><td>15929</td><td><b>9127</b></td><td>11974</td><td>10807</td><td>5249</td></tr><tr><td>2<sup>37</sup> - 25</td><td>7</td><td>70328</td><td>84</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>218775</td><td><b>26241</b></td><td>16098</td><td>29228</td><td>24153</td><td>10405</td></tr><tr><td>2<sup>38</sup> - 45</td><td>3</td><td>100268</td><td>86</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>217145</td><td>43595</td><td>21343</td><td><b>27187</b></td><td>26982</td><td>14897</td></tr><tr><td>2<sup>39</sup> - 7</td><td>1</td><td>174817</td><td>96</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>230235</td><td>28802</td><td>48488</td><td><b>36770</b></td><td>38318</td><td>19677</td></tr><tr><td>2<sup>40</sup> - 87</td><td>1</td><td>266662</td><td>95</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>394908</td><td>49855</td><td><b>80764</b></td><td>66646</td><td>56901</td><td>28016</td></tr><tr><td>2<sup>41</sup> - 21</td><td>3</td><td>205227</td><td>92</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>448887</td><td>52656</td><td>105639</td><td>69940</td><td>62212</td><td><b>27395</b></td></tr><tr><td>2<sup>42</sup> - 11</td><td>5</td><td>557046</td><td>99</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>720206</td><td>93920</td><td>189498</td><td>147651</td><td><b>102116</b></td><td>64309</td></tr><tr><td>2<sup>43</sup> - 57</td><td>7</td><td>198777</td><td>95</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>705224</td><td><b>69021</b></td><td>153095</td><td>95778</td><td>81922</td><td>44112</td></tr><tr><td>2<sup>44</sup> - 17</td><td>7</td><td>307870</td><td>98</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>808057</td><td>131220</td><td>285136</td><td><b>145263</b></td><td>142750</td><td>72964</td></tr><tr><td>2<sup>45</sup> - 55</td><td>1</td><td>3120225</td><td>108</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>2298828</td><td>301730</td><td><b>410169</b></td><td>579449</td><td>404520</td><td>226542</td></tr><tr><td>2<sup>46</sup> - 21</td><td>3</td><td>2759728</td><td>102</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>9075335</td><td>516826</td><td><b>788898</b></td><td>957832</td><td>730020</td><td>382101</td></tr><tr><td>2<sup>47</sup> - 115</td><td>5</td><td>4234340</td><td>108</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>5182631</td><td>650377</td><td>866413</td><td><b>650377</b></td><td>801837</td><td>781907</td></tr><tr><td>2<sup>48</sup> - 59</td><td>5</td><td>2706129</td><td>111</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>6739857</td><td><b>546014</b></td><td>899553</td><td>756358</td><td>651990</td><td>491312</td></tr><tr><td>2<sup>49</sup> - 81</td><td>7</td><td>1239417</td><td>107</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>3582205</td><td>288124</td><td>660449</td><td><b>326050</b></td><td>319641</td><td>252270</td></tr></tbody></table>

The key trend to highlight is that, relatively speaking, SuperSolver gains more advantage over Solver as the primes get larger. This is not as evident for the small primes in Tables 4 and 5 as it is for the larger primes below.

**Remark 5** ($\mathcal{X}(\bar{\mathbb{F}}_p, 2)$ clusters in $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$). An interesting trend to highlight in Table 4 is that the average number of nodes visited in the optimised Delfs-Galbraith walk through $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$ is significantly more than the expected number of elements one would need to select randomly from $S_p^2$ in order to find an element of $S_p$. The reason for this is that components of $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$ cluster together
<table><caption>**Table 5.** The average number of $\mathbb{F}_p$ multiplications used to search for subfield $j$-invariants in Solver and SuperSolver. Further explanation in text.</caption><thead><tr><th rowspan="2">Prime p</th><th rowspan="2">[<sup>#S<sub>p</sub></sup>/<sub>#S<sub>p</sub></sub>]</th><th rowspan="2">F<sub>p</sub>-mults. per step</th><th colspan="2">Solver</th><th colspan="4">SuperSolver</th></tr><tr><th>Fastest L<sub>j</sub>'s [L<sub>(i)</sub> ... L<sub>(v)</sub>]</th><th>L<sub>0</sub></th><th>L<sub>(i)</sub></th><th>L<sub>(ii)</sub></th><th>L<sub>(iii)</sub></th><th>L<sub>(iv)</sub></th><th>L<sub>(v)</sub></th></tr></thead><tbody><tr><td>2<sup>20</sup> - 3</td><td>530</td><td>54</td><td>[L<sub>3</sub>, L<sub>1</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>2</sub>]</td><td>44848</td><td><b>20601</b></td><td>22585</td><td>22235</td><td>23459</td><td>24951</td></tr><tr><td>2<sup>21</sup> - 9</td><td>156</td><td>53</td><td>[L<sub>3</sub>, L<sub>1</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>2</sub>]</td><td>24187</td><td><b>13648</b></td><td>18578</td><td>15453</td><td>18770</td><td>14064</td></tr><tr><td>2<sup>22</sup> - 3</td><td>584</td><td>60</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>1</sub>, L<sub>5</sub>, L<sub>6</sub>]</td><td>52385</td><td>28062</td><td>31962</td><td><b>26410</b></td><td>32555</td><td>38348</td></tr><tr><td>2<sup>23</sup> - 15</td><td>585</td><td>71</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>59691</td><td><b>30508</b></td><td>32883</td><td>39703</td><td>33370</td><td>44556</td></tr><tr><td>2<sup>24</sup> - 3</td><td>1277</td><td>64</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>1</sub>, L<sub>6</sub>]</td><td>112878</td><td><b>53900</b></td><td>62725</td><td>70482</td><td>59206</td><td>73117</td></tr><tr><td>2<sup>25</sup> - 39</td><td>1231</td><td>71</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>1</sub>, L<sub>6</sub>]</td><td>128703</td><td><b>63021</b></td><td>68210</td><td>83333</td><td>94434</td><td>70878</td></tr><tr><td>2<sup>26</sup> - 5</td><td>732</td><td>62</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>1</sub>, L<sub>5</sub>, L<sub>6</sub>]</td><td>85437</td><td>59484</td><td><b>58286</b></td><td>65813</td><td>61261</td><td>62216</td></tr><tr><td>2<sup>27</sup> - 39</td><td>2348</td><td>73</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>251304</td><td>164036</td><td><b>135672</b></td><td>136633</td><td>137780</td><td>185819</td></tr><tr><td>2<sup>28</sup> - 57</td><td>2965</td><td>64</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>1</sub>, L<sub>5</sub>, L<sub>6</sub>]</td><td>631157</td><td>305345</td><td>308003</td><td><b>298049</b></td><td>299314</td><td>351102</td></tr><tr><td>2<sup>29</sup> - 3</td><td>2953</td><td>74</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>326888</td><td>199985</td><td><b>171489</b></td><td>173335</td><td>177986</td><td>235902</td></tr><tr><td>2<sup>30</sup> - 35</td><td>3965</td><td>75</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>412457</td><td>260188</td><td>232753</td><td><b>228089</b></td><td>236360</td><td>301541</td></tr><tr><td>2<sup>31</sup> - 1</td><td>9009</td><td>75</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>1998840</td><td>809040</td><td><b>807306</b></td><td>889210</td><td>871068</td><td>934319</td></tr><tr><td>2<sup>32</sup> - 5</td><td>5142</td><td>75</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>758637</td><td>455571</td><td><b>449889</b></td><td>501335</td><td>474549</td><td>572203</td></tr><tr><td>2<sup>33</sup> - 9</td><td>6638</td><td>77</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>1564701</td><td><b>700390</b></td><td>733705</td><td>900515</td><td>751310</td><td>1153911</td></tr><tr><td>2<sup>34</sup> - 41</td><td>10526</td><td>78</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>1</sub>]</td><td>2537688</td><td>1184024</td><td>1191084</td><td>1467113</td><td><b>1276654</b></td><td>1799292</td></tr><tr><td>2<sup>35</sup> - 31</td><td>117571</td><td>99</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>15272705</td><td><b>5037679</b></td><td>5790782</td><td>6109529</td><td>6396752</td><td>6213090</td></tr><tr><td>2<sup>36</sup> - 5</td><td>29040</td><td>83</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>5244914</td><td>3006618</td><td><b>2913909</b></td><td>2942626</td><td>3099020</td><td>3211580</td></tr><tr><td>2<sup>37</sup> - 25</td><td>70328</td><td>84</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>18322417</td><td><b>4979176</b></td><td>5155517</td><td>7211517</td><td>6950196</td><td>6375918</td></tr><tr><td>2<sup>38</sup> - 45</td><td>100268</td><td>86</td><td>[L<sub>3</sub>, L<sub>7</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>18402937</td><td>8315681</td><td>6856578</td><td><b>6735588</b></td><td>7791309</td><td>9143526</td></tr><tr><td>2<sup>39</sup> - 7</td><td>174817</td><td>96</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>22505327</td><td>9627241</td><td>9879376</td><td><b>9587856</b></td><td>11562406</td><td>12332858</td></tr><tr><td>2<sup>40</sup> - 87</td><td>266662</td><td>95</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>38602102</td><td>16664021</td><td><b>16455546</b></td><td>17377853</td><td>17169885</td><td>17559520</td></tr><tr><td>2<sup>41</sup> - 21</td><td>205227</td><td>92</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>41185437</td><td>17284297</td><td>20890068</td><td>17817383</td><td>18399209</td><td><b>17006068</b></td></tr><tr><td>2<sup>42</sup> - 11</td><td>557046</td><td>99</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>70760036</td><td>31439715</td><td>38704883</td><td>38573868</td><td><b>30864574</b></td><td>40337712</td></tr><tr><td>2<sup>43</sup> - 57</td><td>198777</td><td>95</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>66820000</td><td><b>22863425</b></td><td>30733754</td><td>24686759</td><td>24474388</td><td>27514948</td></tr><tr><td>2<sup>44</sup> - 17</td><td>307870</td><td>98</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>79795521</td><td>43991657</td><td>58381667</td><td><b>38022655</b></td><td>43217829</td><td>45803624</td></tr><tr><td>2<sup>45</sup> - 55</td><td>3120225</td><td>108</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>247697962</td><td>103871110</td><td><b>87674099</b></td><td>156886333</td><td>126109617</td><td>144250917</td></tr><tr><td>2<sup>46</sup> - 21</td><td>2759728</td><td>102</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>923415913</td><td>174816651</td><td><b>163893709</b></td><td>198006728</td><td>205399202</td><td>220786555</td></tr><tr><td>2<sup>47</sup> - 115</td><td>4234340</td><td>108</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>550653552</td><td>222915969</td><td>183895536</td><td><b>175113230</b></td><td>248769348</td><td>272706341</td></tr><tr><td>2<sup>48</sup> - 59</td><td>2706129</td><td>111</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>729589278</td><td><b>188237971</b></td><td>192728454</td><td>205161765</td><td>251091015</td><td>310387211</td></tr><tr><td>2<sup>49</sup> - 81</td><td>1239417</td><td>107</td><td>[L<sub>7</sub>, L<sub>3</sub>, L<sub>5</sub>, L<sub>6</sub>, L<sub>23</sub>]</td><td>385982057</td><td>99186957</td><td>141171527</td><td><b>88278361</b></td><td>99648273</td><td>160633132</td></tr></tbody></table>

in $\mathcal{X}(\mathbb{F}_p, 2)$. Thus, with respect to finding subfield nodes, walks in $\mathcal{X}(\mathbb{F}_p, 2)$ are significantly different from selecting nodes at random from $S_{p^2}$. The types of clusterings in $\mathcal{X}(\mathbb{F}_p, 2)$ depend on the value of $p \bmod 8$ [13, Theorem 2.7], which is why this value is given alongside $p$ in each row. Write $N$ for the ratio between the number of nodes we visited on average (i.e., the bold column) and the number of elements we would expect to draw at random from $S_{p^2}$ before finding one in $S_p$ (i.e., $\#S_{p^2}/\#S_p$). Table 4 shows that (i) when $p \equiv 1 \bmod 4$, we typically see $1 \le N \le 2$; (ii) when $p \equiv 3 \bmod 8$, we typically see $2 \le N \le 3$; and (iii) when $p \equiv 7 \bmod 8$, we often see $N > 3$. For more experimental data illustrating this phenomenon, see [3, §4.3]. In practice, we do not see the $N > 1$ as enough of a reason to incur the significant overhead of walking in $\mathcal{X}(\mathbb{F}_p, \ell)$ for $\ell > 2$ instead. In any case, the method of fast subfield root detection proposed in this paper will work regardless of the $\ell$-isogenies that are used to take steps in a given walk. In fact, if walking in $\mathcal{X}(\mathbb{F}_p, \ell)$ for $\ell > 2$ results in better concrete performance than for $\ell = 2$, the greater cost of taking a step in $\mathcal{X}(\mathbb{F}_p, \ell)$ is likely to increase the size of the set of “fast $\ell$’s” and the relative efficacy of invoking subfield root detection.
**Large primes and optimal node coverage.** Table 6 illustrates the increased efficacy of SuperSolver over Solver as the supersingular isogeny graphs get larger. Recall that we reported some of the results from this table up front in Section 1, namely from the experiments using primes from the isogeny literature. We chose the largest prime below $2^k$ for $k \in \{50, 100, \dots, 800\}$, and started from a pseudo-random $j$-invariant in $S_{p^2} \setminus S_p$ as usual. Since these instances are too large to actually run the full subfield search until it terminates, in each case we ran both Solver and SuperSolver (for the three sets of $\ell$'s that were predicted to perform best during preprocessing) until the number of $\mathbb{F}_p$ multiplications used exceeded $10^8$, and then immediately stopped. The numbers reported in bold in Table 6 are the total number of nodes covered (i.e., both walked onto and inspected) during these walks. For the smallest prime $p = 2^{50} - 27$, SuperSolver covers between 3 and 4 times the number of nodes that Solver does; for the largest prime $p = 2^{800} - 105$, SuperSolver covers between 18 and 19 times the number of nodes. Though primes beyond this size are unlikely to be of cryptographic interest, it is worth pointing out that this trend continues: the larger $p$ grows, the more profitable it becomes to keep adding $\ell$'s in the fast subfield inspection algorithm.

**Storing and accessing the reduced modular polynomials.** The unreduced modular polynomials $\Phi_{\ell}(X, Y) \in \mathbb{Z}[X, Y]$ require a significant amount of storage, but recall that the preprocessing phase immediately reduces all of the coefficients into $\mathbb{F}_p$ to produce $\Phi_{\ell,p}(X, Y) \in \mathbb{F}_p[X, Y]$. This can be done once-and-for-all for a specific prime, and this makes the storage and access of the $\Phi_{\ell,p}(X, Y)$ a non-issue. Storing $\Phi_{\ell,p}(X, Y)$ requires at most $(N_{\ell}^2/2) \cdot \log_2(p)$ bits. For example, the largest $\Phi_{\ell,p}(X, Y)$ for the 250-bit prime above is $\Phi_{13,p}(X, Y)$, which requires the storage of at most $N_{13}^2/2 = 14^2/2 = 98$ elements of $\mathbb{F}_p$, around 3KB. The largest $\Phi_{\ell,p}(X, Y)$ for the 800-bit prime above requires the storage of at most $N_{19}^2/2 = 20^2/2 = 200$ elements of $\mathbb{F}_p$, around 20KB. Any of these would comfortably fit into the L1 cache on a modern CPU.

**Concrete security of the supersingular isogeny problem.** Our SuperSolver suite makes it straightforward to obtain precise estimates on the concrete classical security offered by the general supersingular isogeny problem in $S_{p^2}$, for any prime $p$. Combining a small experiment (like those reported in Table 6) with the expected number of nodes one must cover before reaching a subfield node allows us to obtain accurate counts on the expected number of $\mathbb{F}_p$ multiplications, squarings and additions that must be carried out during a full cryptanalytic attack. It is then a matter of costing these $\mathbb{F}_p$ operations with respect to the appropriate metric, whether that be bit operations, cycle counts, gate counts, or circuit depth.

Take, for example, the 256-bit prime

$$p = 737430436214997974490748205438634569979446953723240325119999999999999999999999$$

underlying SQISign [15] to illustrate how our software can be used to obtain precise security estimates. The precomputation phase of SuperSolver (which takes
a few seconds on input of $p$) reveals that taking an optimised step in $\mathcal{X}(\bar{\mathbb{F}}_p, 2)$ costs 407 multiplications in $\mathbb{F}_p$. Based on this cost, the precomputation further determines that the fastest set of $\ell$'s to proceed with are

$$ \ell \in \{3, 5, 7, 11, 13\}. $$

On average, the combination of this set of $\ell$'s and Algorithm 2 reduces the cost of the subfield search from 407 multiplications in $\mathbb{F}_p$ per node to 58.0 multiplications in $\mathbb{F}_p$ per node (see Table 1). Thus, on average, solving the supersingular isogeny problem costs

$$ 58.0 \times \left( \frac{\#S_{p^2}}{\#S_p} \right) \mathbb{F}_p \text{ multiplications.} $$

Since $p \equiv 7 \mod 12$, we have $\#S_{p^2} = \lfloor p/12 \rfloor + 1$ [27, Theorem V.4.1(c)], and since $p \equiv 7 \mod 8$, $\#S_p$ is exactly the class number of the imaginary quadratic field $\mathbb{Q}(\sqrt{-p})$ [13, Equation 1]. We suppose this class number is $N$, i.e., $\#S_p = N$. Writing $N = 2^k$, where $k$ is correct to 3 decimal places, we would obtain that the average cost of breaking this instance of SQISign is $2^{257.622-k}$ multiplications in $\mathbb{F}_p$.<fnref n="10" />

In Table 7 we give average counts for the cost of breaking the supersingular isogeny problem using SuperSolver for a number of primes underlying either B-SIDH or SQISign.

**Acknowledgements.** Thanks to Sam Frengley, Michael Naehrig, Krijn Reijnders, Benjamin Smith, Greg Zaverucha, and the CRYPTO2022 reviewers for their valuable comments on an earlier version of this paper. We also thank Drew Sutherland for answering our questions about alternative modular functions.

## References

1. G. Adj, D. Cervantes-Vázquez, J. Chi-Domínguez, A. Menezes, and F. Rodríguez-Henríquez. On the cost of computing isogenies between supersingular elliptic curves. In *International Conference on Selected Areas in Cryptography*, pages 322–343. Springer, 2018.

2. G. Adj and F. Rodríguez-Henríquez. Square root computation over even extension fields. *IEEE Transactions on Computers*, 63(11):2829–2841, 2013.

3. S. Arpin, C. Camacho-Navarro, K. Lauter, J. Lim, K. Nelson, T. Scholl, and J. Sotáková. Adventures in supersingular land. *Experimental Mathematics*, pages 1–28, 2021.

4. E. R. Berlekamp. Factoring polynomials over large finite fields. *Mathematics of computation*, 24(111):713–735, 1970.

5. W. Beullens, T. Kleinjung, and F. Vercauteren. CSI-FiSh: Efficient isogeny based signatures through class group computations. In S. D. Galbraith and S. Moriai, editors, *ASIACRYPT 2019*, volume 11921 of *Lecture Notes in Computer Science*, pages 227–247. Springer, 2019.

<footnote n="10">¹⁰For large, cryptographic sized primes $p$, computing class numbers is very computationally expensive. Indeed, a recent class group computation for a 512-bit prime terminated in ≈ 52 core years.</footnote>
6. J. H. Bruinier, K. Ono, and A. V. Sutherland. Class polynomials for nonholomorphic modular functions. *Journal of Number Theory*, 161:204-229, 2016.

7. J. Burdges and L. De Feo. Delay encryption. In A. Canteaut and F. Standaert, editors, *Advances in Cryptology - EUROCRYPT 2021 - 40th Annual International Conference on the Theory and Applications of Cryptographic Techniques, Zagreb, Croatia, October 17-21, 2021, Proceedings, Part I*, volume 12696 of *Lecture Notes in Computer Science*, pages 302-326. Springer, 2021.

8. D. G. Cantor and H. Zassenhaus. A new algorithm for factoring polynomials over finite fields. *Mathematics of Computation*, pages 587-592, 1981.

9. W. Castryck and T. Decru. CSIDH on the surface. In J. Ding and J. Tillich, editors, *Post-Quantum Cryptography - 11th International Conference, PQCrypto 2020, Paris, France, April 15-17, 2020, Proceedings*, volume 12100 of *Lecture Notes in Computer Science*, pages 111-129. Springer, 2020.

10. W. Castryck, T. Decru, and F. Vercauteren. Radical isogenies. In S. Moriai and H. Wang, editors, *Advances in Cryptology - ASIACRYPT 2020 - 26th International Conference on the Theory and Application of Cryptology and Information Security, Daejeon, South Korea, December 7-11, 2020, Proceedings, Part II*, volume 12492 of *Lecture Notes in Computer Science*, pages 493-519. Springer, 2020.

11. C. Costello. B-SIDH: supersingular isogeny Diffie-Hellman using twisted torsion. In *ASIACRYPT*, pages 440-463. Springer, 2020.

12. C. Costello, M. Meyer, and M. Naehrig. Sieving for twin smooth integers with solutions to the Prouhet-Tarry-Escott problem. In *EUROCRYPT*, volume 12696, pages 272-301. Springer, 2021.

13. C. Delfs and S. D. Galbraith. Computing isogenies between supersingular elliptic curves over $\mathbb{F}_p$. *Designs, Codes and Cryptography*, 78(2):425-440, 2016.

14. L. De Feo, D. Jao, and J. Plût. Towards quantum-resistant cryptosystems from supersingular elliptic curve isogenies. *J. Math. Cryptol.*, 8(3):209-247, 2014.

15. L. De Feo, D. Kohel, A. Leroux, C. Petit, and B. Wesolowski. SQISign: compact post-quantum signatures from quaternions and isogenies. In *ASIACRYPT*, pages 64-93. Springer, 2020.

16. S. D. Galbraith, C. Petit, and J. Silva. Identification protocols and signature schemes based on supersingular isogeny problems. *J. Cryptol.*, 33(1):130-175, 2020.

17. M. Hamburg. Fast and compact elliptic-curve cryptography. Cryptology ePrint Archive, Report 2012/309, 2012. https://ia.cr/2012/309.

18. R. A. Horn and C. R. Johnson. *Topics in matrix analysis*. Cambridge university press, 1994.

19. D. Jao, R. Azarderakhsh, M. Campagna, C. Costello, L. De Feo, B. Hess, A. Jalali, B. Koziel, B. LaMacchia, P. Longa, M. Naehrig, J. Renes, V. Soukharev, and D. Urbanik. SIKE: Supersingular Isogeny Key Encapsulation. Manuscript available at sike.org/, 2017.

20. C. Leonardi. *Security Analysis of Isogeny-Based Cryptosystems*. PhD thesis, University of Waterloo, Ontario, Canada, 2020.

21. R. Lidl and H. Niederreiter. *Introduction to finite fields and their applications*. Cambridge University Press, 1994.

22. A. J. Menezes, P. C. van Oorschot, and S. A. Vanstone. *Handbook of applied cryptography*. CRC press, 2018.

23. J-F. Mestre. La méthode des graphes. exemples et applications. In *Proceedings of the international conference on class numbers and fundamental units of algebraic number fields (Katata)*, pages 217-242. Citeseer, 1986.

24. A. K. Pizer. Ramanujan graphs and Hecke operators. *Bulletin of the American Mathematical Society*, 23(1):127-137, 1990.
25. M. Scott. A note on the calculation of some functions in finite fields: Tricks of the trade. *IACR Cryptol. ePrint Arch.*, page 1497, 2020.

26. V. Shoup. *A computational introduction to number theory and algebra*. Cambridge university press, 2009.

27. J. H. Silverman. *The arithmetic of elliptic curves*, volume 106. Springer, 2009.

28. A. V. Sutherland. Modular Polynomials. https://math.mit.edu/~drew/ClassicalModPolys.html. Online; accessed 30 September 2021.

29. A. V. Sutherland. On the evaluation of modular polynomials. *The Open Book Series*, 1(1):531–555, 2013.

30. The Sage Developers. *SageMath, the Sage Mathematics Software System (Version 9.2)*, 2021. https://www.sagemath.org.

31. J. Vélu. Isogénies entre courbes elliptiques. *CR Acad. Sci. Paris Sér. AB*, 273(A238-A241):5, 1971.
<table><caption>**Table 6.** The number of nodes inspected per $10^8$ field multiplications for the largest primes of various bitlengths. The **Solver** column corresponds to optimised Delfs-Galbraith walks in $\mathcal{X}(\mathbb{F}_p, 2)$ – see Section 3. The **SuperSolver** columns correspond to enabling our fast subfield root detection algorithm with the three fastest sets of $\ell$'s (left to right) – see Section 5. Numbers in round brackets are the approximate number of $\mathbb{F}_p$ multiplications per node inspected, as computed during the precomputation phase that determines which sets of $\ell$'s will perform fastest.</caption><thead><tr><th rowspan="2">Prime p</th><th colspan="3">Solver</th><th colspan="3">SuperSolver</th></tr><tr><th>Nodes inspected</th><th>F<sub>p</sub>-mults. per node</th><th>Fastest Sets</th><th>Nodes inspected</th><th>F<sub>p</sub>-mults. per node</th><th>Improv. Factor</th></tr></thead><tbody><tr><td>2<sup>50</sup> - 27</td><td>883,007</td><td>113</td><td>{3,5,7}<br>{3,5}<br>{3,7}</td><td>2,859,221<br>2,736,601<br>2,533,945</td><td>35.0<br>36.5<br>39.4</td><td>2.8 - 3.2x</td></tr><tr><td>2<sup>100</sup> - 15</td><td>443,951</td><td>223</td><td>{3,5,7}<br>{3,5,7,11}<br>{3,5,7,9}</td><td>2,165,681<br>2,121,313<br>1,988,731</td><td>46.0<br>47.1<br>47.1</td><td>4.5 - 4.9x</td></tr><tr><td>2<sup>150</sup> - 3</td><td>317,215</td><td>315</td><td>{3,5,7,9,11}<br>{3,5,7,11}<br>{3,5,7,9}</td><td>1,847,413<br>1,895,201<br>1,776,751</td><td>51.8<br>52.9<br>52.9</td><td>5.6 - 6.0x</td></tr><tr><td>2<sup>200</sup> - 75</td><td>241,987</td><td>415</td><td>{3,5,7,9,11}<br>{3,5,7,9,11,13}<br>{3,5,7,11,13}</td><td>1,700,749<br>1,715,449<br>1,716,767</td><td>56.2<br>56.3<br>58.3</td><td>7.0 - 7.1x</td></tr><tr><td>2<sup>250</sup> - 207</td><td>191,115</td><td>521</td><td>{3,5,7,9,11,13}<br>{3,5,7,9,11}<br>{3,5,7,11,13}</td><td>1,607,145<br>1,561,645<br>1,586,495</td><td>60.1<br>61.1<br>63.0</td><td>8.2 - 8.4x</td></tr><tr><td>2<sup>300</sup> - 153</td><td>164,273</td><td>610</td><td>{3,5,7,9,11,13}<br>{3,5,7,9,11}<br>{3,5,7,9,11,13,17}</td><td>1,531,993<br>1,468,279<br>1,489,991</td><td>63.0<br>65.0<br>65.3</td><td>8.0 - 8.4x</td></tr><tr><td>2<sup>350</sup> - 113</td><td>141,097</td><td>709</td><td>{3,5,7,9,11,13}<br>{3,5,7,9,11,13,17}<br>{3,5,7,9,11}</td><td>1,452,529<br>1,432,345<br>1,372,351</td><td>66.5<br>68.0<br>70.0</td><td>9.7 - 10.3x</td></tr><tr><td>2<sup>400</sup> - 593</td><td>123,649</td><td>809</td><td>{3,5,7,9,11,13}<br>{3,5,7,9,11,13,17}<br>{3,5,7,9,11,13,19}</td><td>1,380,849<br>1,378,991<br>1,339,805</td><td>70.0<br>70.6<br>72.8</td><td>10.8 - 11.2x</td></tr><tr><td>2<sup>450</sup> - 501</td><td>110,407</td><td>906</td><td>{3,5,7,11,13,9,17}<br>{3,5,7,11,13,9}<br>{3,5,7,11,13,9,17,19}</td><td>1,330,891<br>1,317,849<br>1,309,703</td><td>73.2<br>73.3<br>74.8</td><td>11.9 - 12.1x</td></tr><tr><td>2<sup>500</sup> - 863</td><td>97,209</td><td>1032</td><td>{3,5,7,9,11,13,17}<br>{3,5,7,9,11,13,17,19}<br>{3,5,7,9,11,13}</td><td>1,274,503<br>1,266,275<br>1,245,721</td><td>76.4<br>77.3<br>77.5</td><td>12.8 - 13.1x</td></tr><tr><td>2<sup>550</sup> - 5</td><td>90,031</td><td>1111</td><td>{3,5,7,9,11,13,17}<br>{3,5,7,9,11,13,17,19}<br>{3,5,7,9,11,13}</td><td>1,239,501<br>1,238,921<br>1,201,873</td><td>78.6<br>79.0<br>80.3</td><td>13.3 - 13.8x</td></tr><tr><td>2<sup>600</sup> - 95</td><td>81,253</td><td>1230</td><td>{3,5,7,9,11,13,17,19}<br>{3,5,7,9,11,13,17}<br>{3,5,7,9,11,13,19}</td><td>1,200,945<br>1,191,549<br>1,166,297</td><td>81.5<br>81.7<br>83.6</td><td>14.4 - 14.8x</td></tr><tr><td>2<sup>650</sup> - 611</td><td>76,207</td><td>1314</td><td>{3,5,7,9,11,13,17,19}<br>{3,5,7,9,11,13,17}<br>{3,5,7,9,11,13,19}</td><td>1,176,411<br>1,161,061<br>1,137,873</td><td>83.3<br>83.9<br>85.7</td><td>14.9 - 15.4x</td></tr><tr><td>2<sup>700</sup> - 1113</td><td>71,037</td><td>1408</td><td>{3,5,7,9,11,13,17,19}<br>{3,5,7,9,11,13,17}<br>{3,5,7,9,11,13,17,19,23}</td><td>1,148,963<br>1,127,317<br>1,123,125</td><td>85.2<br>86.4<br>87.5</td><td>15.8 - 16.2x</td></tr><tr><td>2<sup>750</sup> - 161</td><td>66,237</td><td>1510</td><td>{3,5,7,9,11,13,17,19}<br>{3,5,7,9,11,13,17,19,23}<br>{3,5,7,9,11,13,17,19,23}</td><td>1,121,045<br>1,093,351<br>1,101,767</td><td>87.4<br>89.0<br>89.3</td><td>16.5 - 16.9x</td></tr><tr><td>2<sup>800</sup> - 105</td><td>62,163</td><td>1610</td><td>{3,5,7,9,11,13,17,19}<br>{3,5,7,9,11,13,17,19,23}<br>{3,5,7,9,11,13,17,19,15}</td><td>1,095,195<br>1,081,825<br>1,008,481</td><td>89.4<br>90.9<br>90.9</td><td>16.2 - 17.6x</td></tr></tbody></table>
<table><caption>**Table 7.** The average number of $\mathbb{F}_p$ multiplications required to solve the supersingular isogeny problem using **SuperSolver**. When $p \equiv 1 \pmod 4$, we assume that $N = 2^k$ is the class number of $\mathbb{Q}(\sqrt{-4p})$, where $k$ correct to 3 decimal places. Otherwise, it is the class number of $\mathbb{Q}(\sqrt{-p})$. As $N$ varies for each prime, we will index $N$ and $k$ by the row in the table, i.e., $N_i = 2^{k_i}$ will be the class number of the $i$-th prime in the table. The number of $\mathbb{F}_p$ multiplications per node using **SuperSolver** is taken from Table 1.</caption><thead><tr><th>Prime p</th><th>p mod 8</th><th>Average number of F<sub>p</sub>-mults. per node</th><th>#S<sub>p2</sub></th><th>#S<sub>p</sub></th><th>Average cost of SuperSolver</th></tr></thead><tbody><tr><td>B-SIDH-p247 [11]</td><td>7</td><td>55.6</td><td>2<sup>242.559</sup></td><td>2<sup>k<sub>1</sub></sup></td><td>2<sup>248.356-k<sub>1</sub></sup></td></tr><tr><td>TwinSmooth-p250 [12]</td><td>1</td><td>59.1</td><td>2<sup>246.220</sup></td><td>2<sup>k<sub>2</sub>-1</sup></td><td>2<sup>252.105-k<sub>2</sub></sup></td></tr><tr><td>SQISign-p256 [15]</td><td>7</td><td>55.7</td><td>2<sup>251.764</sup></td><td>2<sup>k<sub>3</sub></sup></td><td>2<sup>257.564-k<sub>3</sub></sup></td></tr><tr><td>TwinSmooth-p384 [12]</td><td>1</td><td>63.1</td><td>2<sup>379.735</sup></td><td>2<sup>k<sub>4</sub>-1</sup></td><td>2<sup>385.715-k<sub>4</sub></sup></td></tr><tr><td>TwinSmooth-p512 [12]</td><td>5</td><td>69.1</td><td>2<sup>507.896</sup></td><td>2<sup>k<sub>5</sub>-1</sup></td><td>2<sup>514.007-k<sub>5</sub></sup></td></tr></tbody></table>