# PRISM: Simple and Compact Identification and Signatures from Large Prime Degree Isogenies

Andrea Basso, Giacomo Borin, Wouter Castryck, Maria Corte-Real Santos, Riccardo Invernizzi, Antonin Leroux, Luciano Maino, Frederik Vercauteren, Benjamin Wesolowski

► To cite this version:

Andrea Basso, Giacomo Borin, Wouter Castryck, Maria Corte-Real Santos, Riccardo Invernizzi, et al.: PRISM: Simple and Compact Identification and Signatures from Large Prime Degree Isogenies. PKC 2025, May 2025, Røros, Norway. pp.300-332, 10.1007/978-3-031-91826-1_10. hal-05377080v1

HAL Id: hal-05377080

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

Submitted on 21 Nov 2025 (v1), last revised 10 Dec 2025 (v2)

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 Attribution 4.0 International License
PEGASIS: Practical Effective Class Group
Action using 4-Dimensional Isogenies

Pierrick Dartois¹,², Jonathan Komada Eriksen⁵, Tako Boris Fouotsa³, Arthur Herlédan Le Merdy⁴, Riccardo Invernizzi⁵, Damien Robert¹,², Ryan Rueger⁶,⁷, Frederik Vercauteren⁵ and Benjamin Wesolowski⁴

¹ Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, F-33400, Talence, France

² INRIA, IMB, UMR 5251, F-33400, Talence, France

³ EPFL, LASEC, Lausanne, Switzerland

⁴ ENS de Lyon, CNRS, UMPA, UMR 5669, Lyon, France

⁵ KU Leuven, COSIC, Heverlee, Belgium

⁶ IBM Research Europe, Zürich, Switzerland

⁷ Technische Universität München, Germany

**Abstract.** In this paper, we present the first practical algorithm to compute an effective group action of the class group of any imaginary quadratic order $\mathcal{O}$ on a set of supersingular elliptic curves primitively oriented by $\mathcal{O}$. Effective means that we can act with any element of the class group directly, and are not restricted to acting by products of ideals of small norm, as for instance in CSIDH. Such restricted effective group actions often hamper cryptographic constructions, e.g. in signature or MPC protocols.

Our algorithm is a refinement of the Clapoti approach by Page and Robert, and uses 4-dimensional isogenies. As such, it runs in polynomial time, does not require the computation of the structure of the class group, nor expensive lattice reductions, and our techniques allow it to be instantiated with the orientation given by the Frobenius endomorphism. This makes the algorithm practical even at security levels as high as CSIDH-4096. Our implementation in SageMath takes 1.5s to compute a group action at the CSIDH-512 security level, 21s at CSIDH-2048 level and around 2 minutes at the CSIDH-4096 level. This marks the first instantiation of an effective cryptographic group action at such high security levels. For comparison, the recent KLaPoTi approach requires around 200s at the CSIDH-512 level in SageMath and 2s in Rust.

# 1 Introduction

The simplicity and flexibility of the Diffie-Hellman key-agreement protocol [26] has made it ideally suited to create more advanced protocols, such as public-key encryption [29], digital signatures [62], password-authenticated key exchange [36], threshold encryption [25], threshold signatures [35, 34], updatable encryption [10] and many more. These protocols are all based on the computational efficiency and commutativity of group exponentiation, combined with the computational hardness of the discrete logarithm problem.
Due to Shor's algorithm [63], the discrete logarithm problem is not quantum
hard. To construct a post-quantum analogue of the Diffie-Hellman protocol,
we need to remove the underlying group structure. One way to achieve this is
to use a commutative group action instead. The vectorisation and parallelisa-
tion problems, which are analogues of the discrete-logarithm and computational
Diffie-Hellman problems respectively, are assumed to be (superpolynomially, but
sub-exponentially [39]) hard, even for quantum computers. Furthermore, both
problems are equivalently hard for quantum computers [32, 42, 33]. Using group
actions, it is possible to construct post-quantum analogues of the advanced pro-
tocols given above: public-key encryption [43], digital signatures [7], password-
authenticated key-exchange [1], threshold encryption and signatures [24] and
updatable encryption [40].

A group action is said to be *effective* (EGA) if, amongst other conditions, the
action can be computed in polynomial time [2]. That is, if $G$ acts on $X$, we have
a polynomial-time algorithm that computes $gx$ for any $g$ in $G$ and $x$ in $X$. If we
only have a polynomial-time algorithm that can compute $g_i x$ for $\{g_1, \dots, g_n\}$, a
set of (polynomially many) generators of $G$, the action is said to be a *restricted
effective* group action (REGA).

The most prominent instantiation of a REGA is given by the action of the class group of an imaginary quadratic order on a certain set of oriented elliptic curves, with the CSIDH protocol [13] (or its variant CSURF [12]) being the most efficient. In CSIDH, the orientation is given by the Frobenius endomorphism, and the set of oriented curves simply consists of the supersingular curves defined over $\mathbb{F}_p$, up to $\mathbb{F}_p$-isomorphism. Since we can only efficiently evaluate the action of ideals of smooth norm, CSIDH is indeed a REGA and not an EGA. This often results in difficulties in constructing post-quantum variants of protocols based on Diffie-Hellman: many of the previously mentioned examples really require an EGA, while other examples, such as the signature scheme SeaSign [22], resort to some form of rejection sampling to not leak information about the secret (as a workaround for being a REGA).

A three step approach to turn a REGA into an EGA was given by the sig-
nature scheme CSI-FiSh [7]: the first step computes the structure of the class
group (which can be done in quantum polynomial time, but classically takes
subexponential time), the second step computes a (very) short basis of the rela-
tion lattice (which takes exponential time, even quantumly) and the third step
expresses any element as a product of the generators with small exponents, us-
ing the short basis. Step 1 has been done in practice for CSIDH-512, but further
research [6, 9, 51] concludes that using a 512 bit prime for CSIDH does not
achieve NIST level 1 security, with the latter two sources claiming that primes
of more than 2000 bits are required. Since the computation of the class group
is infeasible for such large primes, the primitive SCALLOP [21] and improve-
ments [14, 3] avoided such computation by using a (large prime conductor)
suborder of an order with small discriminant, so that the class group structure
comes for free. This approach however does not fundamentally solve steps 2 and
3 (see https://yx7.cc/blah/2023-04-14.html for an informal, but detailed
discussion of this issue), and the highest level achieved in [3] is equivalent to CSIDH-1536, for which a single group action takes almost 12 minutes using a C++ implementation. This approach is therefore completely hopeless to instan- tiate practical effective group actions at higher security levels (even using a quantum computer).

**Our contributions.** We describe a practical algorithm to compute the group action of any element in the class group of an imaginary quadratic order $\mathcal{O}$ on a set of supersingular elliptic curves over $\overline{\mathbb{F}}_p$ primitively oriented by $\mathcal{O}$. Our algo- rithm works for all orders of discriminant $\Delta_{\mathcal{O}}$ all the way up to $|\Delta_{\mathcal{O}}| \approx p$ and in particular, applies in the CSIDH/CSURF setting, where the orientation is given by Frobenius, and we have $\Delta_{\mathcal{O}} = -p$ or $\Delta_{\mathcal{O}} = -4p$. As such, we obtain the first EGA at security levels equivalent to CSIDH-2048 and CSIDH-4096. We imple- ment the algorithm in SageMath, and show that it achieves a highly practical runtime, which scales well to higher security levels. Not only does our algorithm allow us to easily instantiate security levels far beyond what previous EGA instantiations were able to; the timings of our proof-of-concept SageMath imple- mentation also outperform all other EGA instantiations (even compared to their highly optimized C++/Rust implementations) at corresponding security levels. Our implementation is publicly available at https://github.com/pegasis4d.

**Technical overview.** Our work is based on the Clapoti framework [49], and proceeds by embedding the isogeny we wish to evaluate in an isogeny in dimen- sion 4. In order to evaluate the action given by the ideal $\mathfrak{a} \subseteq \mathcal{O}$ on an elliptic curve $E$ using only one 4-dimensional $2^e$-isogeny, the Clapoti framework requires finding two other integral ideals $\mathfrak{b}, \mathfrak{c}$ equivalent to $\mathfrak{a}$ such that

$$uN(\mathfrak{b}) + vN(\mathfrak{c}) = 2^e, \qquad (1)$$

with the restriction that $u, v$ are integers that can be written as the sum of two squares (in particular, they must be easy to factorise), and where $2^e < p$ (or more accurately, the full torsion group $E[2^e]$ should be defined over $\overline{\mathbb{F}}_{p^{2k}}$ for a small $k$).

Several issues arise when trying to solve the above equation for $|\Delta_{\mathcal{O}}| \approx p$: even without the restriction on $u$ and $v$ being sums of squares, sometimes the equation simply has no solution. Indeed, since the smallest ideals equivalent to $\mathfrak{a}$ are already expected to have norm around $\sqrt{p}$, and a solution is guaranteed to exist only when $N(\mathfrak{b})N(\mathfrak{c}) < 2^e$, this contradicts the requirement $2^e < p$. As such, imposing the extra condition for both $u$ and $v$ to be sums of squares further lowers the probability of finding a solution.

The main insight is that we can derive a related equation which is much easier to solve by exploiting the fact we can always efficiently compute small degree isogenies. We use this insight in two ways: first, we write $\mathfrak{b} = \mathfrak{b}_k\mathfrak{b}_e$ and $\mathfrak{c} = \mathfrak{c}_k\mathfrak{c}_e$ where $\mathfrak{b}_e$ and $\mathfrak{c}_e$ are ideals of smooth norm and compute the isogenies $\varphi_{\mathfrak{b}_e} : E \to E_1$ and $\varphi_{\mathfrak{c}_e} : E \to E_2$, simply using Elkies' algorithm [30] or Vélu's
algorithm [64]. We can then apply the Clapoti framework to $E_1 \times E_2$ instead of $E \times E$, resulting in the equation

$$uN(\mathfrak{b}_k) + vN(\mathfrak{c}_k) = 2^e,$$

where now the product $N(\mathfrak{b}_k)N(\mathfrak{c}_k)$ is typically much smaller than $2^e$ since we have removed the norms of the ideals $\mathfrak{b}_e$ and $\mathfrak{c}_e$.

Second, we use the same insight to relax the conditions on $u$ and $v$: instead of requiring that each can be written as a sum of squares, we again allow to take out smooth factors $g_u$ and $g_v$ such that $u = g_u(x_u^2+y_u^2)$ and $v = g_v(x_v^2+y_v^2)$. Note that integers containing prime factors (to an odd power) which are 3 mod 4 can never be written as a sum of squares, so allowing factors like 3, 7, 11, 19, ... to be taken out of $u$ and $v$ drastically increases the probability that both $u/g_u$ and $v/g_v$ are sums of squares. Applying the above, allows us to solve the (relaxed) norm equation efficiently, even at the highest security levels.

With $u$ and $v$ of this form, we are able to efficiently construct a 2-dimensional $u$-isogeny $\Phi_u : E_1^2 \to E_u^2$ and a 2-dimensional $v$-isogeny $\Phi_v : E_2^2 \to E_v^2$. We then compute a 4-dimensional $2^e$-isogeny $F : E_u^2 \times E_v^2 \to E_a^2 \times E'^2$ that “embeds” $\Phi_u$, $\Phi_v$, $\phi_{\mathfrak{b}_k} : E_1 \to E_a$ and $\phi_{\mathfrak{c}_k} : E_2 \to E_a$, where $E_a := \mathfrak{a} \cdot E$. Using level-2 theta coordinates as in [18], the computation of $F$ is relatively fast. We finally extract the resulting curve $E_a$ from the codomain of $F$. For very high security levels, we remark that it is helpful to relax the conditions on $u, v$ and allow them to be arbitrary positive integers, at the cost of using more dimension-4 isogenies (see Sections 4.2 and 6.4).

**Comparison with related work.** We briefly compare our algorithm with 2 related works: the original Clapoti framework [49] and the recent KLaPoTi algorithm [50] that relies on KLPT.

**Clapoti:** As described above our work is based on the Clapoti framework [49], which presents the first polynomial time algorithm to act with a random ideal $\mathfrak{a} \subseteq \mathcal{O}$ by finding two integral ideals $\mathfrak{b}, \mathfrak{c}$ equivalent to $\mathfrak{a}$ such that

$$uN(\mathfrak{b}) + vN(\mathfrak{c}) = M. \qquad (2)$$

For a provably polynomial time algorithm, $M$ has to be chosen powersmooth, but in order to achieve a practical algorithm, one really is limited to $M = 2^e$ where the $2^e$-torsion is easily accessible, in particular, $2^e < p$. However, the expected number of solutions is $2^e/N(\mathfrak{b})N(\mathfrak{c})$ where $N(\mathfrak{b}), N(\mathfrak{c})$ are roughly $\sqrt{p}$ and thus if $p+1 = f2^e$ even for small $f$, this equation simply will not have a solution. To illustrate this: for the 4000-bit parameter set, directly applying Clapoti will fail in 97% of the cases when $M = 2^e < p$, and this is without the restriction on $u$ and $v$ being sums of squares. If one adds this restriction on $u$ and $v$, the probability to find a solution is close to 0 since the probability that a pair $(u, v)$ each of size $\sqrt{p}$ is a sum of squares is proportional to $1/\log(p)$. Furthermore, note that it is impossible to test all pairs $(u, v)$ since this would involve factorisation, even further lowering the probability by another 2 to 3
orders of magnitude depending on the size of *p*. In conclusion: running Clapoti with $M = 2^e < p$ will almost always fail, and thus is not a practically efficient algorithm.

**KLaPoti:** A recent paper by Panny, Petit and Stopar [50] has also instantiated a primitive KLaPoTi based on Clapoti. KLaPoTi differs from PEGASIS in several aspects. The main ones are the way Equation (1) is solved, the dimension of the isogeny used in the evaluation of the group action and the size of the base prime *p*. KLaPoTi relies on the KLPT algorithm [38] to solve Equation (1), which now requires *u*, *v* to be perfect squares. The fact that *u* and *v* are perfect squares allows to perform the group action computation using dimension-2 isogenies only. Nevertheless, the relatively large size of the solution obtained by this KLPT approach requires $|\Delta|^3 < 2^e$ where $\Delta$ is the discriminant of the order $\mathcal{O}$. This implies that KLaPoTi uses a base prime *p* whose bitsize is 3 times larger than the one used in CSIDH and PEGASIS, and uses an order $\mathcal{O}$ of discriminant $\Delta$ with $|\Delta|^3 < 2^e < p$. This highly affects the performance and the key sizes in KLaPoTi. Even though PEGASIS performs the group action using 4-dimensional isogenies, the fact that $|\Delta| \approx 2^e \approx p$ in PEGASIS enables a group action computation which is more efficient. For example, to achieve similar security to CSIDH-512, KLaPoTi uses a prime *p* of 1536 bits and a single group action computation takes approximately 3 minutes [50, Section 7.1] in SageMath, while PEGASIS uses a prime of about 512 bits and a single group action computation takes approximately 1.5 seconds.

**Organisation of paper.** The rest of the paper is organized as follows. In Section 2 we recall the necessary background material. In Section 3 we describe our whole algorithm for general orientations. In Section 4, we detail an important subroutine, namely how to compute an isogeny of prescribed degree in dimension 2, before we in Section 5 specialise our algorithm for the CSIDH orientation. Finally, in Section 6 we detail some implementation choices before we present our timings.

**Acknowledgements.** This work was supported in part by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (under grant agreement ISOCRYPT - No. 101020788, and grant agreement AGATHA CRYPTOY - No. 101116169); by the Research Council KU Leuven grant C14/24/099 and by CyberSecurity Research Flanders with reference number VR20192203; by the Agence Nationale de la Recherche under grant ANR-22-PNCQ-0002 (HQI); by SNSF Consolidator Grant Crypto-nIs 213766. R. Invernizzi is funded by Research Foundation - Flanders (FWO) under a PhD Fellowship fundamental research (project number 1138925N).

## 2 Preliminaries

For a general introduction to isogenies, especially in the context of cryptography, we refer the reader to [31].
## 2.1 Isogeny class-group action on oriented curves

Classic results [65, Th. 4.5] tell us that the class-group $Cl(O)$ of an imaginary quadratic order $O$ acts freely and transitively on the set $Ord_q(O)$ of ordinary elliptic curves $E/\mathbb{F}_q$ with endomorphism ring $End(E) \cong O$. In [15], generalising the work of [13, 12], Colò and Kohel introduced the framework for a similar action, in which $Cl(O)$ acts on a set of supersingular elliptic curves $E$ which contain $O$ as a subring of $End(E)$. We briefly recount this general construction, together with results from [48].

Let $O = \mathbb{Z}[\alpha]$ be an imaginary quadratic order and $E/\mathbb{F}_{p^n}$ an elliptic curve. We call an injective morphism of rings $\iota: O \to End(E)$ an *O-orientation* and say it is *primitive* if $\iota(\alpha) = \omega$ is a primitive element of $End(E)$ viewed as a lattice, i.e. $\omega$ does not factor through a (non-trivial) scalar multiplication. We call the pair $(E, \iota)$ a (primitively) *O-oriented* elliptic curve.

Let $(E, \iota)$ be primitively *O-oriented* and $\mathfrak{a}$ an integral invertible ideal of $O$ with norm coprime to $p$. Let $\varphi_{\mathfrak{a}}: E \to E_{\mathfrak{a}}$ be the isogeny with kernel

$$E[\mathfrak{a}] = \bigcap_{\sigma \in \mathfrak{a}} \ker(\iota(\sigma))$$

and degree $N(\mathfrak{a})$. Now consider the map $O \to End(E_{\mathfrak{a}}) : \gamma \mapsto \varphi_{\mathfrak{a}}\iota(\gamma)\widehat{\varphi_{\mathfrak{a}}}$. It is an injective morphism of additive groups, but sends $1 \mapsto [\deg(\varphi_{\mathfrak{a}})]$ and so is not a morphism of rings. We fix this by first noting that $\varphi_{\mathfrak{a}}\iota(\gamma)\widehat{\varphi_{\mathfrak{a}}}$ always factors through $[\deg(\varphi_{\mathfrak{a}})]$ and then defining $(\varphi_{\mathfrak{a}})_*(\iota)$ as the map given by $(\varphi_{\mathfrak{a}})_*(\iota)(\gamma) = \varphi_{\mathfrak{a}}\iota(\gamma)\widehat{\varphi_{\mathfrak{a}}}/[\deg(\varphi_{\mathfrak{a}})]$ to obtain a morphism of rings and thus a well-defined *O-orientation*. Indeed, by diagram-chasing one verifies that $\ker(\varphi_{\mathfrak{a}})$ is invariant under $\omega = \iota(\alpha)$, demonstrating that $\varphi_{\mathfrak{a}}\iota(\gamma)\widehat{\varphi_{\mathfrak{a}}}(E[\deg(\varphi_{\mathfrak{a}})]) = \{0\}$ for all $\gamma$ in $\mathbb{Z}[\alpha]$, and that $\varphi_{\mathfrak{a}}\iota(\gamma)\widehat{\varphi_{\mathfrak{a}}}$ factors through $[\deg(\varphi_{\mathfrak{a}})]$. In fact, $(\varphi_{\mathfrak{a}})_*(\iota)$ is a *primitive* orientation [48, Prop. 3.5].

We say that the isogenies $\varphi_{\mathfrak{a}}: (E, \iota) \to (E_{\mathfrak{a}}, (\varphi_{\mathfrak{a}})_*(\iota))$ are *O-oriented*. Likewise, an isomorphism $\varphi: E \to E'$ is called an *O-isomorphism* $(E, \iota) \to (E', \iota')$ if $\iota' = \varphi_*(\iota) := (\gamma \mapsto \varphi\iota(\gamma)\hat{\varphi})$. Using this terminology, we define $SS_p^{\mathrm{pr}}(O)$ to be the *O-isomorphism* classes of primitively *O-oriented* supersingular elliptic curves defined over $\mathbb{F}_{p^n}$. This set is non-empty if and only if $p$ is not split in $K = O \otimes_{\mathbb{Z}} \mathbb{Q}$ and does not divide the conductor of $O$ in $O_K$ [48, Prop. 3.2].

We note that every class $[\mathfrak{a}]$ in the class-group $Cl(O)$ is represented by an invertible integral ideal $\mathfrak{b}$ with norm $N(\mathfrak{b})$ coprime to $p$ [17, Cor. 7.17]. Moreover, it is clear that principal ideals $\gamma O$ correspond to endomorphisms $\varphi_{\gamma O} = \iota(\gamma)$. Now the map

$$\mathrm{Cl}(O) \times \mathrm{SS}_p^{\mathrm{pr}}(O) \to \mathrm{SS}_p^{\mathrm{pr}}(O) \quad [\mathfrak{a}], (E, \iota) \mapsto \mathfrak{a} \cdot (E, \iota) := (E_{\mathfrak{a}}, (\varphi_{\mathfrak{a}})_*(\iota))$$

is well-defined; and when $SS_p^{\mathrm{pr}}(O)$ is non-empty, the action is free with at most two orbits [48, Prop. 3.3, Th. 3.4]. Restricting to one of these orbits, we obtain a free and transitive group action.

The vectorisation problem of this action - given $(E, \iota), (E', \iota')$ compute $[\mathfrak{a}]$ such that $\mathfrak{a} \cdot (E, \iota) = (E', \iota')$ - is thought to be quantum resistant. This property
led to several cryptographic constructions based on ideal class group actions [13, 15, 21, 14] which define (restricted) effective group actions, as mentioned in Section 1. The first one to be introduced was CSIDH (Commutative Supersingular Isogeny Diffie-Hellman)<fnref n="8" /> [13] where $\mathcal{O} = \mathbb{Z}[\sqrt{-p}]$ and the orientation induces an isomorphism $\mathcal{O} \rightsquigarrow \mathrm{End}_{\mathbb{F}_p}(\mathcal{O})$, mapping $\sqrt{-p}$ to the Frobenius endomorphism of $(x, y) \mapsto (x^p, y^p)$ for all $E \in \mathrm{SS}_p^{\mathrm{pr}}(\mathcal{O})$. Even though these schemes are thought to be quantum resistant, the best known quantum attack against the underlying vectorisation problem is a subexponential algorithm due to Kuperberg [39]. For that reason, there is an uncertainty on the parameters required to ensure a sufficient security level for CSIDH (and other class group action schemes). As mentioned in Section 1, conservative sources [6, 9, 51] estimate that a 2000 bit prime $p$ is required to ensure a NIST-1 security level for CSIDH.

## 2.2 Polarised isogenies between abelian varieties

Every abelian variety $A$ has a dual $\hat{A}$ of the same dimension and every isogeny $\varphi: A \to B$ has a dual $\hat{\varphi}: \hat{B} \to \hat{A}$ of the same degree. If there exists a special kind of isogeny $\lambda_A: A \to \hat{A}$ called a *polarisation*, we say $(A, \lambda_A)$ is a *polarised* abelian variety (PAV). When $\lambda_A$ is an isomorphism, we say it is a *principal* polarisation and call $A$ a *principally polarised* abelian variety (PPAV). Elliptic curves are exactly the 1-dimensional PPAVs; they are canonically principally polarised.

**Definition 2.1 (Polarised Isogenies).** Let $\varphi: (A, \lambda_A) \to (B, \lambda_B)$ be an isogeny between PAVs. If there exists an integer $d$ such that $\hat{\varphi}\lambda_B\varphi = [d]\lambda_A$, we say that $\varphi$ is $(\lambda_A, \lambda_B)$-polarised with polarised degree $\mathrm{deg}_p(\varphi) = d$. When $A, B$ are PPAVs we define the polarised dual of $\varphi$ to be $\tilde{\varphi} = \lambda_A^{-1}\hat{\varphi}\lambda_B: (B, \lambda_B) \to (A, \lambda_A)$. Then $\varphi$ is $(\lambda_A, \lambda_B)$-polarised with degree $d$ if $\tilde{\varphi}\varphi = [d]$ in $\mathrm{End}(A)$.

($\lambda_A, \lambda_B$)-polarised isogenies of polarised degree $d$ are characterised by an isotropy condition of their kernel.

**Lemma 2.1.** Let $\varphi: (A, \lambda_A) \to (B, \lambda_B)$ be a $(\lambda_A, \lambda_B)$-polarised isogeny with polarised degree $d$ such that $p \nmid d$ and let $g := \mathrm{dim}(A)$. Then $\mathrm{ker}(\varphi)$ is a subgroup of cardinality $d^g$ of $A[d]$, isotropic with respect to the Weil pairing $e^{[d]\lambda_A}$ associated to the polarisation $[d]\lambda_A: A \to \hat{A}$ (as defined in [41, p. 135]).

Conversely, if $A$ is a PPAV of dimension $g$ and $\varphi: A \to B$ has kernel $\mathrm{ker}(\varphi) \subseteq A[d]$ of cardinality $d^g$ and isotropic for $e^{[d]\lambda_A}$, then there exists a unique principal polarisation $\lambda_B$ on $B$ such that $\varphi: A \to B$ is a $(\lambda_A, \lambda_B)$-polarised isogeny with polarised degree $d$.

*Proof.* Suppose that $\varphi$ is a $(\lambda_A, \lambda_B)$-polarised isogeny with polarised degree $d$. Let $x, y \in \mathrm{ker}(\varphi)$ and $y' \in B$ such that $y = \tilde{\varphi}(y')$ (which exists by surjectivity

<footnote n="8">⁸ Note that CSIDH is based on the CRS protocol [16, 60], which uses the full endomorphism ring of an ordinary curve instead of an orientation of a supersingular curve</footnote>
of $\tilde{\varphi}$ and belongs to $B[d]$ since $\varphi \circ \tilde{\varphi} = [d])$. Then, by the definition of $e^{[d]\lambda_A}$ and the usual properties of the Weil pairing [41, Lemma 16.2],

$$ 
\begin{aligned}
e^{[d]\lambda_A}(x, y) &= e_d(x, \lambda_A(y)) = e_d(x, \lambda_A \circ \tilde{\varphi}(y')) = e_d(x, \tilde{\varphi} \circ \lambda_B(y')) \\
&= e_d(\varphi(x), \lambda_B(y')) = 1,
\end{aligned}
$$

since $\varphi(x) = 0$. This proves that $\ker(\varphi)$ is isotropic.

Conversely, suppose that $\ker(\varphi) \subseteq A[d]$ of cardinality $d^g$ and isotropic for $e^{[d]\lambda_A}$. Then [41, Proposition 16.8] ensures the existence of a polarisation $\lambda_B$ on $B$ such that $\tilde{\varphi} \circ \lambda_B \circ \varphi = [d]\lambda_A$. Since $\deg(\tilde{\varphi}) = \deg(\varphi) = |\ker(\varphi)| = d^g$ and $\deg([d]) = d^{2g}$, it follows that $\deg(\lambda_B) = \deg(\lambda_A) = 1$ by multiplicativity of the degrees so $\lambda_B$ is a principal polarisation. This completes the proof. $\square$

For brevity, we will drop mention of the polarisations when they are clear from context, and say *d-isogeny* to mean an isogeny with polarised degree *d*. In particular, it is implicit that the corresponding polarisation to a PAV denoted $A$ is $\lambda_A$. Between elliptic curves, the principal polarisation being canonical, the dual $\tilde{\varphi}$ and polarised dual $\tilde{\varphi}$ of an isogeny can be identified and will both be denoted by $\tilde{\varphi}$; and the notions of degree and polarised degree are also the same.

We recall a standard result about factoring isogenies between PPAVs, which we will need to recall the Clapoti construction in the more general context of [49, Remark 2.10].

**Lemma 2.2 (Factoring Isogenies between PPAVs).** Let $A, B$ be PPAVs defined over a field $k$ and $\varphi: A \to B$ be a $d$-isogeny between them. For every pair of positive coprime integers $d_1, d_2$ each being coprime to $\text{char}(k)$ satisfying $d_1d_2 = d$, there exists a PPAV $C$, a $d_1$-isogeny $\varphi_1: A \to C$ and a $d_2$-isogeny $\varphi_2: C \to B$, each uniquely defined up to isomorphism, such that $\varphi = \varphi_2 \circ \varphi_1$.

*Proof.* Since $d_1$ and $d_2$ are coprime, we have $\ker(\varphi) = \ker(\varphi)[d_1] \oplus \ker(\varphi)[d_2]$ with $|\ker(\varphi)[d_i]| = d_i^g$ for $i \in \{1, 2\}$. Consider the isogeny $\varphi_1: A \to C$ of kernel $\ker(\varphi)[d_1]$ and $\varphi_2: C \to B'$ the isogeny of kernel $\varphi_1(\ker(\varphi))$. Then by construction $\ker(\varphi_2 \circ \varphi_1) = \ker(\varphi)$ so $\varphi_2 \circ \varphi_1 = \varphi$ up to postcomposition with an isomorphism $B \simeq B'$ by [45, Theorem 4, p. 74].

Since $\varphi$ is a $d$-isogeny, Lemma 2.1 ensures that $\ker(\varphi) \subseteq A[d]$ is isotropic for $e^{[d]\lambda_A}$ so $\ker(\varphi_1) \subseteq A[d_1]$ is isotropic for $e^{[d_1]\lambda_A}$ by the usual properties of the Weil pairing [41, Lemma 16.1] and of cardinality $d_1^g$, as we saw. Using Lemma 2.1 again, we obtain a principal polarisation $\lambda_C$ on $C$ such that $\varphi_1$ is a $(\lambda_A, \lambda_C)$-polarised isogeny of polarised degree $d_1$.

To conclude, we verify that $\varphi_2$ is a $(\lambda_C, \lambda_B)$-polarised isogeny of polarised degree $d_2$. Substituting into $\tilde{\varphi}\lambda_B\varphi = [d_1d_2]\lambda_A$, we obtain $\tilde{\varphi}_1\tilde{\varphi}_2\lambda_B\varphi_2\varphi_1 = [d_2d_1]\lambda_A = [d_2]\tilde{\varphi}_1\lambda_C\varphi_1$ and conclude that $\tilde{\varphi}_2\lambda_B\varphi_2 = [d_2]\lambda_C$. Indeed, isogenies are surjective so $fh = gh$ implies $f = g$ and by forming the dual (twice) we see that $hf = hg$ also implies $f = g$. Note that this last argument also proves the uniqueness of $\varphi_1$ and $\varphi_2$. This completes the proof. $\square$
## 2.3 Kani's lemma and embedding isogenies into higher degree

If $(\varphi_{ji})_{ij}$ is a matrix of polarised isogenies $\varphi_{ij}: A_i \to B_j$ between PPAVs, then its polarised dual is given by $(\tilde{\varphi}_{ij})_{ij}$. Whether or not a $2 \times 2$ matrix of polarised isogenies is a polarised isogeny is the content of Kani's lemma.

**Lemma 2.3 (Kani, [54, Lemma 2.1]).** Let $A_1, A_2, B_1, B_2$ be PPAVs defined over $k$ and let $\varphi_{ij}: A_i \to B_j$ be polarised $d_{ij}$-isogenies. The map

$$ \Phi = \begin{pmatrix} \varphi_{11} & \varphi_{21} \\ \varphi_{12} & \varphi_{22} \end{pmatrix} : A_1 \times A_2 \to B_1 \times B_2 $$

is a polarised isogeny if and only if $d_{11} = d_{22}, d_{12} = d_{21}$ and the corresponding square $S_{\Phi}$

$$ \begin{array}{ccc} A_1 & \xrightarrow{\varphi_{11}} & B_1 \\ -\varphi_{12} & \downarrow & \downarrow \widetilde{\varphi_{21}} \\ B_2 & \xrightarrow{\widetilde{\varphi_{22}}} & A_2 \end{array} $$

commutes. In such cases, $\deg_p(\Phi) = d_{11} + d_{12}$ and we call $\Phi$ the Kani-map corresponding to the Kani-square $S_{\Phi}$. Further, if $d_{11} = d_{22}$ is coprime to $d_{12} = d_{21}$, we call the Kani-square $a$ $(d_{11}, d_{12})$-isogeny diamond. Finally, when $d_{11}$ is coprime to $d_{12}$ and $d = d_{11} + d_{12}$ is coprime to the characteristic of $k$, then

$$ \ker(\Phi) = \{([d_{11}]x, \widetilde{\varphi_{21}}\varphi_{11}(x)) \mid x \in A_1[d]\}. $$

If $\Phi$ can be efficiently evaluated, then so can $\varphi_{ij}$, and $\Phi$ is said to embed $\varphi_{ij}$. In other words, $\Phi$ is an efficient representation of the $\varphi_{ij}$ in the sense of the following definition.

**Definition 2.2.** Let $\mathcal{A}$ be an algorithm. Given a $d$-isogeny $\varphi: A \to B$ between PPAVs of dimension $g$ defined over $\mathbb{F}_q$, an efficient representation of $\varphi$ with respect to $\mathcal{A}$ is some data $D_{\varphi} \in \{0, 1\}^*$ of polynomial size in $g$, $\log(d)$ and $\log(q)$ such that for all $P \in A(\mathbb{F}_{q^k})$, $\mathcal{A}$ with input $(D_{\varphi}, P)$ returns $\varphi(P)$ in polynomial time in $g$, $\log(d)$, $k$ and $\log(q)$.

# 3 Our algorithmic approach

We now describe our algorithm for computing the class group action by a general orientation.

## 3.1 Applying Kani's lemma in dimension 4

The factorisation lemma (Lemma 2.2) gives us a strategy to embed arbitrary isogenies into higher dimension [49, Remark 2.10]. If $\varphi: A \to B$, $\psi: A \to C$
are isogenies between PPAVs with coprime polarised degrees, we can form a Kani-diamond by factoring the composition $\psi\tilde{\phi}$ as $\tilde{\phi}'\psi'$ whereby $\psi': B \to D$ and $\phi': C \to D$ have polarised degrees $\deg_p(\phi') = \deg_p(\phi)$, $\deg_p(\psi) = \deg_p(\psi')$. The kernel of the resulting Kani-map can be computed so long as $\psi\tilde{\phi}$ is known on $B[\deg_p(\phi) + \deg_p(\psi)]$.

When computing the class-group action of $[\mathfrak{a}]$ on a primitively oriented elliptic curve $(E, \iota)$, the isogenies $\phi, \psi$ are taken to be $\hat{\phi}_{\mathfrak{b}}, \hat{\phi}_{\mathfrak{c}}: E_{\mathfrak{a}} \to E$ where $\mathfrak{b}, \mathfrak{c}$ are representatives of $[\mathfrak{a}]$. The composition $\hat{\phi}_{\mathfrak{c}}\phi_{\mathfrak{b}} = \phi_{\tau_{\mathfrak{b}}}$ is an endomorphism on $E$ dictated by $\iota$ and so the kernel of the resulting 2-dimensional Kani-map can be computed so long as the $N(\mathfrak{b})+N(\mathfrak{c})$ torsion is accessible. For practically efficient computation (using [20] for instance), we must find ideals $\mathfrak{b}, \mathfrak{c}$ representing $[\mathfrak{a}]$ such that $N(\mathfrak{b}) + N(\mathfrak{c}) = 2^e$, where $2^{e+2} \mid p+1$ in order to work over $\mathbb{F}_{p^2}$ (the reason we need $2^{e+2}$ to divide $p+1$ is due to technicalities when using level-2 theta coordinates to compute higher-dimensional isogenies). However, finding ideals $\mathfrak{b}, \mathfrak{c}$ satisfying this norm equation is not easy in practice and can only be done with standard techniques when $p$ is much bigger than $|\mathrm{disc}(\mathcal{O})|$ (which is not suitable for CSIDH/CSURF, where $|\mathrm{disc}(\mathcal{O})| = p$ or $4p$).

**Goal:** solve the norm equation:

$$uN(\mathfrak{b}) + vN(\mathfrak{c}) = 2^e \quad (3)$$

with $2^{e+2} \mid p+1$, $[\mathfrak{a}] = [\mathfrak{b}] = [\mathfrak{c}]$, $u, v > 0$, $\gcd(uN(\mathfrak{b}), vN(\mathfrak{c})) = 1$, and such that we can build efficiently computable dimension-d isogenies $\Phi_u: E^d \to A_u, \Phi_v: E^d \to A_v$ of polarised degrees $u, v$.

Indeed, the factoring lemma (Lemma 2.2) ensures that we have a $(uN(\mathfrak{b}), vN(\mathfrak{c}))$-isogeny diamond

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

where $\Phi_{\mathfrak{b}} := \mathrm{diag}(\phi_{\mathfrak{b}}), \Phi_{\mathfrak{c}} := \mathrm{diag}(\phi_{\mathfrak{c}}): E^d \to E_{\mathfrak{a}}^d, \Phi$ is an $uN(\mathfrak{b})$-isogeny and $\Psi$ is a $vN(\mathfrak{c})$-isogeny. Hence, by Kani's Lemma 2.3, we can embed the compositions $\Phi_{\mathfrak{b}} \circ \tilde{\Phi}_{u}$ and $\Phi_{\mathfrak{c}} \circ \tilde{\Phi}_{v}$ into a 2d-dimensional $2^e$-isogeny $F: A_u \times A_v \to E_{\mathfrak{a}}^d \times A$.

Except when $u$ and $v$ are perfect squares as in [50] (which does not make the norm equation much easier to solve), we cannot expect $\Phi_u$ and $\Phi_v$ to be one-dimensional. Using higher dimensions $d \ge 2$ appears necessary because it is difficult to construct (efficient representations of) outgoing 1-dimensional isogenies $\Phi_u, \Phi_v$ of large given degree $u, v$ on $E$ without knowing the endomorphism ring $\mathrm{End}(E)$: the only tools we have at our disposal to produce arbitrary-degree
isogenies all require higher dimensions. In fact, a variant of this problem (construct a large prime degree isogeny on $E$) is conjectured to be cryptographically hard [4, Problem 2].

In any case, Zarhin's trick [28, Th. 11.29, Eq. (4)] would allow us to efficiently construct endomorphisms $\Phi_u, \Phi_v$ of $E^4$ with the prescribed polarised degrees $u, v$. However, we consider computing the corresponding 8-dimensional isogeny $F$ to be too expensive. Instead, we will see in Section 4 how to find solutions when $d = 2$. Consequently, we have to compute a 4-dimensional isogeny $F$, which can be done efficiently with level-2 theta coordinates [18]. As mentioned earlier, the algorithms to compute the $2^e$-isogeny $F$ require $2^{e+2}$-torsion points, so we need $2^{e+2}|p+1$ instead of $2^e|p+1$ to be able to work over $\mathbb{F}_{p^2}$.

Unfortunately, there do not always exist positive integer solutions $u, v$ to $uN(\mathfrak{b}) + vN(\mathfrak{c}) = 2^e$. Indeed, the Frobenius coin problem tells us that we are only guaranteed solutions when $2^e \ge N(\mathfrak{b})N(\mathfrak{c}) - N(\mathfrak{b}) - N(\mathfrak{c})$; on the other hand, when $\mathcal{O} = \mathbb{Z}[\sqrt{-p}]$ (as in CSIDH), the smallest representative of any class $[\mathfrak{a}]$ in $Cl(\mathcal{O})$ is expected to be of size roughly $\sqrt{p}$. So, recalling that $2^{e+2} | p+1$, we see that we rarely have guaranteed solutions $u, v$. The solution we propose to overcome this difficulty is to give some additional freedom on the choice of ideals $\mathfrak{b}, \mathfrak{c}$ by factoring out their smooth part.

## 3.2 Factoring out easy steps first

We will now describe how to compute the action of $[\mathfrak{a}]$ on $(E, \iota)$ in four dimensions. Throughout this section, $\mathfrak{b}, \mathfrak{c}$ are primitive$^9$ representatives of $[\mathfrak{a}]$.

**The idea.** Suppose we can factor $\mathfrak{b} = \mathfrak{b}_e\mathfrak{b}_k, \mathfrak{c} = \mathfrak{c}_e\mathfrak{c}_k$, so that the isogenies $\varphi_{\mathfrak{b}_e}: E \to E_1, \varphi_{\mathfrak{c}_e}: E \to E_2$ are efficiently computable (e.g. smooth degree) and so that we can find efficiently computable isogenies $\Phi_u: E_1^2 \to A_u, \Phi_v: E_2^2 \to A_v$ (e.g. integer-matrix endomorphisms) with polarised degrees $u, v$ satisfying

$$uN(\mathfrak{b}_k) + vN(\mathfrak{c}_k) = 2^e. \qquad (4)$$

Then we can consider the following $(uN(\mathfrak{b}), vN(\mathfrak{c}))$-isogeny diamond

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

<footnote>$^9$ i.e. neither $\mathfrak{b}$ nor $\mathfrak{c}$ is contained in $NO$ for any $N > 1$.</footnote>
where $\Phi_I$ is the 2-dimensional isogeny $\mathrm{diag}(\varphi_I, \varphi_I)$ given by the action of $I \subseteq O$.
Then, Kani's lemma yields a 4-dimensional isogeny

$$ F := \begin{pmatrix} \Phi_{b_k} \circ \tilde{\Phi}_u \Phi_{c_k} \circ \tilde{\Phi}_v \\ -\tilde{\Psi} & \tilde{\Phi} \end{pmatrix} : A_u \times A_v \to E_a^2 \times A, \qquad (5) $$

with kernel

$$ \begin{aligned} \ker(F) &= \{([uN(\mathfrak{b}_k)]x, \Phi_v \tilde{\Phi}_{c_k} \tilde{\Phi}_{b_k} \tilde{\Phi}_u(x)) \mid x \in A_u[2^e]\} \\ &= \{([uN(\mathfrak{b}_k)]x, \Phi_v \tilde{\Phi}_{\bar{c}_k b_k} \tilde{\Phi}_u(x)) \mid x \in A_u[2^e]\} \\ &= \{([N(\mathfrak{b}_k)]\tilde{\Phi}_u(y), \Phi_v \tilde{\Phi}_{\bar{c}_k b_k}(y)) \mid y \in E_1^2[2^e]\}. \end{aligned} \qquad (6) $$

Since $\mathfrak{b}, \mathfrak{c}$ both represent $[\mathfrak{a}]$, $\bar{\mathfrak{c}}\mathfrak{b} = \sigma O$ is principal and we can evaluate $\varphi_{\bar{\mathfrak{c}}\mathfrak{b}} = \iota(\sigma)$.
Having computed $\varphi_{c_e}$ and $\varphi_{b_e}$, we can evaluate

$$ N(\mathfrak{b}_e)N(\mathfrak{c}_e)\varphi_{\bar{\mathfrak{c}}\mathfrak{k}\mathfrak{b}_k} = \varphi_{c_e}\iota(\sigma)\tilde{\varphi}_{\mathfrak{b}_e}. \qquad (7) $$

and compute $\ker(F)$ (more exactly points of $2^{e+2}$-torsion generating $\ker(F)$ as required in [18]). We refer to Section 5.2 for more details on the kernel computation in the case of CSIDH.

**The algorithm.** By construction, $uN(\mathfrak{b}_k) + vN(\mathfrak{c}_k) = 2^e$ has at least as many positive integer solutions $u, v$ as $uN(\mathfrak{b}) + vN(\mathfrak{c}) = 2^e$. Indeed, every ideal appearing as $\mathfrak{b}_k$ can also be chosen as $\mathfrak{b}$ (with $\mathfrak{b}_e = 1$). More relevant, though, is that the norms of the ideals $\mathfrak{b}_k$ are smaller, and so we hope for more $u, v$ solutions.

The ideals $\mathfrak{b}_e, \mathfrak{c}_e$ should be chosen so that the corresponding isogenies $\varphi_{\mathfrak{b}_e}, \varphi_{\mathfrak{c}_e}$ are easy to compute, e.g. of smooth degrees. More precisely, letting $\mathfrak{B}$ be a set of small primes, we permit $\mathfrak{b}_e, \mathfrak{c}_e$ to be a product of ideals lying above the primes in $\mathfrak{B}$. Note that since $\mathfrak{b}_e, \mathfrak{c}_e$ are primitive, these primes will all be split in $O$.

With this requirement we now proceed as follows (we refer to Section 5.1 for details).

1. Obtain a set $S$ of short representatives of $[\mathfrak{a}]$ (e.g. by Lagrange-reduction).

2. Factor each $\mathfrak{b}$ in $S$ as $\mathfrak{b} = \mathfrak{b}_e\mathfrak{b}_k$ maximising the norm of $\mathfrak{b}_e$.

3. Select a pair of ideals $\mathfrak{b} = \mathfrak{b}_e\mathfrak{b}_k$ and $\mathfrak{c} = \mathfrak{c}_e\mathfrak{c}_k$ from $S$ with $N(\mathfrak{b}_k), N(\mathfrak{c}_k)$ coprime, prioritizing those with the smallest $N(\mathfrak{b}_k), N(\mathfrak{c}_k)$.

4. Enumerate the positive integer solutions $u, v$ of $uN(\mathfrak{b}_k) + vN(\mathfrak{c}_k) = 2^e$.

5. Return a solution $u, v$ if it allows for the construction of $\Phi_u, \Phi_v$; else choose a different $\mathfrak{b}, \mathfrak{c}$ pair.

We briefly remark that this technique of splitting the isogeny computation into an easy 1-dimensional, followed by potentially expensive ($d \ge 2$)-dimensional computation also allows us to compute more isogenies $\Phi_u$ of prescribed polarised degree $u$. Indeed, the naïve approach of constructing $\Phi_u$ as a $2 \times 2$ integer-matrix

$$ \begin{pmatrix} x_u - y_u \\ y_u x_u \end{pmatrix} \qquad (8) $$
only permits $u = x_u^2 + y_u^2$ to be a sum of two squares. Composing $\Phi_u$ as the
diagonal of some one-dimensional isogeny of sufficiently smooth degree $g_u$ (e.g.
$\mathfrak{B}$-smooth), and then a matrix as above allows for polarised degrees $u = g_u(x_u^2 +
y_u^2)$ (see Section 4.1).

3.3 Finding the codomain orientation

Once we have computed the 4-dimensional isogeny $F : A_u \times A_v \to E_a^2 \times A$ and extracted $E_a$ from its codomain (following the approach presented in Appendix B in the case of CSIDH), we still have to compute the orientation $\iota_a := (\varphi_a)_*(\iota)$ on $E_a$. In the case of CSIDH, this orientation is naturally given by the Frobenius endomorphism and does not need to be computed. However, this orientation is not that easy to determine in general.

Given a generator $\alpha \in \mathcal{O}$ of $\mathcal{O}$ (i.e. $\mathcal{O} = \mathbb{Z}[\alpha]$), we only have to compute an efficient representation of $\iota_a(\alpha)$, as defined in Definition 2.2. Assuming $|\mathrm{disc}(\mathcal{O})| \le 2^{2f}$ with $2^f | p+1$, we can always find $\alpha$ of norm bounded by $2^{2f}$ and it suffices to evaluate $\iota_a(\alpha)$ on a basis of $E_a[2^f]$ to efficiently represent this endomorphism, e.g. using the higher dimensional interpolation algorithms from [56, Theorem 3.7 & Section 6.4]. Furthermore, if $e+2 \le f$, then knowing $\iota(\alpha)$ on the $2^f$-torsion is sufficient to compute the $2^{e+2}$-torsion lying above $\ker(F)$ given by Equation (6) so to compute the action of any ideal on $E$. The same holds for $E_a$: the knowledge of $\iota_a(\alpha)$ on $E_a[2^f]$ is sufficient to propagate the ideal class group action further.

Now we explain how we compute $\iota_a(\alpha)$ on a basis of $E_a[2^f]$. Using the fact
that $\overline{cb} = \sigma\mathcal{O}$, we obtain that

$$
\iota_a(\alpha) = (\varphi_a)_*(\iota) = \frac{1}{N(\mathfrak{a})} \varphi_a \circ \iota(\alpha) \circ \hat{\varphi}_a = \frac{1}{N(\mathfrak{b})} \varphi_b \circ \iota(\alpha) \circ \hat{\varphi}_b.
$$

Let $(P_1, Q_1)$ and $(P_a, Q_a)$ be two bases of $E_1[2^f]$ and $E_a[2^f]$ respectively. Then
by Eq. (5), we have $F(\Phi_u(P_1, 0), 0) = ([u]\varphi_{b_k}(P_1), 0, *)$ so we can evaluate
$[u]\varphi_{b_k}(P_1)$, hence $\varphi_{b_k}(P_1)$ by multiplying by an inverse of $u$ mod $2^f$. We com-
pute $\varphi_{b_k}(Q_1)$ similarly. Now, with (easy) discrete logarithm computations in
$E_a[2^f]$, we may find $a, b, c, d \in \mathbb{Z}$ such that $\varphi_{b_k}(P_1) = [a]P_a + [b]Q_a$ and $\varphi_{b_k}(Q_1) =$
$[c]P_a + [d]Q_a$. Then, a simple matrix inversion yields

$$
\hat{\varphi}_{b_k}(P_a) = [\delta \mathrm{N}(\mathfrak{b}_k)]([d]P_1 - [b]Q_1) \text{ and } \hat{\varphi}_{b_k}(Q_a) = [\delta \mathrm{N}(\mathfrak{b}_k)](-[c]P_1 + [a]Q_1),
$$

with $\delta$ an inverse of $ad-bc$ mod $2^f$ (which does exist because $N(\mathfrak{b}_k) = \deg(\varphi_{\mathfrak{b}_k})$ is odd, so $\varphi_{\mathfrak{b}_k}$ maps a $2^f$-torsion basis to a $2^f$-torsion basis). So we have computed $(\hat{\varphi}_{\mathfrak{b}_k}(P_a), \hat{\varphi}_{\mathfrak{b}_k}(Q_a))$. Similarly, knowing $\varphi_{\mathfrak{b}_e}$, we can compute the action of its dual $\hat{\varphi}_{\mathfrak{b}_e}$ on the $2^f$-torsion provided that $N(\mathfrak{b}_e)$ is odd. We can impose this condition by removing 2 from the set of allowed primes $\mathfrak{B}$. Hence, we can evaluate $\hat{\varphi}_{\mathfrak{b}_e} \circ \hat{\varphi}_{\mathfrak{b}_k} = \hat{\varphi}_b$ on the basis $(P_a, Q_a)$. We can then evaluate $[N(\mathfrak{b})]_{\iota_a}(\alpha) = \varphi_b \circ \iota(\alpha) \circ \hat{\varphi}_b = \varphi_{b_k} \circ \varphi_{b_e} \circ \iota(\alpha) \circ \hat{\varphi}_{b_e} \circ \hat{\varphi}_{b_k}$ on $(P_a, Q_a)$, using our knowledge of the action of $\iota(\alpha)$ on the $2^f$-torsion, of $\varphi_{b_e}$ and of $\varphi_{b_k}$ via $F$. Since $N(\mathfrak{b}) = N(\mathfrak{b}_e)N(\mathfrak{b}_k)$ is odd, we can invert it modulo $2^f$ and finally obtain $(\iota_a(\alpha)(P_a), \iota_a(\alpha)(Q_a))$, as desired.
# 4 Computing 2-dimensional isogenies of prescribed polarised degree

In this section, we explain how to construct dimension 2 isogenies $\Phi_u, \Phi_v$ of prescribed polarised degree $u, v$, in order to solve Eq. (3) (or more precisely the version in Eq. (4)).

We first treat an easy special case in Section 4.1. Then we explain in Section 4.2 how to treat the case of $u$ (or $v$) an arbitrary positive integer. The drawback of the latter approach is that it involves computing a 4-dimensional $2^e$-isogeny to build $\Phi_u$ in dimension 2.

## 4.1 The simplest case: when the polarised degree is $\mathfrak{B}$-good

Let us first assume that $u$ and $v$ are a sum of two squares. Then, given a curve $E$ and an integer $u = x_u^2 + y_u^2$, the endomorphism

$$M_u = \begin{pmatrix} x_u - y_u \\ y_u - x_u \end{pmatrix} \qquad (9)$$

on $E^2$ is such that $M_u \tilde{M}_u = [u]$ is the multiplication by $u$ on each component. Hence $M_u$ induces an isogeny $\Phi_u$ of polarised degree $u$ as desired.

The probability that a random number $u$ can be written as sum of two squares is roughly $1/\sqrt{\log(u)} \approx 1/\sqrt{\log(p)/2}$. Unfortunately, writing $u$ as above requires the factorisation of $u$. If we instead restrict $u$ to be a prime, the probability goes down to the order of $1/\log(u)$. Since we need both $u$ and $v$ to be of this form, we expect to find a valid pair after $\approx \log(p)^2$ attempts. Looking at the size of the primes we are working with, this approach becomes soon too expensive.

Instead, we can adopt the same strategy as in Section 3.2 to relax our condition and aim for $\mathfrak{B}$-good values $u, v$, for a given set of small primes.

**Definition 4.1.** Let $\mathfrak{B}$ be a set of small primes. We say that $u \in \mathbb{N}^*$ is $\mathfrak{B}$-good when $u$ is of the form $u = g_u(x_u^2 + y_u^2)$ where $x_u, y_u \in \mathbb{N}$ and $g_u \in \mathbb{N}^*$ is a product of primes in $\mathfrak{B}$.

We may then write $u = g_u(x_u^2 + y_u^2)$ and $v = g_v(x_v^2 + y_v^2)$, as above, where $g_u$ and $g_v$ contain all the prime factors of $u, v$ that are in $\mathfrak{B}$. Recall that an integer can be written as a sum of two squares if and only if, in its prime power decomposition, there is no prime power $\ell^k$ such that $k$ is odd and $\ell \equiv 3 \mod 4$ [37, Ch. 2.15, Ex. 10]. We can then restrict $g_u$ and $g_v$ to primes that are 3 mod 4. Notice that small prime numbers are more likely to appear in the factorization of a number, and hence to cause failure of the sum of squares condition. Heuristically, by taking out factors 3, 7 and 11 we already increase the probability of finding a sum of squares by a factor 3. If we choose $\mathcal{O}$ such that these small primes are split in $\mathcal{O}$, we can easily compute isogenies of degrees $g_u$ and $g_v$ via the action of $\mathrm{Cl}(\mathcal{O})$. In particular, in the case of CSIDH, these isogenies are defined over $\mathbb{F}_p$ so they are even more efficient to compute (see Section 5.2).
Given a pair $u, v$ with no factors in $\mathfrak{B}$, we are now left with the task of determining if both numbers can be expressed as sums of two squares. As mentioned above, checking this condition requires factoring, which we cannot afford. We instead perform trial division up to a fixed bound ($10^4$ for the 500 and 1000 parameters, and $10^5$ for all the others), and then test what is left for primality. We observe that this procedure reduces the chances of finding a good pair by 2 to 3 orders of magnitude, as expected.

If the primality test succeeds, and all nonsquare factors are 1 mod 4, we can apply Cornacchia's algorithm to obtain $x_u$ and $y_u$ (respectively $x_v, y_v$). This strategy is further discussed in Section 5.2 in the case of CSIDH.

**Remark 4.1.** It is possible to relax the $\mathfrak{B}$-good condition on $(u, v)$ by allowing $u$ (resp. $v$) to take the form $g_u(ax_u^2 + by_u^2)$ (resp. $g_v(ax_v^2 + by_v^2)$), where $g_u, g_v, a, b \in \mathfrak{B}$, while ensuring that Kani's lemma remains applicable as in Section 3.1. Then, at the cost of computing 4 additional isogenies of degree $a$ and $b$ with Elkies' algorithm, we improve the likelihood that $u$ and $v$ are suitable for computing isogenies of degree $u$ and $v$ in dimension 2. Experiments have shown that allowing $a$ and $b$ to take values in $\{2, 3, 5, 7\}$ increases the probability of success by approximately a factor of 3. Further details on this construction can be found in Appendix C.

## 4.2 The general case

We now treat the case of a general $u, v$. Although it is easy to construct $\Phi_u, \Phi_v$ in dimension 4, we saw in Section 3.1 that this would involve a dimension 8 isogeny for $F$. In this section we explain how to construct them in dimension 2, on $E^2$. For that, we extend the approach of QFESTA [46] from dimension 1 to dimension 2. A drawback is that our approach will be heuristic, although our implementation shows that the heuristic works well in practice.

Let $(E, \iota)$ be an $\mathcal{O}$-oriented supersingular elliptic curve. Let $\Delta = \text{disc}(\mathcal{O})$ be the discriminant of $\mathcal{O}$. First, we show that we can heuristically build (efficient) endomorphisms on $E^2$ of polarised degree $N$, as long as $N \gg |\Delta|$. By construction, we can evaluate any endomorphism $\gamma \in \mathcal{O}$, which gives us a 1-dimensional isogeny $\iota(\gamma)$ of degree $N(\gamma)$. Since $\sqrt{\Delta} \in \mathcal{O}$, the quadratic form $\gamma \mapsto N(\gamma)$ restricts on the suborder $\mathbb{Z}+\sqrt{\Delta}\mathbb{Z}$ of $\mathcal{O}$ to the quadratic form $q(x, y) = x^2+|\Delta|y^2$.

Now let $\gamma_1, \gamma_2 \in \mathcal{O}$ of norms $n_1, n_2$ respectively, and consider the endomorphism

$$ \gamma = \begin{pmatrix} \iota(\gamma_1) & \iota(\overline{\gamma_2}) \\ -\iota(\gamma_2) & \iota(\overline{\gamma_1}) \end{pmatrix} \in \text{End}(E^2). $$

Then, by [56, Lemma 3.2], its polarised dual is

$$ \tilde{\gamma} = \begin{pmatrix} \iota(\overline{\gamma_1}) & -\iota(\overline{\gamma_2}) \\ \iota(\gamma_2) & \iota(\gamma_1) \end{pmatrix}, $$

so that $\tilde{\gamma} \circ \gamma = [n_1+n_2]$ in $\text{End}(E^2)$ and $\gamma$ is an $(n_1+n_2)$-isogeny. In particular, we can build endomorphisms on $E^2$ of polarised degree $N$ as long as $N$ is represented by the quadratic form $q(x, y, z, t) = (x^2 + y^2) + |\Delta|(z^2 + t^2)$.
**Heuristic 1.** Let $N \gg |\Delta|$. Then we can heuristically build an efficient endo-morphism $\gamma \in \text{End}(E^2)$ of polarised degree $N$.

*Justification.* Since $N \gg |\Delta|$, we can sample many small $z, t$ such that $(z^2 + t^2)|\Delta| < N$. Then if $m = N - |\Delta|(z^2 + t^2)$ is a sum of two squares, we can find $\gamma$ by the discussion above. This situation can be efficiently detected whenever $m$ is a prime congruent to 1 mod 4, which heuristically happens with probability $\approx 1/\log(N)$. Since $N \gg |\Delta|$, we can try with many different couples $(z, t)$ until one succeeds. $\square$

In practice, we can use a direct adaptation of the algorithm **RepresentInteger** from [23] to find a solution of $(x^2 + y^2) + |\Delta|(z^2 + t^2) = N$.

Now by Lemma 2.2, if $\gamma$ is an efficiently represented endomorphism of $E^2$ of polarised degree $N = N_1N_2$, with $N_1, N_2$ coprime, then it decomposes uniquely as $\gamma = \mu_2 \circ \gamma_1 = \mu_1 \circ \gamma_2$ where the $\gamma_i, \mu_i$ are isogenies of polarised degree $N_i$. Furthermore, by Kani's lemma (Lemma 2.3) the 4-dimensional isogeny

$$ \Gamma = \left( \begin{array}{cc} \gamma_1 & \widetilde{\mu_2} \\ -\gamma_2 & \widetilde{\mu_1} \end{array} \right) \qquad (10) $$

is of polarised degree $N_1 + N_2$. The kernel of $\Gamma$ is given by $\text{Ker } \Gamma = \{(N_1P, \gamma P) \mid P \in E^g[N_1 + N_2]\}$, so we can efficiently evaluate $\Gamma$ to obtain the 2-dimensional isogeny $\gamma_1$ of polarised degree $N_1$ as long as the $(N_1 + N_2)$-torsion on $E$ is smooth and accessible (i.e. defined over a small extension of $\mathbb{F}_p$). Hence, an efficient representation of $\gamma$ yields an efficient representation of $\gamma_1$ in the sense of Definition 2.2.

For simplicity, and because this matches our implementation choice, we will now restrict to the case where we look for $N_1 + N_2 = 2^e$ a power of two.

**Theorem 4.1.** Assume that $E[2^e] \subseteq E(\mathbb{F}_{p^2})$ (with $e \ge 3$), that $2^{3e/2-1} \gg |\Delta|$, and that $u < 2^{e-1}$ is odd. Then, under Heuristic 1, we can efficiently represent an isogeny $\Phi_u$ of polarised degree $u$ on $E^2$.

*Proof.* First consider $N = N_1N_2$ with $N_1 = u, N_2 = 2^e - u$. Since $u$ is assumed to be odd, $N_1$ is coprime to $N_2$. If $N \gg |\Delta|$, we can use Heuristic 1 to build an endomorphism in dimension 2 of polarised degree $N$, which we split via a 4-dimensional isogeny of polarised degree $2^e$ to get $\Phi_u$.

If this is not the case, then since $2^{3e/2} \gg |\Delta|$ by assumption, this means that $u$ is smaller than $2^f$, with $f = \lceil e/2 \rceil$. We use a two step solution like in [47, 27]. First we let $u_2 = u(2^f - u)$, then $2^f - 1 \le u_2 < 2^{2f-2} \le 2^{e-1}$. Then we let $N = u_2(2^e - u_2)$. By our choice of $f$, $N \ge 2^{3e/2-1}$, so $N \gg |\Delta|$ and we can use Heuristic 1 to build an endomorphism in dimension 2 of polarised degree $N$, which we split a first time (via a 4-dimensional isogeny of polarised degree $2^e$) to get a 2-dimensional isogeny $\Phi_{u_2}$ of polarised degree $u_2$, which we split again (via a 4-dimensional isogeny of polarised degree $2^f$) to get $\Phi_u$. $\square$

*Remark 4.2.*
- To compute $\Gamma$ as defined in Eq. (10), we only need the $m$-torsion to be smooth and accessible on $E$, where $N_1 + N_2 | m^2$, see [58, § 6.3 and Appendix B]. This allows to relax the conditions on Theorem 4.1 to $u < 2^{2e}$ and $2^{3e-1} \gg \Delta$.

- In practice, for our application from Section 3, we have $u, v \approx \sqrt{|\Delta|}$. Furthermore, in the case of CSIDH, we have $\Delta = -p$, and we select $p$ such that $p = c2^f - 1$ for a small cofactor $c$. In particular, we can take $e = f$ in Theorem 4.1, so that $2^e \approx p$. In that case we already have $u(2^e - u) \gg p$, and we do not need the two step solution, nor the more relaxed torsion condition from the item above.

Combining Section 3 with Theorem 4.1, the resulting algorithm to compute the action of a quadratic ideal $\mathfrak{a} \subset \mathcal{O}$ on an $\mathcal{O}$-oriented elliptic curve is very similar to the algorithm used in [5] to convert a quaternionic ideal to an isogeny when the full endomorphism ring of $E$ is known. We make up for the lack of known endomorphisms on $E$ to build isogenies of suitable fixed degree by going up to dimension 2 and working on $E^2$. This doubling of dimension propagates to the rest of the algorithm, and we end up computing $\Phi_u, \Phi_v$ and $E_{\mathfrak{a}}$ via an isogeny of dimension 4, compared to dimension 2 in [5].

This change of dimension has important consequences for the fine-tuning of the algorithm. Indeed, as we have seen in Section 4.1, $\Phi_u$ (resp. $\Phi_v$) are much easier to compute whenever $u$ (resp. $v$) is a sum of two squares, or at least $\mathfrak{B}$-good in the sense of Definition 4.1. In particular, in that case we do not need to use a 4-dimensional isogeny to compute $\Phi_u$.

A similar situation happened in [5]: if $u$ was a sum of two squares we could compute $\Phi_u$ directly without going through a 2-dimensional isogeny. However, experiments showed that spending more time on the norm equation to find $u$ a sum of two squares did not compensate having to do one less 2-dimensional isogeny. The situation is different in our case, because a 4-dimensional $2^e$-isogeny is much slower than a 2-dimensional $2^e$-isogeny (by a factor roughly $\times 8$). Depending on the size of our parameters, experiments show that the optimal choices depend on the size of $\Delta$.

- If $\Delta$ is of moderate size, we look for $\mathfrak{B}$-good $u, v$, and only compute one 4-dimensional isogeny for $E_{\mathfrak{a}}$.

- If $\Delta$ is of large size, we look for only one of $u, v$ to be $\mathfrak{B}$-good, while we let the other one be arbitrary. This time computing $\varphi_{\mathfrak{a}}$ requires two 4-dimensional isogenies. For CSIDH, our experiments in Section 6.4 suggest this strategy becomes faster than the previous one for $p$ of 4000 bits. We remark that for very large $\Delta$, we expect that letting $u, v$ be arbitrary and computing three 4-dimensional isogenies would be the optimal choice.

In the case of CSIDH (detailed in Sections 5 and 6), only the first choice is implemented and the computation stays reasonably efficient even for our biggest parameters ($p$ of 4000 bits) as long as the allowed set of small split primes $\mathfrak{B}$ is well chosen. We refer to Section 6.4 for further discussion on the above trade-off.
## 5 Our algorithm specialized to CSIDH

In this section, we apply our algorithm to the CSIDH setting, or more accurately, the CSURF setting [12]; we fix a prime of the form $p = c2^f - 1$ where $c$ is a small odd integer and $f \gg 1$ so that $p \equiv 7 \pmod 8$, and work with curves oriented by $\mathbb{Z}[\frac{1+\sqrt{-p}}{2}]$, where the orientation is given by sending $\sqrt{-p}$ to the Frobenius endomorphism. Not only does using the orientation given by Frobenius make evaluating the orientation very fast, but by using the maximal order instead of the typical CSIDH suborder $\mathbb{Z}[\sqrt{-p}]$, we are able to execute our whole algorithm over $\mathbb{F}_p$, by choosing a particularly suitable basis of $E[2^e]$, and working with $x$-only arithmetic. This is further detailed in Section 5.2 and Appendix D.

At a high level, the three-step process is summarized as in Algorithm 1. We now describe the involved steps in more detail.

<caption>**Algorithm 1 EVALUATEACTION($E, \mathfrak{a}$)**</caption>

**Input:** $E$ a supersingular elliptic curve over $\mathbb{F}_p$ defined on the surface.
**Input:** $\mathfrak{a} \subseteq \mathbb{Z}[\frac{1+\sqrt{-p}}{2}]$ an ideal.
**Output:** $\mathfrak{a} \star E$
1: Compute $\mathfrak{l} \subseteq \mathbb{Z}[\frac{1+\sqrt{-p}}{2}]$ an ideal of small norm. ▷ Precomputation
2: Set $E' := E$ and $\mathfrak{a}' = \mathfrak{a}$.
3: **fail**, $\mathfrak{b}_e, \mathfrak{b}_k, \mathfrak{c}_e, \mathfrak{c}_k, u, v, e := \text{FINDUV}(\mathfrak{a}')$ ▷ Step 1
4: **if fail then**
5: Compute $\mathfrak{a}' = \mathfrak{a}'\mathfrak{l}$ ▷ Rerandomize the ideal
6: Compute $E' = \bar{\mathfrak{l}} \star E'$ using Elkies algorithm.
7: **run** Line 3 on $E', \mathfrak{a}'$.
8: **end if**
9: **KernelData** := COMPUTE**KernelData**($E', \mathfrak{b}_e, \mathfrak{b}_k, \mathfrak{c}_e, \mathfrak{c}_k, u, v, e$) ▷ Step 2
10: Compute the 4-dimensional isogeny $F$, defined by **KernelData**. ▷ Step 3
11: Extract $E_\mathfrak{a}$ from the codomain of $F$. ▷ See Appendix B.4.
12: **return** $E_\mathfrak{a}$.

### 5.1 Step 1: the norm equation

In this section, we explain how to solve the norm equation (4) sketched in Section 3.2. This subsection is summarized in Algorithm 2. As in Section 4.1, we fix a set of small Elkies primes $\mathfrak{B}$, i.e. of primes that split in $\mathbb{Q}(\sqrt{-p})$. Notice that we can always impose $2 \in \mathfrak{B}$. We begin by creating a list of primitive ideals equivalent to $\mathfrak{a}$. More precisely, we obtain a list of equivalent ideals $\mathfrak{a}^{(i)}$ such that $\mathfrak{a}^{(i)} = \mathfrak{a}_e^{(i)}\mathfrak{a}_k^{(i)}$ and $N(\mathfrak{a}_k^{(i)})$ is a product of primes in $\mathfrak{B}$. We then want to find a pair $(i, j)$ such that for $N_1 = N(\mathfrak{a}_k^{(i)})$ and $N_2 = N(\mathfrak{a}_k^{(j)})$ (note that the smooth part has been removed from the norm), we can find positive integers $u, v$ solving the equation

$$uN_1 + vN_2 = 2^e. \qquad (11)$$
Since generically $N_1N_2$ determines the minimal value for $e$ such that the above equation has a solution, we test pairs in order of increasing $N_1N_2$.

Given such a pair of norms $N_1, N_2$, we first check that they are coprime. Indeed, if $\gcd(N_1, N_2)$ is not a power of 2, Equation (11) has no solutions. Since $2 \in \mathfrak{B}$ and we assume $N_1, N_2$ has no prime factors in $\mathfrak{B}$, this requirement simply becomes $\gcd(N_1, N_2) = 1$.

If $\gcd(N_1, N_2) = 1$, we can find $u_0$ and $v_0$ such that $u_0N_1 + v_0N_2 = 1$ via the extended Euclidean algorithm, and consequently $2^e u_0N_1 + 2^e v_0N_2 = 2^e$. Note that one of $u_0$ and $v_0$ will be negative and that we parametrize all possible solutions to the equation by writing

$$ (2^e u_0 + kN_2)N_1 + (2^e v_0 - kN_1)N_2 = 2^e, $$

for some $k \in \mathbb{Z}$. This also allows to determine the minimal $e$ for which Equation (11) has a solution by enumerating the set of $k$ for which both $u = 2^e u_0 + kN_2$ and $v = 2^e v_0 - kN_1$ are positive.

**Remark 5.1.** The exponent $e$ in Equation (11) determines the number of four-dimensional 2-isogeny steps to be computed, so taking a lower exponent $e$ will be more efficient.

As remarked before and also observed in [5], solving Equation (11) directly for $\mathfrak{a}^{(i)}$ and $\mathfrak{a}^{(j)}$, often fails since $N_1, N_2$ are in $O(\sqrt{p})$ and $2^e < p$, which makes the original Clapoti approach unpractical in our case. Removing factors in $\mathfrak{B}$, makes $N_1$ and $N_2$ much smaller and hence greatly increases the probability of finding solutions to Equation (11): the number of expected solutions grows roughly as the product of all factors removed. This leads to the following tradeoff: the bigger the set of primes $\mathfrak{B}$, the higher the probability of finding a solution, but the more time is spent in the second step computing small degree isogenies. Furthermore, even if the probability of finding a solution for a given choice of $\mathfrak{B}$ is low, and in particular, we fail to find a solution, we can simply rerandomize the starting ideal (by some small ideal) and start over. More details about this approach are given below.

Assuming that we can find enough solution pairs $(u, v)$, for each one we have to determine whether $u$ and $v$ are $\mathfrak{B}$-good, i.e. if we can write $u = g_u(x_u^2 + y_u^2)$ and $v = g_v(x_v^2 + y_v^2)$ with the prime factors of $g_u$ and $g_v$ in $\mathfrak{B}$. Recall that a number can be expressed as a sum of two squares if and only if all its prime factors that are 3 mod 4 appear with an even exponent. The obvious choice is then to set $g_u$ (resp. $g_v$) to be the product of all prime factors that are 3 mod 4 and contained in $\mathfrak{B}$ and dividing $u$ (resp. $v$). To determine whether what is left contains more prime factors 3 mod 4 we would need to factor $u$, which is too expensive in general. As such we simply proceed by trial division up to a small bound, hoping that what is left at the end is a prime that is 1 mod 4 (note that we simply reject when it is 3 mod 4). If everything succeeds, we have found the factorization of $u$, and hence its decomposition as a sum of two squares. The same procedure applies to $v$. How to choose an optimal set $\mathfrak{B}$ is analyzed in greater detail in Section 6.1.
Once we have found a full solution, we proceed to the next step: we will denote the ideals $\mathfrak{a}^{(i)}$, $\mathfrak{a}^{(j)}$ as $\mathfrak{b}$ and $\mathfrak{c}$, so in particular, $N_1 = N(\mathfrak{b}_k)$ and $N_2 = N(\mathfrak{c}_k)$.

**Rerandomization** For some ideals, Equation (11) may not have a solution with $u, v$ of the required shape. As already observed, increasing the bound $\mathfrak{B}$ is always a possible solution, though this quickly becomes expensive.

Another solution to this is to simply rerandomize the starting ideal $\mathfrak{a}$. Namely, if none of the sampled, equivalent ideals gives us a solution for Equation (11), we can multiply $\mathfrak{a}$ by a non-principal ideal $\mathfrak{l}$ for which the action is easy to evaluate (e.g. an ideal above $\mathfrak{l}$, where $\mathfrak{l}$ is the smallest split prime). We obtain a new ideal $\mathfrak{a}' = \mathfrak{a}\mathfrak{l}$ and try to solve Equation (11) for this ideal instead. If we manage to do so, we can simply run our algorithm on $\mathfrak{a}'$, using this solution, starting from the starting curve $E_{\mathfrak{l}} = \mathfrak{l} \star E$ instead, which is easy to compute by the choice of $\mathfrak{l}$. Otherwise, we can rerandomize again by computing $\mathfrak{a}'' = \mathfrak{l}\mathfrak{a}'$, and so on. As soon as Equation (11) has a reasonable chance of being solved, this method takes care of failures quite efficiently, making Step 1 faster at the cost of a few Elkies steps. The average number of rerandomizations for different levels are reported in Table 2.

<caption>**Algorithm 2 FINDUV(\mathfrak{a})**</caption>

**Input:** $\mathfrak{a} \subseteq \mathbb{Z}^{\lceil \frac{1+\sqrt{-p}}{2} \rceil}$ an ideal.
**Output:** Fail, flag indicating failure to find solution.
**Output:** $\mathfrak{b}_e, \mathfrak{b}_k, \mathfrak{c}_e, \mathfrak{c}_k \subseteq \mathbb{Z}^{\lceil \frac{1+\sqrt{-p}}{2} \rceil}$ ideals, such that $\mathfrak{b}_e \mathfrak{b}_k \sim \mathfrak{c}_e \mathfrak{c}_k \sim \mathfrak{a}$.
**Output:** $u, v, e$ positive integers, such that $u \cdot N(\mathfrak{b}_k) + v \cdot N(\mathfrak{c}_k) = 2^e$, for $e < f - 3$.
1: Compute a list of short ideals $\mathfrak{a}^{(i)}$ equivalent to $\mathfrak{a}$
2: Write $\mathfrak{a}^{(i)} = \mathfrak{a}_e^{(i)} \mathfrak{a}_k^{(i)}$, such that $N(\mathfrak{a}_e^{(i)})$ is a product of primes in $\mathfrak{B}$
3: Sort the list by $N(\mathfrak{a}_k^{(i)})$
4: **for** all pairs $(i, j)$ with $1 \le i < j \le 10$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;Set $N_1 = N(\mathfrak{a}_k^{(i)})$, $N_2 = N(\mathfrak{a}_k^{(j)})$
6: &nbsp;&nbsp;&nbsp;&nbsp;Find all pairs $(u, v)$ such that $uN_1 + vN_2 = 2^e$, with $e \le f - 3$
7: &nbsp;&nbsp;&nbsp;&nbsp;**for** each pair $(u, v)$ **do**
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Remove factors of $\mathfrak{B}$ from $u, v$
9: &nbsp;&nbsp;&nbsp;&nbsp;Perform trial division on $u, v$ ▷ Bound depends on security level
10: &nbsp;&nbsp;&nbsp;&nbsp;**if** factors 3 mod 4 are found with odd exponent **then**
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Continue to next pair
12: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
13: &nbsp;&nbsp;&nbsp;&nbsp;**if** remaining part of $u, v$ are both prime **then**
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**return** False, $\mathfrak{a}_e^{(i)}, \mathfrak{a}_k^{(i)}, \mathfrak{a}_e^{(j)}, \mathfrak{a}_k^{(j)}, u, v, e$
15: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
16: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
17: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
18: **return** True, $-, -, -, -, -, -, -$
## 5.2 Step 2: computing kernel points

After determining a suitable pair of ideals $b, c$, our goal is to construct the 4-dimensional $2^e$-isogeny described in Section 3.2. We now detail how to derive the kernel of this isogeny. This subsection is summarized in Algorithm 3.

Let $E$ denote the starting curve on which we want to act with $\mathfrak{a}$, then we obtain the following diagram derived from the isogeny diamond in Section 3.2:

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

In the diagram above, we have:

- $\Psi_1 = \mathrm{diag}(\varphi_{b_e}, \varphi_{b_e})$, where $\varphi_{b_e}$ is the isogeny corresponding to the action of $b_e$; analogously, $\Psi_2 = \mathrm{diag}(\varphi_{c_e}, \varphi_{c_e})$;

- $\Phi_1 = \mathrm{diag}(\varphi_{b_k}, \varphi_{b_k})$ and $\Phi_2 = \mathrm{diag}(\varphi_{c_k}, \varphi_{c_k})$;

- $\Phi_u$ and $\Phi_v$ are isogenies of polarised degree $u$ and $v$ respectively; as explained in Section 4.1, we can write $\Phi_u = M_u \mathrm{diag}(\varphi_u, \varphi_u)$, with $\mathrm{deg}(\varphi_u) = g_u$, $\Phi_v = M_v \mathrm{diag}(\varphi_v, \varphi_v)$ with $\mathrm{deg}(\varphi_v) = g_v$, $M_u$ and $M_v$ being given by Equation (9).

Further, the isogenies $\Phi$ and $\Psi$ are the isogenies of polarised degree $uN_1$ and $vN_2$ completing the square, whose existence is guaranteed by Lemma 2.2. In this setting we can also describe them explicitly, as the following lemma shows.

**Lemma 5.1.** $\Phi$ is of the form $\Phi'_1 \circ \tilde{\Phi}'_u$, where:

- $\Phi'_1 := \mathrm{diag}(\varphi'_{b_k}, \varphi'_{b_k}) : E'_1^2 \to E_v^2$ with $E'_1 := [\mathfrak{b}_k] \cdot E_v$ and $\varphi'_{b_k}$, the isogeny associated to the action of $b_k$ on $E'_1$;

- $\Phi'_u := M_u \mathrm{diag}(\varphi'_u, \varphi'_u) : E'_1^2 \to E'^2$, $M_u$ being given by Equation (9) and $\varphi'_u$ being the isogeny of degree $g_u$ given by the action of the same ideal $I_u$ as $\varphi_u$ (as a product of Elkies primes).

$\Psi$ is of the form $\Phi'_v \circ \tilde{\Phi}'_2$, where:

- $\Phi'_2 := \mathrm{diag}(\varphi'_{c_k}, \varphi'_{c_k}) : E'_2^2 \to E_u^2$ with $E'_2 := [\mathfrak{c}_k] \cdot E_u$ and $\varphi'_{c_k}$, the isogeny associated to the action of $c_k$ on $E'_2$;
– $\Phi'_v := M_v \operatorname{diag}(\varphi'_v, \varphi'_v) : E'_2^2 \to E'^2$, $M_v$ being given by Equation (9) and $\varphi'_v$ being the isogeny of degree $g_v$ given by the action of the same ideal $I_v$ as $\varphi_v$ (as a product of Elkies primes).

In particular, the common codomain of $\tilde{\Phi}$ and $\Psi$ is a product of elliptic curves $E'^2$.

*Proof.* We have to verify that $\Phi'_1 \circ \tilde{\Phi}'_v$ and $\Phi'_v \circ \tilde{\Phi}'_2$ defined above are $N_1 u$ and $N_2 v$-isogenies respectively making the diagram commute, i.e. such that:

$$ \Phi'_1 \circ \tilde{\Phi}'_v \circ \Phi'_v \circ \tilde{\Phi}'_2 = \Phi_v \circ \tilde{\Phi}_2 \circ \Phi_1 \circ \tilde{\Phi}_u. \qquad (12) $$

First, we verify that the composition on the left makes sense i.e. that $\Phi'_u$ and $\Phi'_v$ have the same codomain. By definition, the codomain of $\Phi'_u$ is

$$ E' = [I_u]E'_1 = [I_u \bar{b}_k]E_v = [I_u \bar{b}_k I_v]E_2 = [I_u I_v \bar{b}_k c_e]E $$

and the codomain of $\Phi'_v$ is

$$ E'' = [I_v]E'_2 = [I_v \bar{c}_k]E_u = [I_v \bar{c}_k I_u]E_1 = [I_u I_v \bar{c}_k b_e]E. $$

But $b = b_e b_k$ and $c = c_e c_k$ are both equivalent to $a$. It follows that $\bar{b}_k c_e$ and $\bar{c}_k b_e$ are equivalent, so that $E' \simeq E''$.

By construction, $\Phi'_1 \circ \tilde{\Phi}'_v$ and $\Phi'_v \circ \tilde{\Phi}'_2$ are $N_1 u$ and $N_2 v$-isogenies respectively, so we only have to prove Equation (12). On the one hand, we have:

$$ \Phi'_1 \circ \tilde{\Phi}'_v \circ \Phi'_v \circ \tilde{\Phi}'_2 = \tilde{M}_u \circ M_v \circ \operatorname{diag}(\varphi'_{b_k} \circ \hat{\varphi}'_u \circ \varphi'_v \circ \hat{\varphi}'_{c_k}, \varphi'_{b_k} \circ \hat{\varphi}'_u \circ \varphi'_v \circ \hat{\varphi}'_{c_k}), $$

since $\tilde{M}_u$ and $M_v$ commute with diagonal isogenies. On the other hand:

$$ \Phi_v \circ \tilde{\Phi}_2 \circ \Phi_1 \circ \tilde{\Phi}_u = M_v \circ \tilde{M}_u \circ \operatorname{diag}(\varphi_v \circ \hat{\varphi}_{c_k} \circ \varphi_{b_k} \circ \hat{\varphi}_u, \varphi_v \circ \hat{\varphi}_{c_k} \circ \varphi_{b_k} \circ \hat{\varphi}_u). $$

Both $\varphi'_{b_k} \circ \hat{\varphi}'_u \circ \varphi'_v \circ \hat{\varphi}'_{c_k}$ and $\varphi_v \circ \hat{\varphi}_{c_k} \circ \varphi_{b_k} \circ \hat{\varphi}_u$ correspond to the action of the ideal $\bar{I}_u I_v \bar{c}_k b_k$ on $E_u$ so these isogenies must be equal. Besides, a simple matrix computation ensures that $\tilde{M}_u$ and $M_v$ commute. This proves Equation (12) and the lemma. $\square$

Let $(P, Q)$ be a basis of $E_1[2^{e+2}]$ (recall from Section 3.1 that the theta isogeny algorithm requires the $2^{e+2}$-torsion to compute a $2^e$-isogeny, see Appendix A for more details). To compute the kernel of the isogeny in dimension 4, we need to compute the points

$$ P_u = \varphi_u(P), \quad Q_u = \varphi_u(Q) $$

on $E_u$ and

$$ P_v = \varphi_v \circ \hat{\varphi}_{c_k} \circ \varphi_{b_k}(P), \quad Q_v = \varphi_v \circ \hat{\varphi}_{c_k} \circ \varphi_{b_k}(Q), $$

on $E_v$.
**Isogenies of degree $g_u$ and $g_v$.** The isogenies $\varphi_u$ and $\varphi_v$ with $\deg(\varphi_u) = g_u$ and $\deg(\varphi_v) = g_v$ have no specific requirement other than their degree. Although it is possible to deal with all small primes via Vélu's algorithm [64], it is more efficient to restrict $\mathfrak{B}$ to Elkies primes (i.e. splitting primes in $\mathbb{Q}(\sqrt{-p})$) and we can use the corresponding algorithm [30] to compute the small degree isogenies.

Note that $\varphi_u$ of degree $g_u$ and $\varphi_v$ of degree $g_v$ can be used to construct 2-dimensional isogenies of polarised degree $u = g_u(x_u^2 + y_u^2)$ and $v = g_v(x_v^2 + y_v^2)$ respectively by composition with matrices $M_u, M_v$ described in Section 4.1. However, we do not need to perform that step now, as it can be included in the 4-dimensional computation.

**Working over $\mathbb{F}_p$.** Note that, by definition, the isogenies $\varphi_b$ and $\varphi_c$ are $\mathbb{F}_p$-rational, and as noted above, our choice of primes in $\mathfrak{B}$ imply that $\varphi_u$ and $\varphi_v$ can also be chosen to be defined over $\mathbb{F}_p$. The main difficulty comes from evaluating these isogenies on a basis of $E[2^{e+2}]$, which is not defined over $\mathbb{F}_p$.

We will show that all operations can in fact be carried out over the Kummer line, where we have a rational basis. Let $f$ the largest integer so that $2^f \mid (p+1)$. By our choice of prime, the two eigenvalues of Frobenius on $E[2^{f-1}]$ are $\pm 1$. Thus, we can pick a basis $\langle P, Q \rangle = E[2^{f-1}]$ corresponding to the two eigenvectors of Frobenius; we let $P, Q$ be points of order $2^{f-1}$ satisfying

$$ \pi(P) = P, \quad \pi(Q) = -Q. \qquad (13) $$

Clearly, these points $P, Q$ are $\mathbb{F}_p$-rational on the Kummer line $E/\pm 1$, and together (adjusting by translation by a point of 2-torsion if needed) they form a basis of $E[2^{f-1}]$. This is where we use that we are on the surface. If we were on the floor, picking similar points would not give a basis, since in that case, $[2^{f-2}]P = [2^{f-2}]Q$, the unique rational point of 2-torsion.

In Appendix D.1 we show how to efficiently generate such a basis, based on a technique which is similar to how one usually computes bases of $E[2^f]$ over $\mathbb{F}_{p^2}$.

**Evaluating the ideals.** The rest of this section will be devoted to the computation of $P_v, Q_v$. In particular, we need to evaluate $\hat{\varphi}_{c_k} \circ \varphi_{b_k}$ on a basis of $E_1[2^{e+2}]$. From Equation (7) we know that

$$ \hat{\varphi}_{c_k} \circ \varphi_{b_k} = \frac{1}{N(\mathfrak{b}_e)N(\mathfrak{c}_e)} \varphi_{c_e} \circ \iota(\sigma) \circ \hat{\varphi}_{\mathfrak{b}_e} \qquad (14) $$

with $\overline{cb} = \sigma O$. By construction, $N(\mathfrak{b}_e)$ and $N(\mathfrak{c}_e)$ are both smooth, with all their factors in $\mathfrak{B}$. Let us assume for simplicity that they are also odd; the even case needs a bit more care and is detailed afterwards.

We detail how to evaluate the endomorphism $\iota(\sigma)$. First, assume that $\sigma$ is of the form $a+b\sqrt{-p}$, with $a, b \in \mathbb{Z}$. Explicitly, this is mapped to the endomorphism $[a] + [b]\pi$, where $\pi$ denotes the Frobenius endomorphism. Evaluating such an
endomorphism on our choice of basis is particularly easy: with the same notation as in the previous section, the evaluation is simply given by

$$[a]P + [b]\pi(P) = [a]P + [b]P = [a+b]P,$$
$$[a]Q + [b]\pi(Q) = [a]Q + [b](-Q) = [a-b]Q.$$

Note in particular that we can evaluate this on the Kummer line. Although $\sigma$ is an element of the order $\mathbb{Z}[\frac{1+\sqrt{-p}}{2}]$, and hence is more generally of the form $a + b\sqrt{-p}$, where $a, b$ are in $\mathbb{Z}[\frac{1}{2}]$, i.e. are rational numbers with denominator at most 2, we show in Appendix D.2 how to evaluate this on our basis (on the Kummer line) without having to use division points, which may not be $\mathbb{F}_p$-rational.

Further, note that since we are after the action on the $2^{e+2}$-torsion, we can compute $(N(\mathfrak{b}_e)N(\mathfrak{c}_e))^{-1}$ (mod $2^{e+2}$), to account for the division in Equation (14). Hence, all that remains to evaluate $\hat{\phi}_{\mathfrak{c}_k} \circ \phi_{\mathfrak{b}_k}$ is to describe how we efficiently recover the isogenies $\phi_{\mathfrak{b}_e}$ and $\phi_{\mathfrak{c}_e}$.

These isogenies are (by definition) the product of prime degree Elkies isogenies. More specifically, say that for a given prime $\ell$, we have $v_\ell(N(\mathfrak{b}_e)) = k$, i.e. $\ell$ divides $N(\mathfrak{b}_e)$ exactly $k$ times. Then $\mathfrak{b}_e$ will act with exactly $k$ degree-$\ell$ Elkies isogenies. We can also assume that all these isogenies will be in the same direction, since otherwise $\mathfrak{b}_e$ and consequently $\mathfrak{b}$ would factor through multiplication by $\ell$, and we assumed that $\mathfrak{b}$ was primitive (if it was not, the ideal $\mathfrak{b}/\ell$ would also be equivalent to $\mathfrak{a}$, but with smaller norm).

With Elkies' algorithm, we can return both kernel polynomials of degree $\ell$ determining the two Elkies isogenies of degree $\ell$. To check which polynomial corresponds to the correct isogeny, we check the eigenvalue of Frobenius, as described in the original CSIDH paper [13, Section 3]. In more detail, assume that $(\ell, \sqrt{-p-\lambda}) = \mathfrak{l} \subseteq \mathfrak{b}_e$, i.e. we wish to evaluate the action of $\mathfrak{l}$. Since the roots of the correct kernel polynomial $h(x)$ are precisely the $x$-coordinates of the points of order $\ell$, satisfying $\pi(P) = [\lambda]P$, we simply test this equality symbolically: we have the equality

$$(x^p, y^p) \equiv [\lambda](x, y) \pmod{h(x), y^2 - f(x)},$$

where $f(x)$ is the polynomial coming from the curve equation. If these values are not equal, we know that this kernel polynomial corresponds to the other eigenvalue. Note that when we are doing multiple $\ell$-isogenies, we need only check the eigenvalue of Frobenius for the first $\ell$-isogeny; after that, we can recognise the right kernel polynomial to be the one not corresponding to the previous $j$-invariant.

**Remark 5.2.** If $\mathfrak{b}_e$ and $\mathfrak{c}_e$ share some steps, we can perform those steps (corresponding to the action of the ideal $\mathfrak{b}_e + \mathfrak{c}_e$) in advance, and start evaluating the ideal action from the resulting curve. The advantage is twofold: first, we avoid evaluating the same ideal twice. Second, we do not need to evaluate the common part on points. We can instead start computing $P_u$ and $P_v$ directly from the codomain curve.
**Even norm part** We now explain how to deal with even normed ideals. For this, we write the ideal $\mathfrak{b}_e$ as

$$ \mathfrak{b}_e = \mathfrak{b}_{e'} \mathfrak{b}_2 $$

where $\mathfrak{b}_2$ is an ideal of norm $2^{f_b}$, and correspondingly for $\mathfrak{c}_e$. As explained in Remark 5.2, the shared steps of $\mathfrak{b}_e, \mathfrak{c}_e$ can be computed at the beginning, hence we can assume that $\mathfrak{b}_2$ and $\mathfrak{c}_2$ are coprime (i.e. the corresponding 2-isogenies go in “opposite” directions). Writing $f_{12} = f_b + f_c$, we get that

$$ \hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k} = \frac{1}{N(\mathfrak{b}_{e'})N(\mathfrak{c}_{e'})2^{f_{12}}} \varphi_{\mathfrak{c}_e} \circ \iota(\sigma) \circ \hat{\varphi}_{\mathfrak{b}_e}. $$

The only added difficulty that needs to be accounted for is that we cannot find an inverse to $2^{f_{12}}$ modulo $2^{e+2}$.

We fix this by taking a basis of a larger torsion group. Using the same technique as for the odd normed case, we can thus evaluate

$$ [2^{f_{12}+1}] \circ \hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k} = \frac{1}{N(\mathfrak{b}_{e'})N(\mathfrak{c}_{e'})} \varphi_{\mathfrak{c}_e} \circ \iota(\sigma) \circ \hat{\varphi}_{\mathfrak{b}_e}, $$

and hence we can find the evaluation of $\hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k}$ on $E_1[2^{e+2}]$ by instead starting with a basis of $E[2^{e+2+f_{12}}]$. (We remark that by construction, we have $e+2+f_{12} \le f$.)

**Algorithm 3 COMPUTE KERNELDATA(E, $\mathfrak{b}_e$, $\mathfrak{c}_e$, $\mathfrak{b}_k$, $\mathfrak{c}_k$, u, v, e)**

**Input:** E a supersingular elliptic curve over $\mathbb{F}_p$ defined on the surface.

**Input:** $\mathfrak{b}_e, \mathfrak{b}_k, \mathfrak{c}_e, \mathfrak{c}_k \subseteq \mathbb{Z}^{[1+\sqrt{-p}/2]}$ ideals, such that $\mathfrak{b}_e \mathfrak{b}_k \sim \mathfrak{c}_e \mathfrak{c}_k$.

**Input:** $u, v, e$ positive integers, such that $u \cdot N(\mathfrak{b}_k) + v \cdot N(\mathfrak{c}_k) = 2^e$, for $e < f-3$.

**Output:** KernelData necessary for computing the 4-dimensional isogeny.

1: Compute a basis $P, Q$ of $E[2^{f-1}]$. ▷ See Appendix D to sample special basis.

2: Compute the isogeny $\phi_{\mathfrak{b}_e} : E \to E_1$ using Elkies algorithm.

3: $P_1, Q_1 := \phi_{\mathfrak{b}_e}(P), \phi_{\mathfrak{b}_e}(Q)$.

4: Double $P_1, Q_1$ until they are of order $2^{e+2}$. ▷ Nr. of doublings depends on $\phi_{\mathfrak{b}_e}$.

5: Write $u = g_u(x_u^2 + y_u^2)$ and $v = g_v(x_v^2 + y_v^2)$.

6: Compute a random isogeny $\phi_u : E_1 \to E_u$ of degree $g_u$ using Elkies algorithm.

7: $P_u, Q_u := \phi_u(P_1), \phi_u(Q_1)$.

8: Find a generator $\sigma$ of the principal ideal $\overline{\mathfrak{c}_e \mathfrak{c}_k \mathfrak{b}_e \mathfrak{b}_k}$.

9: Compute the isogeny $\phi_{\mathfrak{c}_e} : E \to E_2$, using Elkies algorithm.

10: Compute $d \equiv N(\mathfrak{c}_{e'})^{-1} (\text{mod } 2^{e+2})$, where $\mathfrak{c}_{e'}$ is the odd normed part of $\mathfrak{c}_e$.

11: $P_2, Q_2 := [d] \circ \phi_{\mathfrak{c}_e} \circ \sigma(P), [d] \circ \phi_{\mathfrak{c}_e} \circ \sigma(Q)$. ▷ See Algorithm 4 for evaluating $\sigma$.

12: Double $P_2, Q_2$ until they are of order $2^{e+2}$. ▷ Nr. of doublings depends on $\phi_{\mathfrak{b}_e}$.

13: Compute a random isogeny $\phi_v : E_2 \to E_v$ of degree $g_v$ using Elkies algorithm.

14: Compute $P_v, Q_v = \phi_v(P_2), \phi_v(Q_2)$.

15: KernelData := $[N(\mathfrak{b}_k), N(\mathfrak{c}_k), g_u, x_u, y_u, g_y, x_v, y_v, P_u, Q_u, P_v, Q_v]$.

16: return KernelData.
### 5.3 Step 3: computing the 4-dimensional isogeny

In this step, we compute the 4-dimensional $2^e$-isogeny isogeny $F : A_u \times A_v \to E_a^2 \times A$ from Section 3.2 (defined in Equation (5)). We can actually be much more explicit than in Section 3.2. From Section 5.2, we obtain that $A_u = E_u^2$, $A_v = E_v^2$ and $A = E'^2$. In particular

$$F = \begin{pmatrix} \Phi_1 \circ \tilde{\Phi}_u & \Phi_2 \circ \tilde{\Phi}_v \\ -\tilde{\Phi}'_v \circ \tilde{\Phi}'_2 & \tilde{\Phi}'_u \circ \tilde{\Phi}'_1 \end{pmatrix} : E_u^2 \times E_v^2 \to E_a^2 \times E'^2.$$

Equation (6) also becomes:

$$ \begin{aligned} \ker(F) &= \left\{ \begin{aligned} &([uN_1]\varphi_u(P), [uN_1]\varphi_u(Q), \\ &\quad [g_u]M_v \cdot M_u^T(\varphi_v \circ \hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k}(P), \varphi_v \circ \hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k}(Q))) \mid P, Q \in E_u[2^e] \end{aligned} \right\} \\ &= \left\{ \begin{aligned} &([N_1]M_u(\varphi_u(P), \varphi_u(Q)), \\ &\quad M_v(\varphi_v \circ \hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k}(P), \varphi_v \circ \hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k}(Q))) \mid P, Q \in E_1[2^e] \end{aligned} \right\}. \end{aligned} \quad (15) $$

To compute $F$, we proceed with theta coordinates of level 2 and the algorithms from [18] (see Appendix B). These algorithms require a basis of an isotropic subgroup (for the Weil pairing) $K' \subset (E_u^2 \times E_v^2)[2^{e+2}]$ such that $[4]K' = \ker(F)$.

Let $(P, Q)$ be a basis of $E_1[2^{e+2}]$, $\zeta := e_{2^{e+2}}(P, Q)$, $\eta_1, \alpha$ be inverses of $N_1, uN_1$ modulo $2^{e+2}$ respectively, $(P_u, Q_u) := \varphi_u(P, Q)$ and $(P_v, Q_v) := \varphi_v \circ \hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k}(P, Q)$ the points computed in Step 2. Then, applying Lemma A.1.(i) to $F$ and the $\zeta$-symplectic basis $(x_1, x_2, y_1, y_2) := ((\varphi_{\mathfrak{b}_k}(P), 0), (0, \varphi_{\mathfrak{b}_k}(P)), ([\eta_1]\varphi_{\mathfrak{b}_k}(Q), 0), (0, [\eta_1]\varphi_{\mathfrak{b}_k}(Q)))$, we obtain a basis

$$T_1 := ([N_1x_u]P_u, [N_1y_u]P_u, [x_v]P_v, [y_v]P_v)$$

$$T_2 := ([-N_1y_u]P_u, [N_1x_u]P_u, [-y_v]P_v, [x_v]P_v)$$

$$T_3 := ((1 - 2^e \alpha)x_u)Q_u, ((1 - 2^e \alpha)y_u)Q_u, [\eta_1 x_v]Q_v, [\eta_1 y_v]Q_v)$$

$$T_4 := ([-(1 - 2^e \alpha)y_u]Q_u, [(1 - 2^e \alpha)x_u]Q_u, [-\eta_1 y_v]Q_v, [\eta_1 x_v]Q_v),$$

of an isotropic subgroup $K' \subset (E_u^2 \times E_v^2)[2^{e+2}]$ such that $[4]K' = \ker(F)$.

To compute $F$ and extract $E_a$ from its codomain, we take as input:

- The points $P_u, Q_u, P_v, Q_v$ computed in Step 2.

- The integers $N_1, N_2, g_u, x_u, y_u, g_v, x_v, y_v, e$ obtained from a solution of Equation (11) obtained in Step 1.

This data is not only useful to compute $T_1, \dots, T_4$ but also gluing and diagonal isogenies located in the beginning of the 2-isogeny chain $F$ that must be handled with care. For this technical reason, one must provide the data above and not
$T_1, \cdots, T_4$ or kernel points of $F$ directly. We refer to Appendix B for a detailed presentation of the computation of $F$ with level-2 theta coordinates (following the approach of [18]) and to Appendix A for a cryptographer friendly introduction to higher-dimensional isogeny computations with theta coordinates.

**Remark 5.3 (Working over $\mathbb{F}_p$).** Using the basis from Equation (13) for our elliptic curves, we can work with rational theta null points all along our dimension 4 isogeny chain. This is one of the reasons to work on the surface: on the floor, the elliptic curve level 2 theta null points are not rational. Our kernel generators have Frobenius-eigenvalue $\pm 1$, so define $\mathbb{F}_p$-rational points on the Kummer variety and thus their theta-coordinates are $\mathbb{F}_p$-rational (recall that level-2 theta coordinates only give an embedding of the Kummer into projective space, not of the whole variety). This allows us to do the full dimension-4 isogeny computation over $\mathbb{F}_p$, with the exception of the gluing isogenies. Indeed, as currently implemented the gluing formulas mix points with eigenvalue 1 and points with eigenvalue $-1$, so we need to temporally work over $\mathbb{F}_{p^2}$. We expect that using more symmetric gluing formulas would allow us to stay over $\mathbb{F}_p$ even for gluing, but we leave that for future work.

# 6 Implementation results

In this section we describe different implementation choices and their relative efficiency tradeoffs. We give concrete timings, and compare them with the current state of the art of group action primitives. Our implementation is publicly available at https://github.com/pegasis4d.

## 6.1 Parameter choices

Since our algorithm relies on higher dimensional representations of isogenies, we start with choosing primes of the form $p = c2^f - 1$, ensuring that the full $2^f$-torsion is defined over $\mathbb{F}_{p^2}$. Furthermore, since the probability that Equation (11) admits a solution, depends heavily on how close $2^f$ is to $p$, we restrict to the case where $c$ is small.

Among the possible choices for $c$ we retain those that maximize the number of small primes that split in $\mathcal{O}$, and in particular, primes that are 3 mod 4, since these will greatly increase the probability of finding $\mathfrak{B}$-good pairs $(u, v)$. Note that prime factors of $c$ will automatically be split.

To choose the set $\mathfrak{B}$, we proceed as follows: we only include split primes in $\mathfrak{B}$, which has 3 advantages: first, we can execute the whole algorithm over $\mathbb{F}_p$, instead of over $\mathbb{F}_{p^2}$, resulting in a major speed-up. Second, it allows us to use Elkies' algorithm to compute small degree isogenies without explicitly constructing kernel points over extensions of $\mathbb{F}_p$. Third: such split primes are the only ones that appear as factors of norms of primitive ideals in $\mathcal{O}$.

The number of primes to include in $\mathfrak{B}$ is mostly a tradeoff between steps 1 and 2: including more primes in $\mathfrak{B}$ makes Equation (11) easier to solve (step 1
becomes faster), but step 2 will have more (and larger) Elkies isogenies to compute. For each parameter set, we started by setting $\mathfrak{B} = \{2, 3\}$ and progressively added primes until computing $\ell$-isogenies would become too expensive. We then selected the set with the best overall timing. The final choices are reported in Table 1. A more detailed discussion of the selection process can be found in Appendix E.

<table><caption>**Table 1.** Parameter sets used in our implementation for different bitsizes. The prime $p$ is of the form $p = c2^f - 1$ and $\mathfrak{B}$ is the set of small split primes used.</caption><thead><tr><th>Parameter set</th><th>f</th><th>c</th><th>β</th></tr></thead><tbody><tr><td>500</td><td>503</td><td>33</td><td>2, 3, 7, 11, 13</td></tr><tr><td>1000</td><td>1004</td><td>15</td><td>2, 3, 5, 7, 11</td></tr><tr><td>1500</td><td>1551</td><td>9</td><td>2, 3, 5, 11</td></tr><tr><td>2000</td><td>2026</td><td>51</td><td>2, 3, 7, 11, 17</td></tr><tr><td>4000</td><td>4084</td><td>63</td><td>2, 3, 7, 11, 17, 19</td></tr></tbody></table>

## 6.2 Solving the norm equation

The first step (Section 5.1) begins with creating a reduced basis of the starting ideal $\mathfrak{a}$, followed by a search for small vectors. We are looking for short vectors of (reduced) norm roughly $\sqrt{p}$, but we allow them to be slightly larger at this point (up to $\log(p)\sqrt{p}$). We then attempt to solve Equation (11) by taking pairs of the shortest vectors. In general, the more pairs, the higher the probability to find a solution. However, the first combinations will be the smallest ones, hence the ones with the highest probability of success. In practice, we observed that taking combinations only involving the 10 shortest vectors, and rerandomizing upon failure, was the most efficient choice for all parameter sets.

Finally, for each pair of $(u, v)$ that we obtain, we have to attempt to write $u/g_u$ and $v/g_v$ as sums of squares. For this, we do trial division on both with all prime factors up to $10^4$ for the 500 and 1000 bits parameters, and up to $10^5$ for all larger parameters. If we find factors 3 mod 4 (with odd exponent) not included in $\mathfrak{B}$, we immediately discard the pair. Otherwise, we test what is left for primality. If they are prime, we have the factorization of $(u, v)$ and hence their decomposition as sum of two squares. Otherwise, we move on to the next pair.

## 6.3 Implementation results

The results of our SageMath 10.5 implementation can be found in Table 2; timings for each step are in seconds, and are obtained by averaging 100 runs on an Intel Core i5-1235U clocked at 4.0 GHz.
<table><caption>**Table 2.** Timings in seconds for SageMath 10.5 PEGASIS implementation to evaluate one group action at different prime sizes, measured on an Intel Core i5-1235U CPU with maximal clock speed of 4.0 GHz. The last column indicates the number of ideal-class rerandomizations required to find a solution in Step 1. These results are averaged over 100 runs. Step 1 is the time used to solve the norm equation, Step 2 is the time used to derive the kernel of the dimension 4 isogeny, and Step 3 is the time used to compute the dimension 4 isogeny.</caption><thead><tr><th rowspan="2">Prime size (bits)</th><th rowspan="2">Prime</th><th colspan="4">Time (s)</th><th rowspan="2">Rerand.</th></tr><tr><th>Step 1</th><th>Step 2</th><th>Step 3</th><th>Total</th></tr></thead><tbody><tr><td>508</td><td>3 · 11 · 2<sup>503</sup> − 1</td><td>0.097</td><td>0.48</td><td>0.96</td><td>1.53</td><td>0.17</td></tr><tr><td>1008</td><td>3 · 5 · 2<sup>1004</sup> − 1</td><td>0.21</td><td>1.16</td><td>2.84</td><td>4.21</td><td>0.07</td></tr><tr><td>1554</td><td>3<sup>2</sup> · 2<sup>1551</sup> − 1</td><td>1.19</td><td>2.85</td><td>6.49</td><td>10.5</td><td>1.53</td></tr><tr><td>2031</td><td>3 · 17 · 2<sup>2026</sup> − 1</td><td>1.68</td><td>8.34</td><td>11.3</td><td>21.3</td><td>0.70</td></tr><tr><td>4089</td><td>3<sup>2</sup> · 7 · 2<sup>4084</sup> − 1</td><td>15.6</td><td>52.8</td><td>53.5</td><td>122</td><td>0.41</td></tr></tbody></table>

We conclude this section with a more detailed comparison with the other available isogeny-based EGAs in the literature. The comparison is summarised in Table 3. The timings for SCALLOP [21] were measured on an Intel i5-6440HQ processor clocked at 3.5 GHz, while the timings for SCALLOP-HD [14] were measured on an Intel Alder Lake CPU core clocked at 2.1 GHz, and the timings for PEARL-SCALLOP [3] were measured on an Intel i5-1038NG7 processor clocked at 2.0 GHz. The timings for the two versions of KLaPoTi [50] were all re-measured on the same hardware setup as the timings presented in Table 2.

PEGASIS is the first instantiation of an EGA at the 2000-bit and 4000-bit security level. In fact, it is even the first time that the full CSIDH group action can be computed at the 1000-bit level. However, as Table 3 shows, PEGASIS also significantly outperforms all earlier works at their security levels. The closest comparison comes with the Rust implementation of KLaPoTi. However, the fact that KLaPoTi was able to achieve a speedup by two orders of magnitude simply by switching from a high-level SageMath implementation to a low-level Rust implementation is also very promising for PEGASIS. Under the assumption that a comparable speedup would be possible for PEGASIS, we see that PEGASIS gives a highly-practical EGA, even at the highest security levels.

One of the main reasons of this efficiency gain is that unlike all the other EGA instantiations, we are able to use the orientation given by Frobenius. This has three key benefits:

- Evaluating the orientation is very efficient,

- We do not need to represent, nor push forward the orientation,

- We can work over $\mathbb{F}_p$.

The price to pay is the need to go up to dimension 4 to compute the action, but as Table 2 shows, it is still reasonably efficient. We stress that our implementation is only a proof of concept, to explore whether level 4000 was feasible
<table><caption>**Table 3.** Comparison of PEGASIS and other effective group actions in the literature. Timings in seconds of different implementations to evaluate one group action at different (rounded) prime sizes. The SCALLOP variants are starred because they were measured on different hardware setups. The results of both KLaPoTi and PEGASIS are averaged over 100 runs and measured on the same hardware.</caption><thead><tr><th rowspan="3">Paper</th><th rowspan="3">Language</th><th colspan="6">Time (s)</th></tr><tr><th colspan="5">Approx. prime size (bits)</th><th rowspan="2">4000</th></tr><tr><th>500</th><th>1000</th><th>1500</th><th>2000</th></tr></thead><tbody><tr><td>SCALLOP [21]*</td><td>C++</td><td>35</td><td>750</td><td></td><td></td><td></td><td></td></tr><tr><td>SCALLOP-HD [14]*</td><td>Sage</td><td>88</td><td>1140</td><td></td><td></td><td></td><td></td></tr><tr><td>PEARL-SCALLOP [3]*</td><td>C++</td><td>30</td><td>58</td><td>710</td><td></td><td></td><td></td></tr><tr><td rowspan="2">KLaPoTi [50]</td><td>Sage</td><td>207</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Rust</td><td>1.95</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PEGASIS (This work)</td><td>Sage</td><td>1.53</td><td>4.21</td><td>10.5</td><td>21.3</td><td>121</td><td></td></tr></tbody></table>

in practice. In particular, we are missing many standard implementation tricks; as a simple example, we work with affine coordinates instead of projective coordinates. Despite this, our timings are very encouraging, and we hope a low level optimised implementation could even be made comparable in efficiency to a state-of-the-art implementation of CSIDH ran as a REGA, like in [11].

## 6.4 Further improvements

As described in Section 4.2, it is possible to force only one among $u/g_u$ and $v/g_v$ to be a sum of squares, and perform the other isogeny in dimension 4. Note that such 4-dimensional isogeny will have only half of the length of the isogeny that we perform at the end, thus being twice as fast (at least, since the isogeny computation is quasi-linear in its length). On the other hand, steps 1 and 2 would both become much faster. This approach hence becomes promising as soon as the 4-dimensional isogeny computation has a cost comparable to step 2, which is the case for the 2000 and 4000 parameter (see Table 2). For such parameters, we computed the first two steps of the algorithm, with $\mathfrak{B} = \{2, 3, 7\}$ in both cases; the results are reported in Table 4. We note that the number of rerandomizations stays comparable whilst the total time for step 1 and step 2 notably decreases.

Estimating the cost of one $u$-isogeny in dimension 4 as half of the cost of step 3 in Table 2, we expect to reach 21.3 seconds for the 2000 parameter and 106.2 seconds for the 4000. This approach is hence comparable to ours in the former case, and even faster in the latter. We leave its implementation as future work. On the other hand, doing both $u$ and $v$ using 4-dimensional isogenies is less efficient at all security levels we considered.
<table><caption>**Table 4.** Time taken to compute Step 1 and Step 2 when solving the norm equation with single sum of squares, measured in wall-clock seconds. The last column gives the number of rerandomizations needed. These results are averaged over 100 runs.</caption><thead><tr><th rowspan="2">Prime (bits)</th><th rowspan="2">Prime</th><th colspan="3">Time (s)</th><th rowspan="2">Rerand.</th></tr><tr><th>Step 1</th><th>Step 2</th><th>Total</th></tr></thead><tbody><tr><td>2031</td><td>3 · 17 · 2<sup>2026</sup> − 1</td><td>0.49</td><td>3.83</td><td>4.32</td><td>0.70</td></tr><tr><td>4089</td><td>3<sup>2</sup> · 7 · 2<sup>4084</sup> − 1</td><td>3.25</td><td>22.8</td><td>26.0</td><td>1.25</td></tr></tbody></table>

## References

[1] M. Abdalla, T. Eisenhofer, E. Kiltz, S. Kunzweiler, and D. Riepel. “Password-Authenticated Key Exchange from Group Actions”. In: *CRYPTO 2022, Part II*. Ed. by Y. Dodis and T. Shrimpton. Vol. 13508. LNCS. Springer, Cham, Aug. 2022, pp. 699–728. DOI: 10.1007/978-3-031-15979-4_24.

[2] N. Alamati, L. De Feo, H. Montgomery, and S. Patranabis. “Cryptographic Group Actions and Applications”. In: *ASIACRYPT 2020, Part II*. Ed. by S. Moriai and H. Wang. Vol. 12492. LNCS. Springer, Cham, Dec. 2020, pp. 411–439. DOI: 10.1007/978-3-030-64834-3_14.

[3] B. Allombert, J.-F. Biasse, J. K. Eriksen, P. Kutas, C. Leonardi, A. Page, R. Scheidler, and M. T. Bagi. *PEARL-SCALLOP: Parameter Extension Applicable in Real Life for SCALLOP*. 2024.

[4] A. Basso, G. Borin, W. Castryck, M. C.-R. Santos, R. Invernizzi, A. Leroux, L. Maino, F. Vercauteren, and B. Wesolowski. *PRISM: Simple And Compact Identification and Signatures From Large Prime Degree Isogenies*. Cryptology ePrint Archive, Report 2025/135. 2025. URL: https://eprint.iacr.org/2025/135.

[5] A. Basso, L. De Feo, P. Dartois, A. Leroux, L. Maino, G. Pope, D. Robert, and B. Wesolowski. *SQIsign2D-West: The Fast, the Small, and the Safer*. Cryptology ePrint Archive, Report 2024/760. 2024. URL: https://eprint.iacr.org/2024/760.

[6] D. J. Bernstein, T. Lange, C. Martindale, and L. Panny. “Quantum Circuits for the CSIDH: Optimizing Quantum Evaluation of Isogenies”. In: *EUROCRYPT 2019, Part II*. Ed. by Y. Ishai and V. Rijmen. Vol. 11477. LNCS. Springer, Cham, May 2019, pp. 409–441. DOI: 10.1007/978-3-030-17656-3_15.

[7] W. Beullens, T. Kleinjung, and F. Vercauteren. “CSI-FiSh: Efficient Isogeny Based Signatures Through Class Group Computations”. In: *ASIACRYPT 2019, Part I*. Ed. by S. D. Galbraith and S. Moriai. Vol. 11921. LNCS. Springer, Cham, Dec. 2019, pp. 227–247. DOI: 10.1007/978-3-030-34578-5_9.
[8] C. Birkenhake and H. Lange. *Complex Abelian Varieties*. Berlin, Heidelberg: Springer Berlin Heidelberg, 2004. ISBN: 978-3-662-06307-1. DOI: 10.1007/978-3-662-06307-1. URL: https://doi.org/10.1007/978-3-662-06307-1.

[9] X. Bonnetain and A. Schrottenloher. "Quantum Security Analysis of CSIDH". In: *EUROCRYPT 2020, Part II*. Ed. by A. Canteaut and Y. Ishai. Vol. 12106. LNCS. Springer, Cham, May 2020, pp. 493-522. DOI: 10.1007/978-3-030-45724-2_17.

[10] C. Boyd, G. T. Davies, K. Gjøsteen, and Y. Jiang. "Fast and Secure Updatable Encryption". In: *CRYPTO 2020, Part I*. Ed. by D. Micciancio and T. Ristenpart. Vol. 12170. LNCS. Springer, Cham, Aug. 2020, pp. 464-493. DOI: 10.1007/978-3-030-56784-2_16.

[11] F. Campos, A. Hellenbrand, M. Meyer, and K. Reijnders. *dCTIDH: Fast & Deterministic CTIDH*. Cryptology ePrint Archive, Paper 2025/107. 2025. URL: https://eprint.iacr.org/2025/107.

[12] W. Castryck and T. Decru. "CSIDH on the Surface". In: *Post-Quantum Cryptography - 11th International Conference, PQCrypto 2020*. Ed. by J. Ding and J.-P. Tillich. Springer, Cham, 2020, pp. 111-129. DOI: 10.1007/978-3-030-44223-1_7.

[13] W. Castryck, T. Lange, C. Martindale, L. Panny, and J. Renes. "CSIDH: An Efficient Post-Quantum Commutative Group Action". In: *ASIACRYPT 2018, Part III*. Ed. by T. Peyrin and S. Galbraith. Vol. 11274. LNCS. Springer, Cham, Dec. 2018, pp. 395-427. DOI: 10.1007/978-3-030-03332-3_15.

[14] M. Chen, A. Leroux, and L. Panny. "SCALLOP-HD: Group Action from 2-Dimensional Isogenies". In: *PKC 2024, Part II*. Ed. by Q. Tang and V. Teague. Vol. 14603. LNCS. Springer, Cham, Apr. 2024, pp. 190-216. DOI: 10.1007/978-3-031-57725-3_7.

[15] L. Colò and D. Kohel. *Orienting supersingular isogeny graphs*. Cryptology ePrint Archive, Report 2020/985. 2020. URL: https://eprint.iacr.org/2020/985.

[16] J.-M. Couveignes. *Hard Homogeneous Spaces*. Cryptology ePrint Archive, Report 2006/291. 2006. URL: https://eprint.iacr.org/2006/291.

[17] D. A. Cox. *Primes of the Form $x^2+ny^2$: Fermat, Class Field Theory, and Complex Multiplication*. 2nd edition. Wiley, 2013, p. 349.

[18] P. Dartois. *Fast computation of 2-isogenies in dimension 4 and cryptographic applications*. Cryptology ePrint Archive, Paper 2024/1180. 2024. URL: https://eprint.iacr.org/2024/1180.

[19] P. Dartois, A. Leroux, D. Robert, and B. Wesolowski. "SQIsignHD: New Dimensions in Cryptography". In: *Advances in Cryptology - EUROCRYPT 2024*. Ed. by M. Joye and G. Leander. Cham: Springer Nature Switzerland, 2024, pp. 3-32.

[20] P. Dartois, L. Maino, G. Pope, and D. Robert. *An Algorithmic Approach to (2,2)-isogenies in the Theta Model and Applications to Isogeny-based Cryptography*. Cryptology ePrint Archive, Paper 2023/1747. https://
eprint.iacr.org/2023/1747. 2023. URL: https://eprint.iacr.org/2023/1747.

[21] L. De Feo, T. B. Fouotsa, P. Kutas, A. Leroux, S.-P. Merz, L. Panny, and B. Wesolowski. "SCALLOP: Scaling the CSI-FiSh". In: *PKC 2023, Part I*. Ed. by A. Boldyreva and V. Kolesnikov. Vol. 13940. LNCS. Springer, Cham, May 2023, pp. 345-375. DOI: 10.1007/978-3-031-31368-4_13.

[22] L. De Feo and S. D. Galbraith. "SeaSign: Compact Isogeny Signatures from Class Group Actions". In: *EUROCRYPT 2019, Part III*. Ed. by Y. Ishai and V. Rijmen. Vol. 11478. LNCS. Springer, Cham, May 2019, pp. 759-789. DOI: 10.1007/978-3-030-17659-4_26.

[23] L. De Feo, D. Kohel, A. Leroux, C. Petit, and B. Wesolowski. "SQISign: Compact Post-quantum Signatures from Quaternions and Isogenies". In: *ASIACRYPT 2020, Part I*. Ed. by S. Moriai and H. Wang. Vol. 12491. LNCS. Springer, Cham, Dec. 2020, pp. 64-93. DOI: 10.1007/978-3-030-64837-4_3.

[24] L. De Feo and M. Meyer. "Threshold Schemes from Isogeny Assumptions". In: *PKC 2020, Part II*. Ed. by A. Kiayias, M. Kohlweiss, P. Wallden, and V. Zikas. Vol. 12111. LNCS. Springer, Cham, May 2020, pp. 187-212. DOI: 10.1007/978-3-030-45388-6_7.

[25] Y. Desmedt and Y. Frankel. "Threshold Cryptosystems". In: *CRYPTO'89*. Ed. by G. Brassard. Vol. 435. LNCS. Springer, New York, Aug. 1990, pp. 307-315. DOI: 10.1007/0-387-34805-0_28.

[26] W. Diffie and M. E. Hellman. "New Directions in Cryptography". In: *IEEE Transactions on Information Theory* 22.6 (1976), pp. 644-654. DOI: 10.1109/TIT.1976.1055638.

[27] M. Duparc and T. B. Fouotsa. *SQIPrime: A dimension 2 variant of SQISignHD with non-smooth challenge isogenies*. Cryptology ePrint Archive, Report 2024/773. 2024. URL: https://eprint.iacr.org/2024/773.

[28] B. Edixhoven, G. van der Geer, and B. Moonen. *Abelian Varieties*. URL: http://van-der-geer.nl/~gerard/AV.pdf.

[29] T. ElGamal. "A Public Key Cryptosystem and a Signature Scheme Based on Discrete Logarithms". In: *CRYPTO'84*. Ed. by G. R. Blakley and D. Chaum. Vol. 196. LNCS. Springer, Berlin, Heidelberg, Aug. 1984, pp. 10-18. DOI: 10.1007/3-540-39568-7_2.

[30] N. D. Elkies. "Elliptic and modular curves over finite fields and related computational issues". In: *Computational Perspectives on Number Theory: Proceedings of a Conference in Honor of A. O. L. Atkin*. Ed. by D. A. Buell and J. T. Teitelbaum. Vol. 7. AMS/IP Studies in Advanced Mathematics. American Mathematical Society, International Press, 1998, pp. 21-76.

[31] L. D. Feo. *Mathematics of Isogeny Based Cryptography*. 2017. arXiv: 1711.04062 [cs.CR]. URL: https://arxiv.org/abs/1711.04062.

[32] S. Galbraith, L. Panny, B. Smith, and F. Vercauteren. *Quantum Equivalence of the DLP and CDHP for Group Actions*. Cryptology ePrint
Archive, Report 2018/1199. 2018. URL: https://eprint.iacr.org/2018/1199.

[33] S. D. Galbraith, Y.-F. Lai, and H. Montgomery. "A Simpler and More Efficient Reduction of DLog to CDH for Abelian Group Actions". In: *PKC 2024, Part II*. Ed. by Q. Tang and V. Teague. Vol. 14603. LNCS. Springer, Cham, Apr. 2024, pp. 36–60. DOI: 10.1007/978-3-031-57725-3_2.

[34] R. Gennaro, S. Jarecki, H. Krawczyk, and T. Rabin. "Robust Threshold DSS Signatures". In: *EUROCRYPT'96*. Ed. by U. M. Maurer. Vol. 1070. LNCS. Springer, Berlin, Heidelberg, May 1996, pp. 354–371. DOI: 10.1007/3-540-68339-9_31.

[35] L. Harn. "Group-oriented (t, n) threshold digital signature scheme and digital multisignature". In: *IEE Proceedings-Computers and Digital Techniques* 141.5 (1994), pp. 307–313.

[36] D. P. Jablon. "Strong password-only authenticated key exchange". In: *ACM SIGCOMM Computer Communication Review* 26.5 (1996), pp. 5–26.

[37] N. Jacobson. *Basic Algebra I*. Second edition. 1984.

[38] D. Kohel, K. Lauter, C. Petit, and J.-P. Tignol. "On the quaternion-isogeny path problem". In: *LMS Journal of Computation and Mathematics* 17.A (2014), pp. 418–432. DOI: 10.1112/S1461157014000151.

[39] G. Kuperberg. "A subexponential-time quantum algorithm for the dihedral hidden subgroup problem". In: *SIAM Journal on Computing* 35.1 (2005), pp. 170–188.

[40] A. Leroux and M. Roméas. "Updatable Encryption from Group Actions". In: *Post-Quantum Cryptography - 15th International Workshop, PQCrypto 2024, Part II*. Ed. by M.-J. Saarinen and D. Smith-Tone. Springer, Cham, June 2024, pp. 20–53. DOI: 10.1007/978-3-031-62746-0_2.

[41] J. S. Milne. "Abelian Varieties". In: *Arithmetic Geometry*. New York, NY: Springer New York, 1986, pp. 103–150. ISBN: 978-1-4613-8655-1. DOI: 10.1007/978-1-4613-8655-1_5.

[42] H. Montgomery and M. Zhandry. "Full Quantum Equivalence of Group Action DLog and CDH, and More". In: *ASIACRYPT 2022, Part I*. Ed. by S. Agrawal and D. Lin. Vol. 13791. LNCS. Springer, Cham, Dec. 2022, pp. 3–32. DOI: 10.1007/978-3-031-22963-3_1.

[43] T. Moriya, H. Onuki, and T. Takagi. "SiGamal: A Supersingular Isogeny-Based PKE and Its Application to a PRF". In: *ASIACRYPT 2020, Part II*. Ed. by S. Moriai and H. Wang. Vol. 12492. LNCS. Springer, Cham, Dec. 2020, pp. 551–580. DOI: 10.1007/978-3-030-64834-3_19.

[44] D. Mumford. "On the equations defining abelian varieties 1". In: *Inventiones mathematicae* 1.4 (1966), pp. 287–354. DOI: 10.1007/BF01389737.

[45] D. Mumford. *Abelian varieties*. Second Edition. Tata Institute of fundamental research studies in mathematics. London: Oxford University Press, 1974, pp. x+279.
[46] K. Nakagawa and H. Onuki. "QFESTA: Efficient Algorithms and Parameters for FESTA Using Quaternion Algebras". In: CRYPTO 2024, Part V. Ed. by L. Reyzin and D. Stebila. Vol. 14924. LNCS. Springer, Cham, Aug. 2024, pp. 75–106. DOI: 10.1007/978-3-031-68388-6_4.

[47] K. Nakagawa and H. Onuki. SQIsign2D-East: A New Signature Scheme Using 2-dimensional Isogenies. Cryptology ePrint Archive, Report 2024/771. 2024. URL: https://eprint.iacr.org/2024/771.

[48] H. Onuki. "On oriented supersingular elliptic curves". In: Finite Fields and Their Applications 69 (2021), p. 101777. ISSN: 1071-5797. DOI: https://doi.org/10.1016/j.ffa.2020.101777. URL: https://www.sciencedirect.com/science/article/pii/S1071579720301465.

[49] A. Page and D. Robert. Introducing Clapoti(s): Evaluating the isogeny class group action in polynomial time. Cryptology ePrint Archive, Report 2023/1766. 2023. URL: https://eprint.iacr.org/2023/1766.

[50] L. Panny, C. Petit, and M. Stopar. "KLaPoTi: An asymptotically efficient isogeny group action from 2-dimensional isogenies". In: IACR Cryptol. ePrint Arch. (2024), p. 1844. URL: https://eprint.iacr.org/2024/1844.

[51] C. Peikert. "He Gives C-Sieves on the CSIDH". In: EUROCRYPT 2020, Part II. Ed. by A. Canteaut and Y. Ishai. Vol. 12106. LNCS. Springer, Cham, May 2020, pp. 463–492. DOI: 10.1007/978-3-030-45724-2_16.

[52] J. Renes. "Computing Isogenies Between Montgomery Curves Using the Action of (0, 0)". In: Post-Quantum Cryptography - 9th International Conference, PQCrypto 2018. Ed. by T. Lange and R. Steinwandt. Springer, Cham, 2018, pp. 229–247. DOI: 10.1007/978-3-319-79063-3_11.

[53] D. Robert. "Theta functions and cryptographic applications". PhD thesis. Université Henri-Poincaré, Nancy 1, France, July 2010. URL: http://www.normalesup.org/~robert/pro/publications/academic/phd.pdf.

[54] D. Robert. Some applications of higher dimensional isogenies to elliptic curves (overview of results). Cryptology ePrint Archive, Report 2022/1704. 2022. URL: https://eprint.iacr.org/2022/1704.

[55] D. Robert. A note on optimising 2ⁿ-isogenies in higher dimension. http://www.normalesup.org/~robert/pro/publications/notes/2023-06-optimising_isogenies.pdf. 2023.

[56] D. Robert. "Breaking SIDH in Polynomial Time". In: EUROCRYPT 2023, Part V. Ed. by C. Hazay and M. Stam. Vol. 14008. LNCS. Springer, Cham, Apr. 2023, pp. 472–503. DOI: 10.1007/978-3-031-30589-4_17.

[57] D. Robert. The geometric interpretation of the Tate pairing and its applications. Cryptology ePrint Archive, Report 2023/177. 2023. URL: https://eprint.iacr.org/2023/177.

[58] D. Robert. On the efficient representation of isogenies (a survey). Cryptology ePrint Archive, Report 2024/1071. 2024. URL: https://eprint.iacr.org/2024/1071.

[59] D. Robert and N. Sarkis. "Computing 2-isogenies between Kummer lines". In: CiC 1.1 (2024), p. 26. DOI: 10.62056/abvua69p1.
[60] A. Rostovtsev and A. Stolbunov. *Public-Key Cryptosystem Based On Iso-
genies*. Cryptology ePrint Archive, Report 2006/145. 2006. URL: https://eprint.iacr.org/2006/145.

[61] M. C.-R. Santos and K. Reijnders. *Return of the Kummer: a Toolbox for Genus-2 Cryptography*. Cryptology ePrint Archive, Paper 2024/948. 2024. URL: https://eprint.iacr.org/2024/948.

[62] C.-P. Schnorr. "Efficient Identification and Signatures for Smart Cards". In: *CRYPTO '89*. Ed. by G. Brassard. Vol. 435. LNCS. Springer, New York, Aug. 1990, pp. 239–252. DOI: 10.1007/0-387-34805-0_22.

[63] P. W. Shor. "Algorithms for Quantum Computation: Discrete Logarithms and Factoring". In: *35th FOCS*. IEEE Computer Society Press, Nov. 1994, pp. 124–134. DOI: 10.1109/SFCS.1994.365700.

[64] J. Vélu. "Isogénies entre courbes elliptiques". In: *Comptes-rendus de l'Académie des Sciences* 273 (1971). Available at https://gallica.bnf.fr, pp. 238–241.

[65] W. C. Waterhouse. "Abelian varieties over finite fields". eng. In: *Annales scientifiques de l'École Normale Supérieure* 2.4 (1969), pp. 521–560. URL: http://eudml.org/doc/81852.

[66] G. Zanon, M. A. Simplício Jr., G. C. C. F. Pereira, J. Doliskani, and P. S. L. M. Barreto. "Faster Isogeny-Based Compressed Key Agreement". In: *Post-Quantum Cryptography - 9th International Conference, PQCrypto 2018*. Ed. by T. Lange and R. Steinwandt. Springer, Cham, 2018, pp. 248–268. DOI: 10.1007/978-3-319-79063-3_12.
# A An introduction to efficient higher-dimensional isogeny computation

In this section, we explain how higher-dimensional isogenies are computed with the theta model. We use the algorithms introduced in [18] for SQIsignHD [19] and the attacks against SIDH. In Appendix A.1, we first introduce some basic theory on theta coordinates in a cryptographer-friendly way. We then briefly explain how to change theta coordinates to obtain a system of theta coordinates adapted to a 2$^e$-isogeny, how such an isogeny is computed and how to obtain theta coordinates from Montgomery (x:z)-coordinates on a product of Montgomery elliptic curves. In Appendix A.2, we give more details on changes of theta coordinates adapted to an isogeny obtained from Kani's lemma and how to decompose the codomain into a product.

## A.1 Some theory on theta functions

We recall some basics on the theta model. We refer to [20, § 2.1] for a cryptographer friendly introduction. Let $(A, \lambda)$ be a principally polarised abelian variety of dimension $g$ over an algebraically closed field $k$. The polarisation $\lambda : A \to \hat{A}$ is induced by an ample divisor $D$ on $A$ as follows $\lambda = \lambda_D : A \to \hat{A}, x \mapsto [t_x^*D - D]$. Without loss of generality, we can assume that $D$ is symmetric $[-1]^*D \sim D$. Since $\lambda$ is a principal polarisation, $\ker(\lambda) = \{0\}$. We define the polarisation of level $n$ as $\lambda \circ [n] = \lambda_{nD} : x \mapsto [t_x^{*n}D - nD]$ whose kernel is $\ker(\lambda \circ [n]) = A[n]$.

**Theta-structures.** Formally, a *theta-structure* $\Theta_n$ of level $n$ on $(A, \lambda \circ [n])$ is an isomorphism between the Heisenberg group $\mathcal{H}(n)$ of level $n$ and the theta-group $\mathcal{G}(nD)$ of $nD$, $\Theta_n : \mathcal{G}(nD) \to \mathcal{H}(n)$, that we do not define here. However, this level of abstraction is not needed in the following. The reader only needs to know that a theta structure $\Theta_n$ induces a system of coordinates $(\theta_i)_{i \in (Z/nZ)^g}$ (where $g = \dim(A))$ called *theta-coordinates*. They define a map into the projective space $A \to \mathbb{P}_{k}^{n^g-1}$, $x \mapsto (\theta_i(x))_{i \in (Z/nZ)^g}$ which is an embedding when $n \ge 3$ [45, p. 163] and induces an embedding of the Kummer variety $A/\pm \hookrightarrow \mathbb{P}_{k}^{2^g-1}$ when $n=2$ and $A$ is not a product [8, Theorem 4.8.1]. In practice, we work in level $n=2$ so theta-coordinates represent points on the Kummer variety *i.e.* points up to a sign.

A polarised abelian variety with a level $n$ theta structure $(A, \lambda \circ [n], \Theta_n)$ is generally too much data to represent on a computer. Usually, we represent it by its associated *theta null point* $(\theta_i(0))_{i \in (Z/nZ)^g}$. When $n=2$ we do not have a guarantee that the theta null point fully determines $(A, \lambda \circ [n], \Theta_n)$. However, this information is enough in practice to perform arithmetic computations on a polarised abelian variety $(A, \lambda \circ [2])$ that is neither a product nor 2-isogenous to a product. These arithmetic operations include doubling $x \mapsto 2x$ and differential addition $x, y, x-y \mapsto x+y$ (see [53, Algorithm 4.4.10]).
**Change of theta-coordinates.** A symmetric theta structure of level $n$ is fully determined by a symplectic basis of $A[2n]$ [44, Remark 3, p. 319]. Such a basis $(S_1, \dots, S_g, T_1, \dots, T_g)$ of $A[2n]$ satisfies:

$$ \forall 1 \le j \le n, \quad e_{2n}(S_i, S_j) = e_{2n}(T_i, T_j) = 1 \quad \text{and} \quad e_{2n}(S_i, T_j) = \zeta^{\delta_{i,j}}, $$

where $e_{2n}$ is the $2n$-th Weil pairing and $\zeta \in k$ is a primitive $2n$-th root of unity. Such a basis is called a $\zeta$-symplectic basis. In the following, we shall always assume that theta structures are symmetric.

A change of symplectic basis induces a change of theta structure, hence a change of theta coordinates. Namely, let $\Theta_n$ and $\Theta'_n$ be two theta-structures and respectively induced by $\zeta$-symplectic basis $\mathcal{B}$ and $\mathcal{B}'$ of $A[2n]$. Let $M \in \mathrm{Sp}_{2g}(\mathbb{Z}/n\mathbb{Z})$ be the change of basis matrix from $\mathcal{B}$ to $\mathcal{B}'$. Let $(\theta_i)_i$ and $(\theta'_i)_i$ be the systems of theta-coordinates associated to $\Theta_n$ and $\Theta'_n$ respectively. Then, there is an explicit matrix $N(M, \zeta) \in \mathrm{GL}_n(k)$ such that $(\theta'_i)_i = N(M, \zeta) \cdot (\theta_i)_i$ that can be computed as follows.

**Theorem A.1.** [18, Theorem 12] Let us write the symplectic change of basis matrix from $\mathcal{B}$ to $\mathcal{B}'$ by $g \times g$ blocks:

$$ M := \begin{pmatrix} A & C \\ B & D \end{pmatrix} $$

Then, there exists $i_0 \in (\mathbb{Z}/n\mathbb{Z})^g$ and $\lambda \in k^*$ such that for all $i \in (\mathbb{Z}/n\mathbb{Z})^g$,

$$ \theta'_i = \lambda \sum_{j \in (\mathbb{Z}/n\mathbb{Z})^g} \zeta^{\langle i|j \rangle - \langle Ai+Cj+2i_0|Bi+Dj \rangle} \theta_{Ai+Cj+i_0}, $$

where $\langle i|j \rangle := \sum_{k=1}^g i_k j_k$ is the usual scalar product. We can choose any value of $i_0 \in (\mathbb{Z}/n\mathbb{Z})^g$ such that

$$ \sum_{j \in (\mathbb{Z}/n\mathbb{Z})^g} \zeta^{-\langle Cj+2i_0|Dj \rangle} \theta_{i_0+Cj} \neq 0. $$

**Isogeny computations.** Let $F : A \to B$ be a $d$-isogeny. If $\mathcal{B} = (S_1, \dots, S_g, T_1, \dots, T_g)$ is a $\zeta$-symplectic basis of $A[4d]$ adapted to $F$ i.e. such that $\ker(F) = [4]\langle T_1, \dots, T_g \rangle$, then $\mathcal{C} = ([d]F(S_1), \dots, [d]F(S_g), F(T_1), \dots, F(T_g))$ is a $\zeta^d$-symplectic basis of $B[4]$ [19, Theorem 56]. In particular, it means that if we work with theta coordinates $(\theta_i^A)_{i \in (\mathbb{Z}/2\mathbb{Z})^g}$ induced by the symmetric theta structure associated to $[d]\mathcal{B}$ on the domain, then the codomain theta coordinates $(\theta_i^B)_{i \in (\mathbb{Z}/2\mathbb{Z})^g}$ will be induced by the symmetric theta structure associated to $\mathcal{C}$.

When $d=2$, [18, § 4.2] provides formulas to compute the codomain theta null point $(\theta_i^B(0))_{i \in (\mathbb{Z}/2\mathbb{Z})^g}$ from the theta coordinates of $T_1, \dots, T_g$. The codomain theta null point $(\theta_i^B(0))_i$ is then sufficient to evaluate $F$ i.e. to compute $(\theta_i^B(F(x)))_i$ given $(\theta_i^B(x))_i$ for all $x \in A(k)$. Hence the codomain theta null point not only represents the codomain $B$ but also the $2$-isogeny $F$ itself. When
we say that we compute $F$, we usually say that we compute its codomain theta null point (and relevant information trivially derived from it). Note the main difference with Vélu's formulas: 8-torsion points $T_1, \dots, T_g$ are necessary for this computation; 4 or 2-torsion points would not be sufficient.

When $d = 2^e$, and $F = f_e \circ \dots \circ f_1$ is a chain of 2-isogenies, a choice of symplectic basis $\mathcal{B}$ adapted to $F$ (i.e. such that $\ker(F) = [4]\langle T_1, \dots, T_g \rangle$) induces a choice of level 2 theta-structure on the codomain of $f_i$ associated to $[2^{e-i}]f_i \circ \dots \circ f_1([2^i]S_1, \dots, [2^i]S_g, T_1, \dots, T_g)$ for all $i \in [1; e]$, so that we can propagate the formulas from [18, § 4.2] along the chain.

When $F$ has been obtained via Kani's lemma, it is defined on a product of abelian varieties (e.g. a product of elliptic curves) so the first steps of the chain are *gluing isogenies* $f: A_1 \times A_2 \to B$ that need to be handled with extra care since the formulas are different to compute them. *Splitting isogenies* along the chain $f: A \to B_1 \times B_2$ also require some attention because point duplications $x \mapsto 2x$ are more difficult on their domain and codomain than generic ones. In particular, we need to detect these gluings and splittings in advance. This will be the main focus of Appendix B in our special case of interest (when $F$ is the isogeny from Section 5.3).

As in dimension 1, divide and conquer strategies can be applied to minimize the number of duplications and isogeny evaluations required to compute a 2-isogeny chain $F = f_e \circ \dots \circ f_1$. These strategies can also be adapted to take into account the specific constraints imposed by gluing and splitting isogenies and their relative cost. With 4-dimensional computations in mind, more details on these strategies can be found in [18, Appendix E].

**Product theta-structures.** When we compute isogenies obtained via Kani's lemma, we start and terminate on a product of abelian varieties which are naturally equipped with a *product theta-structure*. If $A = A_1 \times A_2$ and $(A_1, \lambda_1), (A_2, \lambda_2)$ are principally polarised abelian varieties, we consider the product polarisation $\lambda := \lambda_1 \times \lambda_2 : (x, y) \mapsto (\lambda_1(x), \lambda_2(x))$. If $\Theta_n^{A_i}$ is a level $n$ theta-structure on $(A_i, \lambda_i \circ [n])$ induced by a $\zeta$-symplectic basis $\mathcal{B}_i = (S_1^{(i)}, \dots, S_{g_i}^{(i)}, T_1^{(i)}, \dots, T_{g_i}^{(i)})$ of $A[2n]$ and of associated system of theta-coordinates $(\theta_j^{A_i})_{j \in (\mathbb{Z}/n\mathbb{Z})^{g_i}}$ for $i \in \{1, 2\}$, then we can define the *product theta-structure* $\Theta_n := \Theta_n^{A_1} \times \Theta_n^{A_2}$ on $(A, \lambda)$ as follows. It is the level $n$ theta-structure induced by the $\zeta$-symplectic basis of $A[2n]$:

$$ \mathcal{B} = \mathcal{B}_1 \times \mathcal{B}_2 := ((S_1^{(1)}, 0), \dots, (S_{g_1}^{(1)}, 0), (0, S_1^{(2)}), \dots, (0, S_{g_2}^{(2)}), \\ (T_1^{(1)}, 0), \dots, (T_{g_1}^{(1)}, 0), (0, T_1^{(2)}), \dots, (0, T_{g_2}^{(2)})). \quad (16) $$

The theta-coordinates $(\theta_j^A)_{j \in (\mathbb{Z}/n\mathbb{Z})^{g_1+g_2}}$ associated to $\Theta_n = \Theta_n^{A_1} \times \Theta_n^{A_2}$ are given as follows:

$$ \forall (x, y) \in A_1(k) \times A_2(k), (j_1, j_2) \in (\mathbb{Z}/n\mathbb{Z})^{g_1+g_2}, \quad \theta_{j_1, j_2}^A(x, y) = \theta_{j_1}^{A_1}(x)\theta_{j_2}^{A_2}(y). \quad (17) $$
In the context of a 2<sup>e</sup>-isogeny computation F : A → B (n = 2<sup>e+2</sup>), we usually start with a product theta-structure induced by 𝒰 and we have to compute a change of theta-coordinates induced by a change of symplectic basis from 𝒰 to 𝒰' = (S<sub>1</sub>, ..., S<sub>g</sub>, T<sub>1</sub>, ..., T<sub>g</sub>) such that ker(F) = [4]⟨T<sub>1</sub>, ..., T<sub>g</sub>⟩.

The codomain is also often a product B = B<sub>1</sub> × B<sub>2</sub> but the level 2 theta-structure on the codomain induced by F is associated to C = ([2<sup>e</sup>]F(S<sub>1</sub>), ..., [2<sup>e</sup>]F(S<sub>g</sub>), F(T<sub>1</sub>), ..., F(T<sub>g</sub>)) which is not a product theta-structure. A change of theta-coordinates is necessary to recover the product theta-coordinates on B. Once this has been done, we obtain the *theta-null point* of the product (θ<sup>B</sup><sub>j<sub>1</sub>,j<sub>2</sub></sub>(0))<sub>j<sub>1</sub>,j<sub>2</sub></sub> = (θ<sup>B<sub>1</sub></sup><sub>j<sub>1</sub></sub>(0) · θ<sup>B<sub>2</sub></sup><sub>j<sub>2</sub></sub>(0))<sub>j<sub>1</sub>,j<sub>2</sub></sub> from which we can infer the theta-null points (θ<sup>B<sub>i</sub></sup><sub>j</sub>(0))<sub>j</sub> of each component B<sub>i</sub>, which can then be used to determine the polarised abelian variety (B<sub>i</sub>, λ<sub>B<sub>i</sub></sub> ◦ [2]). In Appendix A.2, we shall see how to compute the necessary change of theta-coordinates for this splitting operation when F is obtained from Kani's lemma and apply this method to our case of interest (when F is the isogeny from Section 5.3) in Appendix B.4.

**Theta-structures on Montgomery curves.** In our context, we work with isogenies between products of elliptic curves. These elliptic curves are equipped with a Montgomery model over $\mathbb{F}_{p^2}$ (and even $\mathbb{F}_p$ for CSIDH) so we need to convert this model into a level 2 theta model. Conversely, we also need to translate a level 2 theta model obtained after the splitting operation mentioned above into a Montgomery model. The focus of this paragraph will be to explain these conversions.

If E is a Montgomery curve, then Montgomery (x : z)-coordinates and level 2 theta coordinates (θ₀ : θ₁) both induce *Kummer lines* i.e. isomorphisms E/± ↦ P¹. Hence, a change of coordinates between theta and Montgomery is a homography. In [55, Chapter 7, Appendix A], formulas were introduced to express such a homography when the theta structure is symmetric and associated to a 4-torsion basis of specific shape.

Let E be an elliptic curve in the Montgomery model and (P, Q) be a basis of E[4] such that Q := (-1 : 1) (in (x : z) coordinates). Let us write P := (r : s). Then, we may define a level 2 theta-structure on E with theta-null point (a : b) := (r + s : r - s). The conversion map from Montgomery to theta coordinates is then (x : z) ↦ (a(x - z) : b(x + z)).

Conversely, if (a : b) is the theta-null point, the conversion map from theta to
Montgomery coordinates is (θ₀, θ₁) ↦ (aθ₁ + bθ₀ : aθ₁ - bθ₀). Via this map, the
2-torsion theta points (b : a), (a : -b) and (b : -a) are mapped to the 2-torsion
Montgomery points (a² + b² : a² - b²), (0 : 1) and (a² - b² : a² + b²) respectively.
Identifying the 2-torsion, we can then obtain the Montgomery equation of E:

$$
By^2 = x(x - \alpha)(x - 1/\alpha) = x^3 + Ax^2 + x,
$$

with α := (a² + b²)/(a² - b²) and A := -(α + 1/α) = -2(a⁴ + b⁴)/(a⁴ - b⁴). To identify B, we need to know a full point of 4-torsion (x : y : z) on E. This will be important in our context to distinguish the result of the CSIDH class group action E<sub>a</sub> from its quadratic twist over F<sub>p</sub> (see Appendix B.4).
As we have seen previously, a fixed Montgomery curve $E$ has several level 2 theta structures related by symplectic matrices between 4-torsion basis. Using Theorem A.1, we can then obtain theta structures which are not induced by a 4-torsion basis $(P, Q)$ with $(x(Q) : z(Q)) \neq (-1 : 1)$. Conversely, from different theta structures on $E$, we obtain 6 different possible Montgomery coefficients $A$ from the conversion formula $(\theta_0, \theta_1) \mapsto (a\theta_1 + b\theta_0 : a\theta_1 - b\theta_0)$.

**Proposition A.1.** Let $(a : b)$ be a level 2 theta null point on a Montgomery curve $E$ associated to a symmetric theta structure. Then all possible Montgomery coefficients of $E$ are of the form:

$$ -2\frac{a^4 + b^4}{a^4 - b^4}, \quad -\frac{a^4 + 6a^2b^2 + b^4}{2(a^3b + b^3a)}, \quad \zeta_4\frac{a^4 - 6a^2b^2 + b^4}{2(a^3b - b^3a)}, \\ 2\frac{a^4 + b^4}{a^4 - b^4}, \quad \frac{a^4 + 6a^2b^2 + b^4}{2(a^3b + b^3a)}, \quad -\zeta_4\frac{a^4 - 6a^2b^2 + b^4}{2(a^3b - b^3a)}. $$

where $\zeta_4$ is a root of $-1$. They correspond respectively to theta null points:

$$ (a : b), \quad (a + b : a - b), \quad (a + \zeta_4 b : a - \zeta_4 b), \\ (b : a), \quad (a - b : a + b), \quad (a - \zeta_4 b : a + \zeta_4 b). $$

*Proof.* Let $(P, Q)$ be a basis of $E[4]$ inducing a symmetric theta structure of theta null point $(a : b)$ and let $\zeta_4 := e_4(P, Q)$. Then, enumerating all symplectic matrices of $Sp_2(Z/4Z)$, we obtain all $\zeta_4$-symplectic basis $(P', Q')$ of $E[4]$ inducing symmetric level 2 theta structures of theta null points $(a' : b')$ obtained from $(a : b)$ by Theorem A.1. We can then compute all corresponding Montgomery coefficients $A' := -2(a'^4 + b'^4)/(a'^4 - b'^4)$. $\square$

## A.2 Theta structures on the domain and codomain of an isogeny obtained from Kani's lemma

In the following lemma, we give an explicit description of a symplectic basis adapted to an isogeny obtained from Kani's lemma (in point (i)). By Theorem A.1, this can be used to compute a change of theta coordinates from product theta coordinates on the domain before starting the isogeny computation. We also describe a symplectic basis associated to a product level 2 theta structure on the codomain (in point (ii)) and its relation with the symplectic basis naturally induced by $F$ via [19, Theorem 56] (in point (iii)). Using Theorem A.1, the resulting symplectic change of basis matrix can be used to compute the product theta null point and product theta coordinates on the codomain. In particular, we can "extract" components of the codomain.

**Lemma A.1.** Let $a$ and $b$ be odd and coprime integers and $d := a+b$ not divisible by $\text{char}(k)$. Consider an $(a, b)$-isogeny diamond (as defined in Lemma 2.3)
between PPAVs of dimension g:

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

and the associated Kani d-isogeny:

$$F := \begin{pmatrix} \varphi & \tilde{\psi}' \\ -\psi & \tilde{\varphi}' \end{pmatrix} : A \times B' \to B \times A'$$

Let $(x_1, \dots, x_g, y_1, \dots, y_g)$ be a $\zeta$-symplectic basis of $B[4d]$. Let $\alpha$ and $\beta$ be modular inverses of a and b modulo 4d respectively.

(i) For all $i \in [1; g]$, we denote:

$$S_i := ([-\alpha]\tilde{\varphi}(y_i), 0), \quad S_{i+g} := (0, [\beta]\psi'(x_i)),$$

$$T_i := (\tilde{\varphi}(x_i), \psi'(x_i)), \quad T_{i+g} := ([1 - \alpha d]\tilde{\varphi}(y_i), \psi'(y_i)).$$

Then $(S_1, \dots, S_{2g}, T_1, \dots, T_{2g})$ is a $\zeta$-symplectic basis of $(A \times B')[4d]$ adapted to $F$ i.e. such that $\ker(F) = [4]\langle T_1, \dots, T_{2g} \rangle$.

(ii) For all $i \in [1; g]$, we denote:

$$U_i := ([d]x_i, 0), \quad U_{i+g} := (0, [d]\psi \circ \tilde{\varphi}(x_i)),$$

$$V_i := ([d]y_i, 0), \quad V_{i+g} := (0, [d\alpha\beta]\psi \circ \tilde{\varphi}(y_i)).$$

Then $(U_1, \dots, U_{2g}, V_1, \dots, V_{2g})$ is a product $\zeta^d$-symplectic basis of $(B \times A')[4]$.

(iii) The $\zeta^d$-symplectic basis $([d]F(S_1), \dots, [d]F(S_{2g}), F(T_1), \dots, F(T_{2g}))$ of $(B \times A')[4]$ naturally induced by $F$ via [19, Theorem 56] is related to $(U_1, \dots, U_{2g}, V_1, \dots, V_{2g})$ by the following formulas. For all $i \in [1; g]$:

$$[d]F(S_i) = -V_i + [b]V_{i+g}, \quad [d]F(S_{i+g}) = U_i + [\beta]U_{i+g},$$

$$F(T_i) = U_i, \quad F(T_{i+g}) = [b]V_{i+g},$$

so that:

$$U_i = F(T_i), \quad U_{i+g} = [bd]F(S_{i+g}) - [b]F(T_i),$$

$$V_i = F(T_{i+g}) - [d]F(S_i), \quad V_{i+g} = [\beta]F(T_{i+g}).$$

*Proof. (i)* For all principally polarised abelian varieties $C$ and $D$, we denote $\lambda_C$ the principal polarisation on $C$ and $\lambda_C \times \lambda_D$ the product (principal) polarisation on $C \times D$. Then, since $(x_1, \dots, x_g, y_1, \dots, y_g)$ is a $\zeta$-symplectic basis of $A[4d]$, we have the following results for the 4d-th Weil-pairing associated to the product polarisation $\lambda_A \times \lambda_{B'}$:
$$e_{4d}^{\lambda_A \times \lambda_{B'}}(S_i, S_j) = e_{4d}^{\lambda_A}([-\alpha]\tilde{\varphi}(y_i), [-\alpha]\tilde{\varphi}(y_j)) = e_{4d}^{\lambda_A}(y_i, y_j)^{a\alpha} = 1$$

$$e_{4d}^{\lambda_A \times \lambda_{B'}}(S_i, S_{j+g}) = e_{4d}^{\lambda_A}([-\alpha]\tilde{\varphi}(y_i), 0)e_{4d}^{\lambda_{B'}}(0, [\beta]\psi'(x_j)) = 1$$

$$e_{4d}^{\lambda_A \times \lambda_{B'}}(T_i, T_j) = e_{4d}^{\lambda_A}(\tilde{\varphi}(x_i), \tilde{\varphi}(x_j))e_{4d}^{\lambda_{B'}}(\psi'(x_i), \psi'(x_j)) = e_{4d}^{\lambda_A}(x_i, x_j)^{a+b} = 1$$

$$e_{4d}^{\lambda_A \times \lambda_{B'}}(T_i, T_{j+g}) = e_{4d}^{\lambda_A}(\tilde{\varphi}(x_i), [1 - \alpha d]\tilde{\varphi}(y_j))e_{4d}^{\lambda_{B'}}(\psi'(x_i), \psi'(y_j)) \\ = e_{4d}^{\lambda_A}(x_i, y_j)^{a(1-\alpha d)+b} = e_{4d}^{\lambda_A}(x_i, y_j)^0 = 1$$

$$e_{4d}^{\lambda_A \times \lambda_{B'}}(S_i, T_j) = e_{4d}^{\lambda_A}([-\alpha]\tilde{\varphi}(y_i), \tilde{\varphi}(x_j))e_{4d}^{\lambda_{B'}}(0, \psi'(x_j)) = e_{4d}^{\lambda_A}(x_j, y_i)^{a\alpha} = \zeta^{\delta_{i,j}}$$

$$e_{4d}^{\lambda_A \times \lambda_{B'}}(S_i, T_{j+g}) = e_{4d}^{\lambda_A}([-\alpha]\tilde{\varphi}(y_i), [1 - \alpha d]\tilde{\varphi}(y_j))e_{4d}^{\lambda_{B'}}(0, \psi'(y_j)) \\ = e_{4d}^{\lambda_A}(y_i, y_j)^{-(1-\alpha d)\alpha} = 1$$

$$e_{4d}^{\lambda_A \times \lambda_{B'}}(S_{i+g}, T_j) = e_{4d}^{\lambda_A}(0, \tilde{\varphi}(x_j))e_{4d}^{\lambda_{B'}}([\beta]\psi'(x_i), \psi'(x_j)) = e_{4d}^{\lambda_A}(x_i, x_j)^{b\beta} = 1$$

$$e_{4d}^{\lambda_A \times \lambda_{B'}}(S_{i+g}, T_{j+g}) = e_{4d}^{\lambda_A}(0, [1 - \alpha d]\tilde{\varphi}(y_j))e_{4d}^{\lambda_{B'}}([\beta]\psi'(x_i), \psi'(y_j)) \\ = e_{4d}^{\lambda_A}(x_i, y_j)^{b\beta} = e_{4d}^{\lambda_A}(x_i, y_j) = \zeta^{\delta_{i,j}},$$

for all $i, j \in [1; g]$. Hence, $(S_1, \dots, S_{2g}, T_1, \dots, T_{2g})$ is a $\zeta$-symplectic basis of $(A \times B')[4d]$. Since $([4]x_1, \dots, [4]x_g, [4]x_1, \dots, [4]y_g)$ generates $A[d]$, we clearly have by Lemma 2.3:

$$[4]\langle T_1, \dots, T_{2g} \rangle = \{([a]x, \psi' \circ \varphi(x)) \mid x \in A[d]\} = \ker(F).$$

This proves (i).

(ii) $\mathcal{B} := (U_1, \dots, U_{2g}, V_1, \dots, V_{2g})$ is the product $\mathcal{B}_1 \times \mathcal{B}_2$ with

$$\mathcal{B}_1 := ([d]x_1, \dots, [d]x_g, [d]y_1, \dots, [d]y_g)$$

$$\mathcal{B}_2 := ([d]\psi \circ \tilde{\varphi}(x_1), \dots, [d]\psi \circ \tilde{\varphi}(x_g), [d\alpha\beta]\psi \circ \tilde{\varphi}(y_1), \dots, [d\alpha\beta]\psi \circ \tilde{\varphi}(y_g)).$$

And we have for all $i, j \in [1; g]$:

$$e_4^{\lambda_B}([d]x_i, [d]x_j) = e_{4d}^{\lambda_B}(x_i, x_j)^d = 1$$
$$e_4^{\lambda_B}([d]y_i, [d]y_j) = e_{4d}^{\lambda_B}(y_i, y_j)^d = 1$$
$$e_4^{\lambda_B}([d]x_i, [d]y_j) = e_{4d}^{\lambda_B}(x_i, y_j)^d = \zeta^{d\delta_{i,j}},$$

so $\mathcal{B}_1$ is a $\zeta^d$-symplectic basis of $B[4]$. Similarly, we verify that $\mathcal{B}_2$ is a $\zeta^d$-symplectic basis of $A'[4]$. It follows that $\mathcal{B} = \mathcal{B}_1 \times \mathcal{B}_2$ is a product $\zeta^d$-symplectic basis of $(B \times A')[4]$, as desired.

(iii) The fact that $([d]F(S_1), \cdots, [d]F(S_{2g}), F(T_1), \cdots, F(T_{2g}))$ is a $\zeta^d$-symplectic basis of $(B \times A')[4]$ follows from [19, Theorem 56]. It remains to compute for all $i \in [1; g]$:

$$[d]F(S_i) = [d]F([-\alpha]\tilde{\varphi}(y_i), 0) = ([-d\alpha]\varphi \circ \tilde{\varphi}(y_i), [d\alpha]\psi \circ \tilde{\varphi}(y_i))$$
$$= ([-da\alpha]y_i, [d\alpha]\psi \circ \tilde{\varphi}(y_i)) = -V_i + [b]V_{i+g}$$

$$[d]F(S_{i+g}) = [d]F((0, [\beta]\psi'(x_i))) = ([d\beta]\tilde{\psi}' \circ \psi'(x_i), [d\beta]\tilde{\varphi}' \circ \psi'(x_i))$$
$$= ([d\beta b]x_i, [d\beta]\psi \circ \tilde{\varphi}(x_i)) = ([d]x_i, [d\beta]\psi \circ \tilde{\varphi}(x_i)) = U_i + [\beta]U_{i+g}$$

$$F(T_i) = F(\tilde{\varphi}(x_i), \psi'(x_i)) = (\varphi \circ \tilde{\varphi}(x_i) + \tilde{\psi}' \circ \psi'(x_i), -\psi \circ \tilde{\varphi}(x_i) + \tilde{\varphi}' \circ \psi'(x_i))$$
$$= ([a+b]x_i, -\psi \circ \tilde{\varphi}(x_i) + \psi \circ \tilde{\varphi}(x_i)) = ([d]x_i, 0) = U_i$$

$$F(T_{i+g}) = F([1 - \alpha d]\tilde{\varphi}(y_i), \psi'(y_i))$$
$$= ([1 - \alpha d]\varphi \circ \tilde{\varphi}(y_i) + \tilde{\psi}' \circ \psi'(y_i), -[1 - \alpha d]\psi \circ \tilde{\varphi}(y_i) + \tilde{\varphi}' \circ \psi'(y_i))$$
$$= ([a(1 - \alpha d) + b]y_i, [\alpha d - 1 + 1]\psi \circ \tilde{\varphi}(y_i)) = (0, [d]\psi \circ \tilde{\varphi}(y_i)) = [b]V_{i+g},$$

where we used twice the fact that $\psi' \circ \varphi = \varphi' \circ \psi$ implies that $\tilde{\varphi}' \circ \psi' = \psi \circ \tilde{\varphi}$. The inverse equations:

$$U_i = F(T_i), \quad U_{i+g} = [bd]F(S_{i+g}) - [b]F(T_i),$$

$$V_i = F(T_{i+g}) - [d]F(S_i), \quad V_{i+g} = [\beta]F(T_{i+g}),$$

follow immediately. This completes the proof. $\square$

## B Implementing the 4-dimensional isogeny computation for the CSIDH group action in the simplest case

In this Section, we explain how to compute the 4-dimensional isogeny $F: E_u^2 \times E_v^2 \to E_a^2 \times E_a'^2$ from Section 5.3. We compute $F$ as a chain of $e$-2-isogenies, using the algorithms from [18, § 4]. The general approach is very similar to SQIsignHD verification and SIDH attacks described in [18, § 5 and Appendix B].
## B.1 Where are the gluings and splittings?

Because the formulas for gluing isogenies are different, we need to handle gluings with special care and locate them in advance in the 2-isogeny chain. We have the following result, very similar to [18, Lemma 23]:

**Lemma B.1.** *Without loss of generality, we can assume that $2|y_u$ and $2|y_v$. Let $m := v_2(x_v y_u - x_u y_v)$. Then $F$ is of the form $F = G \circ F_{m+1} \circ F_m \circ \cdots \circ F_1$, where*

$$E_u^2 \times E_v^2 \xrightarrow{F_1} A_1^2 \xrightarrow{F_2} \cdots \xrightarrow{F_{m-1}} A_{m-1}^2 \xrightarrow{F_m} A_m^2 \xrightarrow{F_{m+1}} B$$

is a chain of 2-isogenies, the $A_i$ are principally polarised abelian surfaces and $B$ is a principally polarised abelian variety of dimension 4. For all, $i \in [2; m]$, $F_i$ is a diagonal isogeny $\mathrm{Diag}(f_i, f_i)$ with $f_i: A_{i-1} \to A_i$, a 2-dimensional 2-isogeny and $F_1: (R_1, S_1, R_2, S_2) \mapsto (f_1(R_1, R_2), f_1(S_1, S_2))$ with $f_1: E_u \times E_v \to A_1$, a gluing 2-isogeny.

*Besides, the 2-dimensional $2^m$-isogeny $\Phi := f_m \circ \cdots \circ f_1$ has kernel:*

$$\ker(\Phi) = \{([N_1 x_u] \varphi_u(P), [g_u(x_u x_v + y_u y_v)] \varphi_v \circ \hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k}(P)) \mid P \in E_1[2^m]\}.$$

*Proof.* By Eq. (15), we have:

$$\ker(F) = \{([uN_1]\varphi_u(P), [uN_1]\varphi_u(Q), [g_u]M_v \cdot M_u^T(\varphi_v \circ \hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k}(P), \varphi_v \circ \hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k}(Q)) \mid P, Q \in E_u[2^e]\},$$

so the first $m$ isogenies of the 2-isogeny chain have kernel:

$$\begin{align*}
\ker(F_m \circ \cdots \circ F_1) &= [2^{e-m}] \ker(F) \\
&= \{([uN_1]\varphi_u(P), [uN_1]\varphi_u(Q), \\
&\qquad [g_u(x_u x_v + y_u y_v)]\varphi_v \circ \hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k}(P) \\
&\qquad + [g_u(x_v y_u - x_u y_v)]\varphi_v \circ \hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k}(Q), \\
&\qquad - [g_u(x_v y_u - x_u y_v)]\varphi_v \circ \hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k}(P) \\
&\qquad + [g_u(x_u x_v + y_u y_v)]\varphi_v \circ \hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k}(Q)) \mid P, Q \in E_1[2^m]\} \\
&= S(K \times K),
\end{align*}$$

with $S: E_u^2 \times E_v^2 \to (E_u \times E_v)^2$, $(R_1, S_1, R_2, S_2) \mapsto (R_1, R_2, S_1, S_2)$ and

$$K := \{([N_1 x_u] \varphi_u(P), [g_u(x_u x_v + y_u y_v)] \varphi_v \circ \hat{\varphi}_{\mathfrak{c}_k} \circ \varphi_{\mathfrak{b}_k}(P)) \mid P \in E_1[2^m]\}.$$

It follows that $(F_m \circ \cdots \circ F_1) \circ S^{-1}$ is a diagonal $2^m$-isogeny

$$\mathrm{Diag}(\Phi, \Phi) : (E_u \times E_v)^2 \to A_m^2$$

with $\ker(\Phi) = K$, so the $F_i$ and $F$ have the desired form. $\square$

<figure><img src="image_6.png" /></figure>
Locating splittings in the isogeny chain computation is also important. Indeed, on abelian varieties 2-isogenous to a product of abelian varieties, the generic doubling procedure [18, Algorithm 8] may fail because of zero divisions. Doublings are involved to obtain kernel information at every step of the isogeny chain but we have some flexibility on the abelian variety where these doublings take place along the chain, depending on the computational strategy we use. As in [18, Appendix E], we use computational strategies that do not involve doublings on codomains of splitting isogenies. Using the structure of the polarised dual $\tilde{F}$, we obtain the

**Lemma B.2.** Assume that $2|y_u, y_v$ and let $m := v_2(x_v y_u - x_u y_v)$ as in Lemma B.1. Then $F$ is of the form $F = F_e \circ \cdots F_{e-m} \circ G'$, where

$$B' \xrightarrow{F_{e-m}} A_{e-m}^2 \xrightarrow{F_{e-m+1}} A_{e-m+1}^2 \xrightarrow{F_{e-m+2}} \cdots \xrightarrow{F_{e-m+2}} A_{e-1}^2 \xrightarrow{F_e} E_a^2 \times (E')^2$$

is a chain of 2-isogenies, the $A_i$ are principally polarised abelian surfaces and $B'$ is a principally polarised abelian variety of dimension 4. For all, $i \in [e - m + 1; e]$, $F_i$ is a diagonal isogeny $\text{Diag}(f_i, f_i)$ with $f_i : A_{i-1} \to A_i$, a 2-dimensional 2-isogeny and $F_e := S' \circ \text{Diag}(f_e, f_e)$ with $f_1 : A_{e-1} \to E_a \times E'$, a splitting 2-isogeny and $S' : (E_a \times E')^2 \to E_a^2 \times E'^2$, $(R_1, R_2, S_1, S_2) \mapsto (R_1, S_1, R_2, S_2)$.

*Proof.* The polarised dual $\tilde{F} : E_a^2 \times (E')^2 \to E_u^2 \times E_v^2$ is given by:

$$\tilde{F} = \begin{pmatrix} \Phi_u \circ \tilde{\Phi}_1 & -\Phi'_v \circ \tilde{\Phi}'_2 \\ -\Phi_v \circ \tilde{\Phi}_2 & \Phi'_u \circ \tilde{\Phi}'_1 \end{pmatrix}$$

and its kernel is:

$$
\begin{align*}
\ker(\tilde{F}) &= \{([uN_1]P, [uN_1]Q, -\tilde{\Phi}'_2 \circ \tilde{\Phi}'_v \circ \Phi_u \circ \tilde{\Phi}_1(P, Q)) \mid P, Q \in E_a[2^e]\} \\
&= \{([uN_1]P, [uN_1]Q, -M_v^T \cdot M_u(\tilde{\varphi}'_2 \circ \tilde{\varphi}'_v \circ \varphi_u \circ \tilde{\varphi}_1(P), \\
& \qquad \tilde{\varphi}'_2 \circ \tilde{\varphi}'_v \circ \varphi_u \circ \tilde{\varphi}_1(Q))) \mid P, Q \in E_a[2^e]\} \\
&= \{([uN_1]\varphi_u(P), [uN_1]\varphi_u(Q), \\
& \qquad - [x_u x_v + y_u y_v]\tilde{\varphi}'_2 \circ \tilde{\varphi}'_v \circ \varphi_u \circ \tilde{\varphi}_1(P) \\
& \qquad - [x_u y_v - x_v y_u]\tilde{\varphi}'_2 \circ \tilde{\varphi}'_v \circ \varphi_u \circ \tilde{\varphi}_1(Q), \\
& \qquad [x_u y_v - x_v y_u]\tilde{\varphi}'_2 \circ \tilde{\varphi}'_v \circ \varphi_u \circ \tilde{\varphi}_1(P) \\
& \qquad - [x_u x_v + y_u y_v]\tilde{\varphi}'_2 \circ \tilde{\varphi}'_v \circ \varphi_u \circ \tilde{\varphi}_1(Q)) \mid P, Q \in E_a[2^e]\}
\end{align*}
$$

It follows that

$$\ker(\tilde{F})[2^m] = S'(K' \times K'),$$

with:

$$K' = \{([uN_1]P, -[x_u x_v + y_u y_v]\tilde{\varphi}'_2 \circ \tilde{\varphi}'_v \circ \varphi_u \circ \tilde{\varphi}_1(P)) \mid P \in E_a[2^m]\}.$$

We then easily conclude that $F$ has the desired form. $\square$

<figure><img src="image_7.png" /></figure>
By Lemma B.1, to compute $F$, we first have to compute a 2-dimensional 2-isogeny chain $\Phi := f_m \circ \dots \circ f_1 : E_u \times E_v \to A_m$, then a 4-dimensional gluing of abelian surfaces $F_{m+1} : A_m^2 \to B$ and $e-m-1$ generic 4-dimensional 2-isogenies. Except for the first $m$ steps and splitting at the end, we do not have a way of detecting products of abelian varieties located elsewhere in the 2-isogeny chain. However, they can appear with very low probability $\tilde{O}(1/p)$ so we can assume there are none and treat the $e-m-1$ last 2-isogenies as generic ones.

Hence, we proceed as follows to compute $F$:

1. We compute the 2-dimensional isogeny $\Phi := f_m \circ \dots \circ f_1 : E_u \times E_v \to A_m$ (Appendix B.2).

2. We compute the gluing isogeny $F_{m+1} : A_m^2 \to B$ (Appendix B.3).

3. We compute the last $e-m-1$ 2-isogenies of the chain as generic isogenies with a computational strategy avoiding doublings $m$ steps before the end, as in [18, Algorithm 21].

4. We compute a change of theta-coordinates on the codomain $E_a^2 \times E'^2$ to obtain the product theta-structure and we retrieve $E_a$ (Appendix B.4).

In the following, we assume we are given the following inputs defined in Section 5.3:

- The points $P_u, Q_u, P_v, Q_v$;

- The integers $N_1, N_2, g_u, x_u, y_u, g_v, x_v, y_v, e$ obtained from a solution of Eq. (11);

where $(P_u, Q_u)$ and $(P_v, Q_v)$ are images of a basis $(P, Q)$ of $E_1[2^{e+2}]$ via the $g_u$-isogeny $\varphi_u : E_1 \to E_u$ and the $g_v N_1 N_2$-isogeny $\varphi_v \circ \hat{\varphi}_{\epsilon_k} \circ \varphi_{b_k} : E_1 \to E_v$. In the following, we denote $\zeta := e_{2^{e+2}}(P, Q)$, so that $e_{2^{e+2}}(P_u, Q_u) = \zeta^{g_u}$ and $e_{2^{e+2}}(P_v, Q_v) = \zeta^{g_v N_1 N_2}$.

## B.2 The 2-dimensional part

When computing the gluing isogeny $f_1 : E_u \times E_v \to A_1$, we first convert the Montgomery $(x:z)$-coordinates on $E_u \times E_v$ into product theta coordinates. Let $(R_u, S_u)$ and $(R_v, S_v)$ be basis of $E_u[4]$ and $E_v[4]$ respectively such that $x(S_u) = x(S_v) = -1$ and $e_4(R_u, S_u) = e_4(R_v, S_v) = \zeta^{2^e}$. By Appendix A.1, we can then convert Montgomery $(x:z)$-coordinates on $E_u$ and $E_v$ into level 2 theta coordinates associated to symmetric theta structures induced by $(R_u, S_u)$ and $(R_v, S_v)$ respectively. Then, the resulting product level 2 theta structure on $E_u \times E_v$ is induced by the $\zeta^{2^e}$-symplectic basis of $(E_u \times E_v)[4]$:

$$ \mathcal{B}_0 := ((R_u, 0), (0, R_v), (S_u, 0), (0, S_v)). $$

However, $\mathcal{B}_0$ is not adapted to $f_1$ in the sense of Appendix A.1, so we need to find a (non-product) symplectic basis adapted to $f_1$ and to compute the associated change of theta coordinates with Theorem A.1. We can then compute the gluing isogeny $f_1$ in level 2 theta coordinates using [20, Algorithm 8] and generic isogenies $f_2, \dots, f_m$ using [20, Algorithm 7].
Consider the intermediate $\zeta^{2^e}$-symplectic basis $\mathcal{B}_1$ of $(E_u \times E_v)[4]$, where:

$$ \mathcal{B}_1 := (([2^e]P_u, 0), (0, [2^e]P_v), ([2^e\gamma_u]Q_u, 0), (0, [2^e\gamma_v\eta_1\eta_2]Q_v)), $$

and $\gamma_u, \gamma_v, \eta_1, \eta_2$ are inverses of $g_u, g_v, N_1, N_2$ modulo $2^{e+2}$. Given the expression of the kernel of $\Phi := f_m \circ \dots \circ f_1$ from Lemma B.1, we can easily express the symplectic change of basis matrix $M_1$ from $\mathcal{B}_1$ to $[2^m]\mathcal{B}_2$, where $\mathcal{B}_2$ is a $\zeta^{2^{e-m}}$-symplectic basis of $(E_u \times E_v)[2^{m+2}]$ adapted to $\Phi$ (given by Lemma B.3). This matrix $M_1$ is described in Lemma B.4. We can also express the change of basis from $\mathcal{B}_0$ to $\mathcal{B}_1$:

$$ \begin{pmatrix} a_u & 0 & c_u & 0 \\ 0 & a_v & 0 & c_v \\ b_u & 0 & d_u & 0 \\ 0 & b_v & 0 & d_v \end{pmatrix}, $$

where $[2^e]P_u = [a_u]R_u + [b_u]S_u$, $[2^e\gamma_u]Q_u = [c_u]R_u + [d_u]S_u$, $[2^e]P_v = [a_v]R_v + [b_v]S_v$ and $[2^e\gamma_v]Q_v = [c_v]R_v + [d_v]S_v$. The change of basis from $\mathcal{B}_0$ to $[2^m]\mathcal{B}_2$ is then $M_2 := M_0 \cdot M_1$ to which we can directly apply Theorem A.1 to obtain the desired change of theta coordinates.

**Lemma B.3.** We keep the notations from Appendix B.1 and Lemma B.1. Let $\zeta := e_{2^{e+1}}(P, Q)$, $s_{uv} := x_ux_v + y_uy_v$, $\delta_{uv} := x_vy_u - x_uy_v$, $\eta_1, \eta_2, \gamma_u, \gamma_v, \mu, \sigma_{uv}$ be inverses of $N_1, N_2, g_u, g_v, u = g_u(x_u^2 + y_u^2)$, $s_{uv}$ modulo $2^{m+2}$ respectively. Then, $\mathcal{B}_2 := (J_1, J_2, K_1, K_2)$ given by:

$$ J_1 := (-[2^{e-m}\eta_1\mu\gamma_u]Q_u, 0), \quad J_2 := (0, [2^{e-m}\gamma_u\gamma_v\eta_1\eta_2\sigma_{uv}]P_v) $$

$$ K_1 := ([2^{e-m}N_1u]P_u, [2^{e-m}g_u s_{uv}]P_v), $$

$$ K_2 := ([2^{e-m}(N_1u + g_u g_v \mu N_2 \delta_{uv}^2)]Q_u, [2^{e-m}g_u s_{uv}]Q_v) $$

is a $\zeta^{2^{e-m}}$-symplectic basis of $(E_u \times E_v)[2^{m+2}]$ adapted to $\Phi := f_m \circ \dots \circ f_1 : E_u \times E_v \to A_m$ i.e. such that $\ker(\Phi) = \langle [4]K_1, [4]K_2 \rangle$.

*Proof.* By Lemma B.1, we have:

$$ \ker(\Phi) = \langle ([2^{e+2-m}N_1u]P_u, [2^{e+2-m}g_u s_{uv}]P_v), ([2^{e+2-m}N_1u]Q_u, [2^{e+2-m}g_u s_{uv}]Q_v) \rangle, $$

so that $\ker(\Phi) = \langle [4]K_1, [4]K_2 \rangle$ because $2^{e+2-m}\delta_{uv}^2 \equiv 0 \mod 2^{e+2}$. Besides, we clearly have

$$ e_{2^{m+2}}(J_1, J_2) = e_{2^{m+2}}(J_1, K_2) = e_{2^{m+2}}(J_2, K_1) = 1 $$

and

$$ \begin{aligned} e_{2^{m+2}}(K_1, K_2) &= e_{2^{m+2}}([2^{e-m}N_1u]P_u, [2^{e-m}g_u s_{uv}]P_v), \\ & \quad ([2^{e-m}(N_1u + g_u g_v \mu N_2 \delta_{uv}^2)]Q_u, [2^{e-m}g_u s_{uv}]Q_v)) \\ &= e_{2^{m+2}}([2^{e-m}]P_u, [2^{e-m}]Q_u)^{N_1u(N_1u + g_u g_v \mu N_2 \delta_{uv}^2)} \end{aligned} $$
$$
\begin{align*}
& \cdot e_{2m+2}([2^{e-m}]P_v, [2^{e-m}]Q_v)g_u^2 s_{uv}^2 \\
&= e_{2e+2}(P_u, Q_u)2^{e-m}N_1 u (N_1 u + g_u g_v \mu N_2 \delta_{uv}^2) \\
& \cdot e_{2e+2}(P_v, Q_v)2^{e-m}g_u^2 s_{uv}^2 \\
&= e_{2e+2}(P, Q)2^{e-m}(g_u N_1 u (N_1 u + g_u g_v \mu N_2 \delta_{uv}^2) + g_v N_1 N_2 g_u^2 s_{uv}^2) \\
&= \zeta^{2^{e-m}}(g_u N_1^2 u^2 + g_v N_1 N_2 g_u^2 (\delta_{uv}^2 + s_{uv}^2)) \\
&= \zeta^{2^{e-m}}(g_u N_1^2 u^2 + g_v N_1 N_2 g_u^2 (x_u^2 + y_u^2) (x_v^2 + y_v^2)) \\
&= \zeta^{2^{e-m}}g_u u N_1 (u N_1 + v N_2) = \zeta^{2^{2e-m}} = 1
\end{align*}
$$

Indeed, Eq. (11) implies that $2^m \le |\delta_{uv}| \le x_u^2 + y_u^2 + x_v^2 + y_v^2 \le u + v = O(2^{e/2})$, so that $m = O(e/2)$ and $2^{2e-m} \equiv 0 \mod 2^{e+2}$. Finally,

$$
\begin{align*}
e_{2m+2}(J_1, K_1) &= e_{2m+2}((-2^{e-m}\eta_1\mu\gamma_u)Q_u, 0), ([2^{e-m}N_1u]P_u, [2^{e-m}g_u s_{uv}]P_v)) \\
&= e_{2e+2}(Q_u, P_u)^{-2^{e-m}\eta_1\mu\gamma_u N_1 u} = e_{2e+2}(P_u, Q_u)^{2^{e-m}\gamma_u} \\
&= e_{2e+2}(P, Q)^{2^{e-m}\gamma_u g_u} = \zeta^{2^{e-m}},
\end{align*}
$$

and

$$
\begin{align*}
e_{2m+2}(J_2, K_2) &= e_{2m+2}((0, [2^{e-m}\gamma_u\gamma_v\eta_1\eta_2\sigma_{uv}]P_v), \\
&\quad ([2^{e-m}(N_1u + g_u g_v \mu N_2 \delta_{uv}^2)]Q_u, [2^{e-m}g_u s_{uv}]Q_v)) \\
&= e_{2e+2}(P_v, Q_v)^{2^{e-m}\gamma_u\gamma_v\eta_1\eta_2\sigma_{uv}g_u s_{uv}} \\
&= e_{2e+2}(P, Q)^{2^{e-m}g_v N_1 N_2 \gamma_u \gamma_v \eta_1 \eta_2 \sigma_{uv} g_u s_{uv}} \\
&= \zeta^{2^{e-m}},
\end{align*}
$$

so $(J_1, J_2, K_1, K_2)$ is indeed a $\zeta^{2^{e-m}}$-symplectic basis of $(E_u \times E_v)[2^{m+2}]$. $\square$

**Lemma B.4.** *The change of basis matrix from $\mathcal{B}_1$ to $[2^m]\mathcal{B}_2$ is:*

$$
M_1 := \begin{pmatrix}
0 & 0 & N_1 u & 0 \\
0 & \gamma_u \gamma_v \eta_1 \eta_2 \sigma_{uv} & g_u s_{uv} & 0 \\
-\eta_1 \mu & 0 & 0 & g_u N_1 u \\
0 & 0 & 0 & g_u g_v N_1 N_2 s_{uv}
\end{pmatrix}.
$$

*Proof.* It follows immediately from Lemma B.3 and the fact that $m \ge 1$ so that $4|\delta_{uv}^2$. $\square$

## B.3 The 4-dimensional gluing

Let $H := F_m \circ \cdots \circ F_1 : E_u^2 \times E_v^2 \to A_m^2$ as in Lemma B.1. Since $H = \text{Diag}(\Phi, \Phi) \circ S$ with $S : (R_1, S_1, R_2, S_2) \mapsto (R_1, R_2, S_1, S_2)$, if we assume we have computed $\Phi$ in Appendix B.2, we can evaluate $H$ on any point. Images of $H$ are expressed
in product level 2 theta coordinates on $A_m^2$, where the theta structure on $A_m$ is associated to the $\zeta^{2^e}$-symplectic basis of $A_m$[4] induced by $\Phi$, and given by:

$$ \mathcal{C}_0 := ([2^m]\Phi(J_1), [2^m]\Phi(J_2), \Phi(K_1), \Phi(K_2)). $$

Hence, the product theta structure on $A_m^2$ is induced by $\mathcal{C}_0 \times \mathcal{C}_0$, as defined in Eq. (16).

In general, this basis is not adapted to $F_{m+1}: A_m^2 \to B$. Such a basis is given by:

$$ \mathcal{C}_1 := ([2^e]H(S_1), \dots, [2^e]H(S_4), [2^{e-m}]H(T_1), \dots, [2^{e-m}]H(T_4)), $$

where $\mathcal{C} := (S_1, \dots, S_4, T_1, \dots, T_4)$ is a $\zeta$-symplectic basis of $(E_u^2 \times E_v^2)[2^{e+2}]$ adapted to $F$. Applying Lemma A.1 to $(x_1, x_2, y_1, y_2) := ((\varphi_{\mathfrak{b}_k}(P), 0), (0, \varphi_{\mathfrak{b}_k}(P)), ([\eta_1]\varphi_{\mathfrak{b}_k}(Q), 0), (0, [\eta_1]\varphi_{\mathfrak{b}_k}(Q)))$, we obtain:

$$ S_1 := ([-\mu\eta_1 x_u]Q_u, [-\mu\eta_1 y_u]Q_u, 0, 0), \quad S_2 := ([\mu\eta_1 y_u]Q_u, [-\mu\eta_1 x_u]Q_u, 0, 0), $$

$$ S_3 := (0, 0, [\nu\eta_2 x_v]P_v, [\nu\eta_2 y_v]P_v), \quad S_4 := (0, 0, [-\nu\eta_2 y_v]P_v, [\nu\eta_2 x_v]P_v), $$

$$ T_1 := ([N_1 x_u]P_u, [N_1 y_u]P_u, [x_v]P_v, [y_v]P_v), $$

$$ T_2 := ([-N_1 y_u]P_u, [N_1 x_u]P_u, [-y_v]P_v, [x_v]P_v), $$

$$ T_3 := ([1 - 2^e \mu\eta_1)x_u]Q_u, [(1 - 2^e \mu\eta_1)y_u]Q_u, [\eta_1 x_v]Q_v, [\eta_1 y_v]Q_v), $$

$$ T_4 := ([-(1 - 2^e \mu\eta_1)y_u]Q_u, [(1 - 2^e \mu\eta_1)x_u]Q_u, [-\eta_1 y_v]Q_v, [\eta_1 x_v]Q_v), $$

where $\mu, \nu, \eta_1, \eta_2$ are inverses of $u, v, N_1, N_2$ modulo $2^{e+2}$ respectively.

To compute $F_{m+1}$, we need to compute the symplectic change of basis from $\mathcal{C}_0 \times \mathcal{C}_0$ to $\mathcal{C}_1$ and the associated change of theta coordinates. We can then apply the algorithms from [18, § 4.3] to compute $F_{m+1}$ with level 2 theta coordinates.

**Lemma B.5.** *The change of basis matrix from $\mathcal{C}_0 \times \mathcal{C}_0$ to $\mathcal{C}_1$ is*

$$ M_3 := \begin{pmatrix} g_u x_u - g_u y_u & 0 & 0 & 0 & 0 & \mu_2 & \mu_3 \\ 0 & 0 & \lambda x_v - \lambda y_v & \mu_1 y_u & \mu_1 x_u & 0 & 0 \\ g_u y_u & g_u x_u & 0 & 0 & 0 & -\mu_3 & \mu_2 \\ 0 & 0 & \lambda y_v & \lambda x_v & -\mu_1 x_u & \mu_1 y_u & 0 \\ 0 & 0 & 0 & 0 & \mu x_u & -\mu y_u & 0 \\ 0 & 0 & 0 & 0 & 0 & \eta_1 x_v \gamma_u \sigma_{uv} & -\eta_1 y_v \gamma_u \sigma_{uv} \\ 0 & 0 & 0 & 0 & \mu y_u & \mu x_u & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \eta_1 y_v \gamma_u \sigma_{uv} & \eta_1 x_v \gamma_u \sigma_{uv} \end{pmatrix}, $$

where $\lambda := \nu g_u g_v N_1 s_{uv}$, $\mu_1 := \mu g_u^2 g_v s_{uv} N_1 N_2 \delta_{uv} / 2^m$,

$$ \mu_2 := \frac{\delta_{uv}}{2^m} g_u \sigma_{uv} (N_1 u y_u + g_v x_v N_2 \delta_{uv}) \quad \text{and} $$

$$ \mu_3 := \frac{\delta_{uv}}{2^m} g_u \sigma_{uv} (N_1 u x_u - g_v y_v N_2 \delta_{uv}). $$
*Proof.* We have:

$$ 
\begin{aligned}
[2^e]H(S_1) &= (\Phi([-2^e \mu \eta_1 x_u]Q_u, 0), \Phi([-2^e \mu \eta_1 y_u]Q_u, 0)) \\
&= ([g_u x_u]\Phi([-2^e \eta_1 \mu \gamma_u]Q_u, 0), [g_u y_u]\Phi([-2^e \eta_1 \mu \gamma_u]Q_u, 0)) \\
&= ([g_u x_u][2^m]\Phi(J_1), [g_u y_u][2^m]\Phi(J_1))
\end{aligned}
$$

$$ 
\begin{aligned}
[2^e]H(S_2) &= (\Phi([2^e \mu \eta_1 y_u]Q_u, 0), \Phi([-2^e \mu \eta_1 x_u]Q_u, 0)) \\
&= ([g_u y_u][2^m]\Phi(J_1), [g_u x_u][2^m]\Phi(J_1))
\end{aligned}
$$

$$ 
\begin{aligned}
[2^e]H(S_3) &= (\Phi(0, [2^e \nu \eta_2 x_v]P_v), \Phi(0, [2^e \nu \eta_2 y_v]P_v)) \\
&= ([\nu x_v g_u g_v N_1 s_{uv}]\Phi(0, [2^e \gamma_u \gamma_v \eta_1 \eta_2 \sigma_{uv}]P_v), \\
&\qquad [\nu y_v g_u g_v N_1 s_{uv}]\Phi(0, [2^e \gamma_u \gamma_v \eta_1 \eta_2 \sigma_{uv}]P_v)) \\
&= ([\lambda x_v][2^m]\Phi(J_2), [\lambda y_v][2^m]\Phi(J_2))
\end{aligned}
$$

$$ 
\begin{aligned}
[2^e]H(S_4) &= (\Phi(0, [-2^e \nu \eta_2 y_v]P_v), \Phi(0, [2^e \nu \eta_2 x_v]P_v)) \\
&= ([\lambda y_v][2^m]\Phi(J_2), [\lambda x_v][2^m]\Phi(J_2)).
\end{aligned}
$$

We also have:

$$ [2^{e-m}]H(T_1) = (\Phi([2^{e-m}N_1x_u]P_u, [2^{e-m}x_v]P_v), \Phi([2^{e-m}N_1y_u]P_u, [2^{e-m}y_v]P_v)) $$

with:

$$ 
\begin{aligned}
\Phi([2^{e-m}N_1x_u]P_u, [2^{e-m}x_v]P_v) &= [\mu x_u]\Phi([2^{e-m}N_1u]P_u, [2^{e-m}g_u s_{uv}]P_v) \\
&\quad + \Phi(0, [2^{e-m}(x_v - \mu x_u g_u s_{uv})]P_v) \\
&= [\mu x_u]\Phi(K_1) + \Phi(0, [2^{e-m}\mu g_u y_u(x_v y_u - x_u y_v)]P_v) \\
&= [\mu x_u]\Phi(K_1) \\
&\quad + [\mu g_u^2 g_v y_u \delta_{uv} s_{uv} N_1 N_2]\Phi(0, [2^{e-m}\gamma_u \gamma_v \eta_1 \eta_2 \sigma_{uv}]P_v) \\
&= [\mu_1 y_u][2^m]\Phi(J_2) + [\mu x_u]\Phi(K_1)
\end{aligned}
$$

and:

$$ 
\begin{aligned}
\Phi([2^{e-m}N_1y_u]P_u, [2^{e-m}y_v]P_v) &= [\mu y_u]\Phi(K_1) \\
&\quad + \Phi(0, [2^{e-m}(y_v - \mu y_u g_u s_{uv})]P_v) \\
&= [\mu y_u]\Phi(K_1) + \Phi(0, [-2^{e-m}\mu g_u x_u(x_v y_u - x_u y_v)]P_v) \\
&= [\mu y_u]\Phi(K_1) \\
&\quad - [\mu g_u^2 g_v x_u \delta_{uv} s_{uv} N_1 N_2]\Phi(0, [2^{e-m}\gamma_u \gamma_v \eta_1 \eta_2 \sigma_{uv}]P_v) \\
&= -[\mu_1 x_u][2^m]\Phi(J_2) + [\mu y_u]\Phi(K_1),
\end{aligned}
$$

so that:

$$ [2^{e-m}]H(T_1) = ([\mu_1 y_u][2^m]\Phi(J_2) + [\mu x_u]\Phi(K_1), -[\mu_1 x_u][2^m]\Phi(J_2) + [\mu y_u]\Phi(K_1)) $$
And we also have

$$ [2^{e-m}]H(T_2) = (\Phi([-2^{e-m}N_1y_u]P_u, [-2^{e-m}y_v]P_v), \Phi([2^{e-m}N_1x_u]P_u, [2^{e-m}x_v]P_v)) \\ = ([\mu_1x_u][2^m]\Phi(J_2) - [\mu y_u]\Phi(K_1), [\mu_1y_u][2^m]\Phi(J_2) + [\mu x_u]\Phi(K_1)) $$

Similarly, we compute

$$ [2^{e-m}]H(T_3) = (T_{3,1}, T_{3,2}), $$

with:

$$ \begin{aligned} T_{3,1} &:= \Phi([2^{e-m}(1 - 2^e \mu \eta_1)x_u]Q_u, [2^{e-m}\eta_1 x_v]Q_v) \\ &= \Phi([2^{e-m}x_u]Q_u, [2^{e-m}\eta_1 x_v]Q_v) \\ &= [\eta_1 x_v \gamma_u \sigma_{uv}] \Phi([2^{e-m}(N_1 u + g_u g_v \mu N_2 \delta_{uv}^2)]Q_u, [2^{e-m}g_u s_{uv}]Q_v) \\ &\quad + \Phi([2^{e-m}(x_u - \eta_1 x_v \gamma_u \sigma_{uv})(N_1 u + g_u g_v \mu N_2 \delta_{uv}^2)])Q_u, 0) \\ &= [\eta_1 x_v \gamma_u \sigma_{uv}] \Phi(K_2) \\ &\quad + \Phi([-2^{e-m} \delta_{uv} \eta_1 \sigma_{uv} \mu(N_1 u y_u + g_v x_v N_2 \delta_{uv})]Q_u, 0) \\ &= [\eta_1 x_v \gamma_u \sigma_{uv}] \Phi(K_2) \\ &\quad + [\delta_{uv} g_u \sigma_{uv}(N_1 u y_u + g_v x_v N_2 \delta_{uv})] \Phi([-2^{e-m} \eta_1 \mu \gamma_u]Q_u, 0) \\ &= [\mu_2][2^m] \Phi(J_1) + [\eta_1 x_v \gamma_u \sigma_{uv}] \Phi(K_2) \end{aligned} $$

and:

$$ \begin{aligned} T_{3,2} &:= \Phi([2^{e-m}(1 - 2^e \mu \eta_1)y_u]Q_u, [2^{e-m}\eta_1 y_v]Q_v) \\ &= \Phi([2^{e-m}y_u]Q_u, [2^{e-m}\eta_1 y_v]Q_v) \\ &= [\eta_1 y_v \gamma_u \sigma_{uv}] \Phi([2^{e-m}(N_1 u + g_u g_v \mu N_2 \delta_{uv}^2)]Q_u, [2^{e-m}g_u s_{uv}]Q_v) \\ &\quad + \Phi([2^{e-m}(y_u - \eta_1 y_v \gamma_u \sigma_{uv})(N_1 u + g_u g_v \mu N_2 \delta_{uv}^2)])Q_u, 0) \\ &= [\eta_1 y_v \gamma_u \sigma_{uv}] \Phi(K_2) \\ &\quad + \Phi([2^{e-m} \delta_{uv} \eta_1 \sigma_{uv} \mu(N_1 u x_u - y_v g_u g_v N_2 \delta_{uv})]Q_u, 0) \\ &= [\eta_1 y_v \gamma_u \sigma_{uv}] \Phi(K_2) \\ &\quad - [\delta_{uv} g_u \sigma_{uv}(N_1 u x_u - g_u g_v y_v N_2 \delta_{uv})] \Phi([-2^{e-m} \eta_1 \mu \gamma_u]Q_u, 0) \\ &= -[\mu_3] \Phi(J_1) + [\eta_1 y_v \gamma_u \sigma_{uv}] \Phi(K_2), \end{aligned} $$

so that:

$$ [2^{e-m}]H(T_3) = ([\mu_2][2^m]\Phi(J_1) + [\eta_1 x_v \gamma_u \sigma_{uv}]\Phi(K_2), \\ - [\mu_3][2^m]\Phi(J_1) + [\eta_1 y_v \gamma_u \sigma_{uv}]\Phi(K_2)). $$

Finally, we have

$$ \begin{aligned} [2^{e-m}]H(T_4) &= (\Phi([-1 - 2^e \mu \eta_1)y_u]Q_u, [-\eta_1 y_v]Q_v), \\ &= (\Phi([1 - 2^e \mu \eta_1)x_u]Q_u, [\eta_1 x_v]Q_v)) = (-T_{3,2}, T_{3,1}) \\ &= ([\mu_3][2^m]\Phi(J_1) - [\eta_1 y_v \gamma_u \sigma_{uv}]\Phi(K_2), \\ &\quad [\mu_2][2^m]\Phi(J_1) + [\eta_1 x_v \gamma_u \sigma_{uv}]\Phi(K_2)). \end{aligned} $$

This completes the proof. □

<figure><img src="image_8.png" /></figure>
## B.4 The splitting and identification of the result curve $E_a$

Once we have computed $F : E_u^2 \times E_v^2 \to E_a^2 \times E'^2$ as a 2-isogeny chain, the codomain we obtain does not have a product theta structure from which we can obtain the theta null point of $E_a$, using Eq. (17). The theta structure we obtain is induced by the $\zeta^{2^e}$-symplectic basis of $(E_a^2 \times E'^2)[4]$ given by:

$$ \mathcal{D}_1 := ([2^e]F(S_1), \dots, [2^e]F(S_4), F(T_1), \dots, F(T_4)), $$

while a product theta structure of $E_a^2 \times E'^2$ is induced by the $\zeta^{2^e}$-symplectic basis $\mathcal{D}_0 := (U_1, \dots, U_4, V_1, \dots, V_4)$, where:

$$ U_1 := ([2^e]\varphi_{b_k}(P), 0, 0, 0), \quad U_2 := (0, [2^e]\varphi_{b_k}(P), 0, 0), $$

$$ U_3 := (0, 0, [2^e]\Psi \circ \tilde{\Phi}(\varphi_{b_k}(P), 0)), \quad U_4 := (0, 0, [2^e]\Psi \circ \tilde{\Phi}(0, \varphi_{b_k}(P))), $$

$$ V_1 := ([2^e \eta_1] \varphi_{b_k}(Q), 0, 0, 0), \quad V_2 := (0, [2^e \eta_1] \varphi_{b_k}(Q), 0, 0), $$

$$ V_3 := (0, 0, [2^e \mu \nu \eta_1^2 \eta_2] \Psi \circ \tilde{\Phi}(\varphi_{b_k}(Q), 0)), \quad V_4 := (0, 0, [2^e \mu \nu \eta_1^2 \eta_2] \Psi \circ \tilde{\Phi}(0, \varphi_{b_k}(Q))), $$

given by Lemma A.1, where $\eta_1, \eta_2, \mu, \nu$ are modular inverse of $N_1, N_2, u, v$ modulo 4 respectively. By Lemma A.1, we also obtain the

**Lemma B.6.** The change of basis matrix from $\mathcal{D}_1$ to $\mathcal{D}_0$ is given by

$$ M_4 := \begin{pmatrix} 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & vN_2 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & vN_2 & 0 & 0 & 0 & 0 \\ 1 & -vN_2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & -vN_2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & \nu\eta_2 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & \nu\eta_2 \end{pmatrix}. $$

Hence, using the matrix $M_4$ above, the non-product theta null point of $E_a^2 \times E'^2$ and the change of coordinates formulas from [18, Theorem 12], we can compute the product theta null point of $E_a^2 \times E'^2$ given by:

$$ (\theta_{i_1, i_2, i_3, i_4}^{E_a^2 \times E'^2}(0))_{i_1, i_2, i_3, i_4 \in \mathbb{Z}/2\mathbb{Z}} = (\theta_{i_1}^{E_a}(0) \cdot \theta_{i_2}^{E_a}(0) \cdot \theta_{i_3}^{E'}(0) \cdot \theta_{i_4}^{E'}(0))_{i_1, i_2, i_3, i_4 \in \mathbb{Z}/2\mathbb{Z}}. \quad (18) $$

We can extract the projective theta null point $(a : b) := (\theta_0^{E_a}(0) : \theta_1^{E_a}(0))$ by finding $i_3, i_4 \in \mathbb{Z}/2\mathbb{Z}$ such that the $\theta_{0,0,i_3,i_4}^{E_a^2 \times E'^2}(0) \neq 0$ and set $a := \theta_{0,0,i_3,i_4}^{E_a^2 \times E'^2}(0)$ and $b := \theta_{1,0,i_3,i_4}^{E_a^2 \times E'^2}(0)$. This theta null point $(a : b)$ gives 6 different possibilities for the Montgomery coefficient $A$ of

$$ E_a : By^2 = x^3 + Ax^2 + x. $$

Those possibilities are presented in Proposition A.1. All of these possible coefficients are not defined over $\mathbb{F}_p$ as $E_a$ should. This removes 2 to 4 possibilities out of 6.
However, we still cannot distinguish $E_a$ from its quadratic twist i.e. identify $B$ from the theta null point $(a : b)$. To circumvent this difficulty, we impose $B = 1$ (or any quadratic residue over $\mathbb{F}_p$) and require that 4-torsion points of the form $(\pm 1 : * : 1)$ are defined in $E_a(\mathbb{F}_p)$ i.e. that $A+2$ and $A-2$ are quadratic residues over $\mathbb{F}_p$. This lifts all ambiguity since a Montgomery curve is determined by a point of 4-torsion.

In practice, to identify the right Montgomery coefficient $A$ of $E_a$, we pick a 4-torsion point $T_u \in E_u(\mathbb{F}_p)$, which does exist since $E_u$ is defined over $\mathbb{F}_p$, as the result of an ideal action on $E$ (see the proof of Lemma 5.1). Since

$$F(T_u, 0, 0, 0) = ([x_u]\varphi_{b_k} \circ \hat{\varphi}_u(T_u), *),$$

and both $\varphi_{b_k}$ and $\varphi_u$ are associated to ideal actions, the point $T_a := [x_u]\varphi_{b_k} \circ \hat{\varphi}_u(T_u)$ is $\mathbb{F}_p$-rational. Because $N_1 = N(a)$ and $g_u = \deg(\varphi_u)$ are odd and $x_u$ can be assumed to be odd (swapping $x_u$ and $y_u$ if necessary), $T_a$ is also a point of 4-torsion. We can then compute several change of theta coordinates until $T_a$ is mapped to one of the $\mathbb{F}_p$-rational theta points $(1:0)$ or $(0:1)$ that are mapped to the desired $\mathbb{F}_p$-rational Montgomery points $(-1:1)$ and $(1:1)$ respectively. Once this is done, we identify a Montgomery coefficient $A$ from the new theta null point, following Proposition A.1.

To lift $(\pm 1 : 1)$ in $E_a(\mathbb{F}_p)$ given by $y^2 = x^3 + Ax^2 + x$, we need $A \pm 2$ to be a quadratic residue over $\mathbb{F}_p$. If this is not the case, then $E_a$ is given by $-y^2 = x^3 + Ax^2 + x$ (in which $(\pm 1 : 1)$ is $\mathbb{F}_p$-rational). But there is an $\mathbb{F}_p$-rational isomorphism $(x, y) \mapsto (-x, y)$ from $-y^2 = x^3 + Ax^2 + x$ to $y^2 = x^3 - Ax^2 + x$ so we can change $A$ into $-A$ in this case.

In practice, on the Montgomery curve $E_u$, we take $T_u$ the point with $x$-coordinate $x(T_u) = 1$. By the same reasoning as above, replacing $A_u$ by $-A_u$ if necessary, we can assume that $T_u$ is rational. The level 2 theta null point we obtain in the end for $E_a$ implicitly contains the $x$-coordinate of the image of $T_u$, so we do not even need to push the 4-torsion point $T_u$ explicitly.

## C Relaxing the $\mathfrak{B}$-good condition

This section outlines an approach to relax the $\mathfrak{B}$-good condition of Section 4.1 to improve the probability of finding $(u, v)$ such that isogenies of degree $u$ and $v$ can be efficiently computed in dimension 2.

Instead of looking for a pair $(u, v)$ of $\mathfrak{B}$-good integers, we search for $(u, v)$ such that $u$ takes the form $u = g_u(ax_u^2 + by_u^2)$ and $v$ the form $v = g_v(ax_v^2 + by_v^2)$, where $x_u, y_u, x_v, y_v \in \mathbb{N}$ and $g_u, g_v, a, b \in \mathbb{N}^*$ are products of primes lying in $\mathfrak{B}$. Note that the $\mathfrak{B}$-good case is the case $a = b = 1$.

Before discussing the advantages and drawbacks of introducing additional freedom through $a$ and $b$, let us describe how the Kani diagram in dimension 4 for computing $E_a$ (cf. Section 3.1) is reshaped by this change. In this construction, we assume that $\mathfrak{B}$ is a set of Elkies primes. As we shall see, this assumption is crucial to ensure the well-definedness of our Kani diagram.
Moreover, since both **b** and **c** are equivalent to **a**, we have that $E_3 = [I_a I_b] \cdot E_a = E_4$. We denote by $E'$ the curve $[I_a I_b] \cdot E_a$.

Finally, using Lemma 2.2, there exist a $g_u(ax_u^2+by_u^2)N(c_k)$-isogeny $\Psi: A_u \to A$ and a $g_v(ax_v^2+by_v^2)N(b_k)$-isogeny $\Phi: A \to A_v$ such that $\Phi \circ \Psi = \Phi_v \circ \tilde{\Phi}_{c_k} \circ \Phi_{b_k} \circ \tilde{\Phi}_u$. Therefore, we have the following Kani-square

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

By Lemma 2.3 and as $g_u(ax_u^2+by_u^2)N_1 + g_v(ax_v^2+by_v^2)N_2 = 2^e$, this isogeny-diamond leads to an $2^e$-isogeny $F: A_u \times A_v \to E_a \times E' \times A$ with kernel

$$ \ker(F) = \{([g_u(ax_u^2 + by_u^2)N_1](x), \Phi_v \circ \tilde{\Phi}_{c_k} \circ \Phi_{b_k} \circ \tilde{\Phi}_u(x) | x \in A_u[2^e])\}. $$

This kernel can also be expressed from points of $E_1 \times E_{1,ab}$ as

$$ \ker(F) = \{([N_1]\Phi_u(P, Q), \Phi_v \circ \tilde{\Phi}_{c_k} \circ \Phi_{b_k}(P, Q)) | (P, Q) \in E_1 \times E_{1,ab}[2^e]\}. $$

From there, one can compute this $2^e$-isogeny $F$ using the same method as in Section 3.1.

**Assessment of the relaxation.** We now investigate the benefits and drawbacks of considering pairs $(u, v)$ of such a form. First, the main additional cost is the computation of the *a*-isogenies $\varphi_{1,a}, \varphi'_{1,a}, \varphi_{2,a}, \varphi'_{2,a}$ and *b*-isogenies $\varphi_{1,b}, \varphi'_{1,b}, \varphi_{2,b}, \varphi'_{2,b}$. Since half of them are simply pushforwards of the other half, they can essentially be computed with two calls to the Elkies' algorithm for the *a*-isogenies and two calls for the *b*-isogenies. On the other hand, one may expect that the probability of an integer being representable as $ax^2 + by^2$, with $a, b \in \{1, 2, 3, 5, 7\}$, is at least three times greater than the probability of being representable by a sum of two squares. The table of experimental results in Figure 1 supports this assumption.

Let us briefly explain why the set $\{2, 3, 5, 7\}$ is an interesting choice for $\mathfrak{B}$ regarding $a, b$.

We first focus on forms $x^2+by^2$ where $b$ is prime. Our discussion relies mainly on Proposition C.1 and on the fact that an integer whose prime factors are all represented by such a quadratic form can itself be represented by this form.

**Proposition C.1.** [17, Corollary 2.6] Let $b$ be an integer and $p$ be an odd prime such that $\gcd(b, p) = 1$. Then $p$ is represented by a primitive form of discriminant $-4b$ if and only if its Legendre symbol $(-b/p)$ is equal to 1.
<table><thead><tr><th>Bit size</th><th>32 bits</th><th>64 bits</th><th>128 bits</th></tr></thead><tbody><tr><td>Success rate when a, b ∈ {1}</td><td>0.16</td><td>0.12</td><td>0.09</td></tr><tr><td>Success rate when a, b ∈ {1, 2, 3, 7}</td><td>0.48</td><td>0.37</td><td>0.28</td></tr><tr><td>Success rate when a, b ∈ {1, 2, 3, 5, 7}</td><td>0.54</td><td>0.42</td><td>0.32</td></tr></tbody></table>

**Fig. 1.** Success rates for representing an integer as $ax^2 + by^2$, with $x, y \in \mathbb{N}$, depending on the possible values for $a$ and $b$. Each rate is computed from $10^4$ integers drawn uniformly at random. Forms where $a \neq 1$ and $b = 7$ were not considered.

Note that the quadratic form $x^2+by^2$ has discriminant $-4b$. Hence, the best choices for $b$ to maximize the probability that primes are represented by $x^2+by^2$ are those for which the class group of discriminant $-4b$ has cardinality 1. Indeed, by Proposition C.1, if $(b/p) = 1$, then $p$ is represented by a primitive form of discriminant $4b$. When there is a unique class in this group, we have that $x^2+by^2$ is equivalent to this primitive form. Therefore, $x^2+by^2$ also represents $p$.

Since an integer has a probability of one-half of being a quadratic residue modulo a prime $p$, we expect that $x^2 + by^2$ represents $p$ with probability one-half. In addition, since the set of even discriminants of class number 1 is $-4, -8, -12, -16, -28$, there are only 4 such forms interesting to us: $x^2 + y^2$, $x^2 + 2y^2$, $x^2 + 3y^2$ and $x^2 + 7y^2$. We cover all of them with the set of primes $\{2, 3, 7\}$.

The second-best choices for $b$ are those for which the corresponding class group has cardinality 2, such as $b = 5$ and $b = 13$. The probability that a prime is representable by $x^2 + by^2$ is then halved for such $b$. However, from [17], we have the following equivalences

<table><tr><td>$p = x^2 + y^2$</td><td>iff $p = 1 \mod 4$,</td></tr><tr><td>$p = x^2 + 5y^2$</td><td>iff $p = 1, 9 \mod 20$,</td></tr><tr><td>$p = x^2 + 13y^2$</td><td>iff $p = 1, 9, 17, 25, 29, 49 \mod 52$.</td></tr></table>

Hence, if $p$ is represented by $x^2 + 5y^2$ or $x^2 + 13y^2$, it is also represented by $x^2 + y^2$. Nevertheless, adding 5 to the list of prime factors gives access to the forms $2x^2 + 5x^2$, $x^2 + 10x^2$ and so on. This explains why, in Fig. 1, adding 5 improves the probability of success.

Other experiments have shown that adding a quadratic form from a class group with cardinality 3 or more, such as $x^2+11y^2$, has only a negligible impact. This is why, we did not consider forms $ax^2 + by^2$ where $a \neq 1$ and $b = 7$, nor forms where $a$ or $b$ are composite numbers. All these cases lead to class groups of cardinality greater than 3. Additionally, the next prime after 13 that yields a class group of cardinality 2 is 37. Therefore, to achieve the most significant improvement with the smallest set of primes, it seems reasonable to consider only primes up to 7.
# D Algorithms for working over $\mathbb{F}_p$

In this appendix, we give the necessary algorithms for working over $\mathbb{F}_p$.

## D.1 Fast sampling of a basis given by the eigenspaces

Let $p \equiv 7 \pmod 8$, and let $E/\mathbb{F}_p$ be a supersingular elliptic curve on the surface, i.e. oriented by $\mathbb{Z}[\frac{1+\sqrt{-p}}{2}]$. Recall that for $2^f \parallel (p+1)$, we have

$$E(\mathbb{F}_p)[2^f] = \mathbb{Z}/2\mathbb{Z} \times \mathbb{Z}/2^{f-1}\mathbb{Z},$$

and that on $E[2^{f-1}]$, the Frobenius has eigenvalues $\pm 1$. Our goal is to efficiently sample a basis $\langle P, Q \rangle = E[2^{f-1}]$, such that $x(P), x(Q) \in \mathbb{F}_p$. To do this, we first study the three non-trivial 2-torsion points on $E$. Denote them as follows:

- $T_0$ is the point corresponding to the descending 2-isogeny, i.e. $E/\langle T_0 \rangle$ is primitively oriented by $\mathbb{Z}[\sqrt{-p}]$

- $T_{-1}$ is the point corresponding to the horizontal 2-isogeny, such that there does not exist a point $T' \in E(\mathbb{F}_p)$ satisfying $[2]T' = T_{-1}$.

- $T_1$ is the point corresponding to the horizontal 2-isogeny, such that there does exist a point $T' \in E(\mathbb{F}_p)$ satisfying $[2]T' = T_1$.

Although these points could easily be distinguished by looking at their division points, we show that they can more efficiently be distinguished by their Tate-pairings.

**Lemma D.1.** With the notation above, the (reduced) Tate-pairings are given as follows:

- $e_{T,2}(T_{-1}, T_{-1}) = e_{T,2}(T_{-1}, T_0) = e_{T,2}(T_{-1}, T_1) = 1,$

- $e_{T,2}(T_0, T_{-1}) = e_{T,2}(T_0, T_0) = -1$ and $e_{T,2}(T_0, T_1) = 1,$

- $e_{T,2}(T_1, T_{-1}) = e_{T,2}(T_1, T_0) = -1$, and $e_{T,2}(T_1, T_1) = 1.$

*Proof.* Recall (see e.g. [57, 61, Lemma 5]) that for $T_i$ and $\phi_i : E \to E/\langle T_i \rangle$, we have that $e_{T,2}(T_i, R)$ is 1 if there exists a rational point $R'$ such that $\hat{\phi}_i(R') = R$, and $-1$ otherwise. The lemma can then be proved as follows:

Note first that we certainly have $e_{T,2}(T_i, T_1) = 1$, since there exists a point $T'$ such that $[2]T' = \hat{\phi}_i \circ \phi_i(T') = T_1$.

For $T_0$, the corresponding isogeny $\phi_0$ is descending, and hence the $2^f$-torsion on $E/\langle T_0 \rangle$ is cyclic. Thus, $T_0$ is not in the rational image of $\hat{\phi}_0$, and since we know that $e_{T,2}(T_0, T_1) = 1$, we get that $e_{T,2}(T_0, T_0) = e_{T,2}(T_0, T_{-1}) = -1$.

For $T_{-1}$, we have that $\ker(\hat{\phi}_{-1}) = \langle \phi_1(T_1) \rangle$. Since $f \ge 3$, there are rational points $R' \in E/\langle T_{-1} \rangle$ (of order 4) satisfying $\hat{\phi}_{-1}(R') = T_1$, and since $E/\langle T_{-1} \rangle[2]$ is $\mathbb{F}_p$-rational, the other 2-torsion points are mapped to $T_{-1}$ and $T_0$. Hence $e_{T,2}(T_{-1}, T_i) = 1$.

Finally, for $T_1$, we have that $\ker(\hat{\phi}_1) = \langle \phi_1(T_0) \rangle = \langle \phi_1(T_{-1}) \rangle$. Since $T_0$ and $T_{-1}$ are both in $E(\mathbb{F}_p) \setminus [2]E(\mathbb{F}_p)$, its clear that $e_{T,2}(T_1, T_0) = e_{T,2}(T_1, T_{-1}) = -1$. (This also follows by the non degeneracy of the Tate pairing and the fact that $e_{T,2}(T_1, T_1) = 1$.) $\square$
Thus, given a curve and the three non-trivial 2-torsion points, we can by Lemma D.1 easily detect which point corresponds to which index. Together with the following lemma, this gives a particularly fast algorithm for sampling a basis of the form that we want.

**Lemma D.2.** Let $p \equiv 7 \pmod 8$, let

$$E/\mathbb{F}_p : y^2 = g(x)$$

be an elliptic curve with $\mathrm{End}_{\mathbb{F}_p}(E) = \mathbb{Z}[\frac{1+\pi}{2}]$, and let $f = v_2(p+1)$, so that there exists a point on $E(\mathbb{F}_p)$ of exact order $2^{f-1}$. Denote by $T_0, T_{-1}, T_1$ the three non-trivial 2-torsion points, with the indexes as above. Let $x_P \in \mathbb{F}_p$. Then

- $x_P$ lifts to a point $(x_P, y_P) = P \in E(\mathbb{F}_p)$ with $2^{f-1} | \mathrm{ord}(P)$ if and only if $(x_P - x(T_{-1})) \in \mathbb{F}_p^* \setminus \mathbb{F}_p^{*2}$ and $g(x_P) \in \mathbb{F}_p^{*2}$,

- $x_P$ lifts to a point $(x_P, y_P) = P \in E^t(\mathbb{F}_p)$ with $2^{f-1} | \mathrm{ord}(P)$ if and only if $(x_P - x(T_1)) \in \mathbb{F}_p^{*2}$ and $g(x_P) \in \mathbb{F}_p^* \setminus \mathbb{F}_p^{*2}$,

where

$$E^t/\mathbb{F}_p : by^2 = g(x)$$

denotes a quadratic twist of $E$ over $\mathbb{F}_p$. (We remark that we can take $b = -1$ because $-1$ is not a square over $\mathbb{F}_p$.)

*Proof.* For the first part, note that it is obvious that $P \in E(\mathbb{F}_p) \setminus E[2]$ if and only if $g(x_P) \in \mathbb{F}_p^{*2}$. Further, let $R \in E(\mathbb{F}_p)$ be any element of exact order $2^{f-1}$. We can write $P = [u]R + [b]T_0 + S$ with $S$ in [2] $E(\mathbb{F}_p)$. Now $2^{f-1}$ divides the order of $P$ if and only if $u$ is odd. The statement then follows from

$$e_{T,2}(T_{-1}, P) = e_{T,2}(T_{-1}, R)^u e_{T,2}(T_{-1}, T_0)^b = (-1)^u,$$

where $e_{T,2}(T_{-1}, R) = -1$ by combining Lemma D.1 and the non-degeneracy of the Tate pairing.

Next, write $\tau : E \to E^t$ for any twisting isomorphism leaving $x$ fixed, and write $T_{-1}^t, T_0^t, T_1^t$ for the non-trivial 2-torsion points on $E^t$, with indexes same as before. The second part of the lemma then follows from the first after noting that $\tau(T_1) = T_{-1}^t$, and that for any $T_{-1}^t \neq Q \in E^t(\mathbb{F}_p)$, the reduced Tate pairing $e_2(T_{-1}^t, Q)$ is given by the residue class of $\frac{x(Q)-x(T_{-1}^t)}{((x-x(T_{-1}^t))/(y/x)^2)(0_E)} = (x(Q) - x(T_{-1}^t))/b$. $\square$

Note that this gives a particularly effective algorithm for sampling a basis of the form we want.

*Remark D.1 (Choosing a canonical representative for the Montgomery coefficient).* Let $E/\mathbb{F}_q$ be an elliptic curve. Recall (see e.g. [57, Example 5.12]) that if $T \in E[2](\mathbb{F}_q)$, $E$ can be put in a Montgomery form $by^2 = x^3 + Ax^2 + x$ with $T$ sent to $(0,0)$ if and only if there exists a rational cyclic group $G \subset E[4]$ of order
4 that contains $T$. Equivalently, $T$ has trivial self Tate pairing. If that is the case, there are always exactly two different such subgroups, which give the two coefficients $A, -A$. As explained in Appendix B.4, since $-1$ is not a square over $\mathbb{F}_p$, then if $b = -1$, switching to $-A$ allows us to always work with Montgomery models with $b = 1$.

There are always 6 possible Montgomery models over $\bar{\mathbb{F}}_p$ (two for each of the three 2-torsion points). However, in our situation of $E/\mathbb{F}_p$ on the surface, only $T_1$ and $T_{-1}$ give a rational Montgomery model. Indeed, $T_0$ induces the descending isogeny to the floor, which cannot be extended further to a rational 4-isogeny. Imposing the condition that $b = 1$, we only have two possible Montgomery coefficients for our elliptic curves.

In our implementation, we pick a deterministic one based on the lexicographic order. An alternative choice would be to (for instance) pick the choice that sends $T_{-1}$ to $(0,0)$. A further optimisation would be to represent the Montgomery curve not by the $A$ coefficient, but by the $x$-coordinate $\alpha$ of $T_0$ (for instance), from which we recover $A$ as $-\alpha - 1/\alpha$. With this convention, we would be able to pick up immediately the three points $T_1, T_0, T_{-1}$ without a square root computation. For the class group action, since we first compute the theta coordinate of $E_\alpha$, which encodes all the 2-torsion on $E_\alpha$, we can easily output $\alpha_\alpha$ rather than $A_\alpha$. Explicitly, if $(a : b)$ is a rational theta null point on $E_\alpha$, then $\alpha_\alpha = (b^2 - a^2)/(a^2 + b^2)$ is the $x$-coordinate of $T_0$ when $E_\alpha$ is put into the Montgomery form such that the 2-torsion theta point $(-a : b)$ (which has trivial self pairing, hence is either $T_1$ or $T_{-1}$) is sent to $(0 : 1)$.

## D.2 Evaluating endomorphisms with denominator

In this section, we show how to evaluate endomorphisms with denominator, by using Tate-pairings instead of division points. We first give the general algorithm, before noting that our situation is particularly simple.

<caption>**Algorithm 4 EvalEndomorphismWithDenominator($\alpha, P$)**</caption>

**Input:** $\alpha \in \mathbb{Z}^{[\frac{1+\pi}{n}]}$, where $\pi$ denotes the $q$-power Frobenius endomorphism, a point $P \in E(\mathbb{F}_q)$, where $\mathrm{End}_{\mathbb{F}_q}(E) = \mathbb{Z}^{[\frac{1+\pi}{n}]}$
**Output:** $\alpha(P)$
1: Write $\alpha = \frac{a+b\pi}{n}$, with $n$ of minimal denominator.
2: Compute a basis $\langle T_1, T_2 \rangle = E[n]$.
3: Compute the pairings $\mu_0 = e_{W,n}(T_1, T_2)$, $\mu_1 = e_{T,n}(T_1, P)$, $\mu_2 = e_{T,n}(T_2, -P)$.
4: Compute $r, s$ so that $\mu_r^* = \mu_2$ and $\mu_0^* = \mu_1$.
5: Set $T = [r]T_1 + [s]T_2$.
6: **return** $[b]T + [\frac{a+b}{n}]P$.

**Proposition D.1.** *Algorithm 4 is correct, and returns $\alpha(P)$ without using division points.*
*Proof.* Since $\alpha - \frac{a+b}{n} = b\frac{\pi-1}{n}$, and $\frac{a+b}{n} \in \mathbb{Z}$, we see that evaluating $\alpha(P)$ reduces to some simple scalar multiplications and additions after evaluating $\frac{\pi-1}{n}(P)$.

Notice that $E[n] \subseteq E(\mathbb{F}_q)$, since $\pi - 1$ factors through multiplication by $n$. Thus, we can sample a basis $\langle T_1, T_2 \rangle = E[n]$ over $\mathbb{F}_q$. Let $\frac{\pi-1}{n}(P) = T$. We have that $[n]T = (\pi - 1)(T) = 0_E$, since $T \in E(\mathbb{F}_q) = \ker(\pi - 1)$, thus $T$ can be written as $T = [r]T_1 + [s]T_2$ for some $r, s$. By the non-degeneracy of the Weil-pairing, we can recover $r, s$ by computing the evaluations $e_{W,n}(T_1, T_2), e_{W,n}(T_1, T), e_{W,n}(T_2, T)$. Finally, by [57, Equation 15], we have that

$$e_{W,n}(T_i, T) = e_{T,n}(T_i, P),$$

which finishes the proof. $\square$

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

In the CSIDH/CSURF supersingular case, the maximal denominator $n$ that can occur is $n = 2$, which is the situation we are in in Section 5.2. Further, we only ever need to evaluate the endomorphisms in multiples of our chosen basis points $P, Q$, where $P, Q$ are points of order $2^{f-1}$ for $2^f \mid p + 1$

Notice that in this case, the algorithm becomes particularly easy; assuming we aim to compute $\alpha(R)$ for $\alpha = \frac{a+\pi b}{n}$, we are in one of two situations:

- If $R = [m]P$ or $R = [m]Q$ for some $m \equiv 0 \pmod{2}$, we see that the pairings $\mu_1 = \mu_2 = 1$, and hence $\alpha(R) = [\frac{a+b}{2}]R$.

- If $R = [m]P$ or $R = [m]Q$ for some $m \equiv 1 \pmod{2}$, $\frac{\pi-1}{2}(R)$ is instead given by the unique 2-torsion point $T_i$ with $e_{T,2}(T_i, R) = 1$. But we already know that $e_{T,2}(T_{-1}, P) = e_{T,2}(T_1, Q) = -1$ by Lemma D.2, hence we see that a single Legendre test is enough to determine the correct $T_i$.

*Example D.1.* We of course have $\frac{\pi-1}{2}(T_1) = 0$, and by Lemma D.2, $\frac{\pi-1}{2}(T_0) = \frac{\pi-1}{2}(T_{-1}) = T_{-1}$ on $E$, and $\frac{\pi-1}{2}(T_0) = \frac{\pi-1}{2}(T_1) = T_1$ on $E^t$. If we sample $P$ above $T_1$, then if $e_{T,2}(T_0, P) = 1$, $\frac{\pi-1}{2}(P) = T_0$, whereas if $e_{T,2}(T_0, P) = -1$, $\frac{\pi-1}{2}(P) = T_1$.

Further, since translating by $T \in E[2]$ is still well defined on the Kummer line $E/\pm 1$, we see that evaluating our endomorphisms essentially only consists of a single scalar multiplication, and potentially a translation by $T_i$, which can all be done on the Kummer line.

*Remark D.2.* We can further refine our basis sampling to always output points such that $e_{T,2}(T_0, P) = e_{T,2}(T_0, Q) = 1$. By the discussion above, in this case we will always know that the potential translation by $T_i$ is by $T_i = T_0$. Sampling such a basis is easy to do using Lemma D.2: we want that $x_P - x(T_0)$ is a square, i.e. $x_P = x(T_0) + u^2$, and since $Q$ is on the twist, that $x_Q - x(T_0)$ is not a square, i.e. $x_Q = x(T_0) - u^2$. Then by the non degeneracy of the Tate pairing, we have $e_{T,2}(T_{-1}, P) = -1$ and $e_{T,2}(T_1, Q) = -1$, so by bilinearity we also have $e_{T,2}(T_1, P) = -1$ and $e_{T,2}(T_{-1}, Q) = -1$.
### D.3 Using the bottom of the volcano

We describe an alternative strategy to sample a basis on $E$ by using the descending isogeny $\phi_{T_0}: E \to E_0$. Recall that $T_0$ is the unique 2-torsion point on $E$ with non trivial self Tate pairing.

On $E_0$, we are at the bottom of the 2-isogeny volcano, so the 2-Sylow is cyclic, generated by a rational point $P_0$ of order $2^f$. The image of $P_0$ by $\hat{\phi}_{T_0}$ is then a point $P$ of order $2^{f-1}$ on $E$. By construction of $P$, we have $e_{T,2}(T_0, P) = 1$, while $e_{T,2}(T_{-1}, P) = e_{T,2}(T_1, P) = -1$ because $P$ is of full order on $E$.

Doing the same construction on $E_0^t$, we get a point $Q = \hat{\phi}_{T_0}(Q_0) \in E^t$, and $(P, Q)$ is an appropriate basis (as in Remark D.2) to which to apply the evaluation algorithm of Appendix D.2.

We saw on Remark D.1 that there are 4 possible choices of Montgomery coefficient for $E$ (or 2 if we restrict to $b=1$). Also, because the 2-torsion is rational on $E$, then $A+2$ being a square is equivalent to $A-2$ being a square, in which case $(1:1)$ and $(-1:1)$ lift to a rational point of 4-torsion on $E: y^2 = x^3+Ax^2+x$. In the other case, then $-A+2$ and $-A-2$ are squares, and $(1:1), (-1:1)$ lift to $E^t: y^2 = x^3-Ax^2+x$.

By contrast, on $E_0$, there is only one rational point of 2-torsion $T'_0$ (which generates the kernel of $\hat{\phi}_{T_0}$), which has trivial self Tate pairing, so we can find a Montgomery form for $E_0$ with $T'_0$ sent to $(0,0)$. Then $E_0$ is isomorphic to $by^2 = x^3 + A_0x^2 + x$, and replacing $A_0$ with $-A_0$ if necessary, we can assume $b=1$. In particular, there is only one choice for $A_0$ to put $E_0$ in Montgomery form $E_0: y^2 = x^3 + A_0x^2+x$. Furthermore, either $A_0+2$ is a square (so $(1,1)$ lifts to $E_0$) and then $A_0-2$ is not a square (so $(-1,1)$ does not lift), so $-A_0+2$ is a square and $-A_0-2$ is not a square. Conversely, $A_0+2$ not being a square is equivalent to $A_0-2$ being a square, $-A_0+2$ not a square, $-A_0-2$ a square.

We see that $A_0$ gives a convenient way to represent not only $E_0$, but also $E$: we have $T'_0 = (0,0)$, $\phi_{T'_0}$ the canonical ascending isogeny $E_0 \to E$, whose dual $\phi_{T_0}$ has kernel generated by $T_0$, and also that if $A_0+2$ is a square then the image of $(1:1)$ by $\phi_{T'_0}$ gives $T_1$, and the image of $(-1:1)$ gives $T_{-1}$ (and conversely if $A_0+2$ is not a square).

A last advantage is that we can efficiently sample $P_0 \in E_0(\mathbb{F}_p), Q_0 \in E^t(\mathbb{F}_p)$ of full order $2^f$ by adapting the entangled basis algorithm of [66]. We first sample $x(P_0)$ a non square in $\mathbb{F}_p$ of the form $-A_0/(1-u^2)$ until we find that $x(P_0)$ lifts to a rational point $P_0$ on $E_0$. Since the Tate pairing of $T'_0$ with $P_0$ is non trivial by construction, $P_0$, multiplied by the cofactor, has exact order $2^f$. Now set $x(Q_0) = -A_0 - x(P_0) = -u^2x(P_0)$. Then $x(Q_0)$ is a square, and $x(Q_0)^2 + A_0x(Q_0) + x(Q_0) = x(P_0)^2 + A_0x(P_0) + x(P_0) = y(P_0)^2/x(P_0)$ is not a square. So $x(Q_0)$ lifts to a point on the twist $-y^2 = x^3+A_0x^2+x$, with $y(Q_0) = uy(P_0)$. And the non reduced Tate pairing of $T'_0$ and $Q_0$ is given by the class of $-x(Q_0)$, which is not trivial. Hence $Q_0$, multiplied by the cofactor, also has exact order $2^f$ on $E^t$. We can precompute a table of elements such that $-1/(1-u^2)$ is not a square and another table of elements such that $-1/(1-u^2)$ is a square. Then for the entangled basis generation, depending on whether $A_0$ is a square or not, we use the first or the second table.
Remark D.3 (The horizontal isogenies of degree $2^{f-1}$). On $E$, which is at the crater, the ideal (2) splits as $(2) = \mathfrak{p}_1\mathfrak{p}_{-1}$. We have $2^f \| \pi^2 - 1$, and $\mathfrak{p}_1^{f-1}\mathfrak{p}_{-1} \mid \pi - 1$, $\mathfrak{p}_{-1}^{f-1}\mathfrak{p}_1 \mid \pi + 1$.

Let $(P_0, Q_0)$ be our basis sampled on $E_0$ as above, $P = \hat{\phi}_{T_0}(P_0)$, $Q = \hat{\phi}_{T_0}(Q_0)$ the points of order $2^{f-1}$ on $E$. Then since $T'_0$ has non trivial pairing with $P_0$, then by [57, Example 5.16], the isogeny generated by $P_0$ lands on the bottom of the volcano, hence $\langle P \rangle$ is not the kernel of the isogeny associated to $\mathfrak{p}_1^{f-1}$, as can also be seen from the fact that it has non trivial pairing with $T_1$. However, if $P' = P + T_{-1}$ or $P' = P + T_0$, then $e_{T,2}(T_1, P') = 1$ by Lemma D.1, and so $\langle P' \rangle$ is the kernel of $\mathfrak{p}_1^{f-1}$. A similar argument shows that $E[\mathfrak{p}_{-1}^{f-1}] = \langle Q + T_0 \rangle$.

Remark D.4 (The 2-isogenies between $E$ and $E_0$). As illustrated above, it is convenient to use the 2-isogenies $\phi_{T_0}$ and $\phi_{T'_0}$ to move between $E$ and $E_0$. We also would like to have $E$ in Montgomery form, and $E_0$ too for the entangled basis generation. To put a curve $E$ in Montgomery form, we need the $x$-coordinate of a point of 4-torsion $T' \in E$ such that $x(T') \in \mathbb{F}_p$. In particular, if we start with $E_0$ in Montgomery form, and take the quotient by $T'_0$, we easily find $E$ in Legendre form (see for instance [59, Example B.1]), but finding a Montgomery form for $E$ requires either a square root or to push the 8-torsion point $2^{f-3}P_0$ (or $2^{f-3}Q_0$) to $E$ via $\phi_{T'_0}$.

On the other hand, if we start from $E$ in Montgomery form and we are given $\alpha = x(T_0)$, then we can directly find $E_0$ in Montgomery form using the formulas from [52, Proposition 2]. Hence it may be useful to use $\alpha$ to represent $E$ (or even use the level 2 theta null point), as already suggested in Remark D.1, to save a square root.

## E Parameter choices

In this section we give more details about the choice of the set of primes to include in $\mathfrak{B}$. Clearly, adding more primes to $\mathfrak{B}$ improves the time of Step 1 by making the values of $N_1$ and $N_2$ in Equation (11) smaller. On the other hand, we are on average required to compute more and isogenies of larger degree in Step 2. If we start from a very small set of primes, say $\mathfrak{B} = \{2, 3\}$, and keep adding the smallest possible prime to this set, we hence expect the time of Step 1 to decrease, up to a point where most expensive operations are eventually checking the validity of a $u, v$ pair instead of finding one. Since on average we need to check a constant number of pairs $u, v$, at some point adding primes to $\mathfrak{B}$ will not improve Step 1, while still slowing down Step 2. This is exactly what we observe in Figure 2, where we compare the times for different choices of $\mathfrak{B}$ in the 500 parameter set. In each different column, all the (Elkies) primes up to the last prime reported are included in $\mathfrak{B}$. It is then clear that the best choice is $\mathfrak{B} = \{2, 3, 7, 11, 13\}$. Adding 29 or bigger primes has essentially no impact on Step 1, while negatively affecting Step 2. Similar experiments were performed for all security levels, and the results in Table 1 were derived accordingly.
<figure><img src="image_11.png" /><figcaption>**Fig. 2.** Time in seconds for Steps 1 and 2 for different choices of $\mathfrak{B}$ in the 500 parameter set.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 3.** Success rates of Step 1 for different choices of $\mathfrak{B}$.</figcaption></figure>

A more detailed analysis is shown in Figure 3. Looking again at the 500 parameter set, on the *x* axis we have the last prime included in $\mathfrak{B}$, but this time including also non-Elkies primes. So for instance the column 7 corresponds
to $\mathfrak{B} = \{2,3,5,7\}$, while in Figure 2 it corresponds to $\mathfrak{B} = \{2,3,7\}$, 5 being non-Elkies. Every dot represents a random ideal passed to Step 1; cyan dots correspond to successful runs, while gray dots correspond to failures. For each dot, the $y$ axis indicates the time taken to run Step 1. The green and red dots show the average time for the given set, where green means that the denoted prime is Elkies, and red that it is non-Elkies. Here we see that adding non-Elkies primes to $\mathfrak{B}$ does not help significantly, since they cannot appear as factors in the ideals. Moreover, as already discussed, they would significantly slow down Step 2 forcing us to work over extension fields.