Finding Orientations of Supersingular Elliptic Curves and Quaternion Orders

Orientations of supersingular elliptic curves encode the information of an endomorphism of the curve. Computing the full endomorphism ring is a known hard problem, so one might consider how hard it is to find one such orientation. We prove that access to an oracle which tells if an elliptic curve is $\mathfrak{O}$-orientable for a fixed imaginary quadratic order $\mathfrak{O}$ provides non-trivial information towards computing an endomorphism corresponding to the $\mathfrak{O}$-orientation. We provide explicit algorithms and in-depth complexity analysis. We also consider the question in terms of quaternion algebras. We provide algorithms which compute an embedding of a fixed imaginary quadratic order into a maximal order of the quaternion algebra ramified at $p$ and $\infty$. We provide code implementations in Sagemath which is efficient for finding embeddings of imaginary quadratic orders of discriminants up to $O(p)$, even for cryptographically sized $p$.


Introduction
Isogeny-based cryptography is a relatively new branch of post-quantum cryptography which is based on hard problems supposedly intractable even for quantum computers.The underlying hard problems were first introduced publicly in 2006 by the hash-function proposal of Charles-Goren-Lauter [11], and the works of Couveignes [14] and Rostovtsev-Stolbunov [48].Since then, this field has blossomed with the introductions of new schemes such as SIDH [29] (now broken by [9,37,47]), CSIDH [10], and SQISign [19].The hardness of all isogeny-based schemes is based on some variant of the path finding problem, which asks to find an isogeny between two given supersingular elliptic curves.The quaternion analogue of this hard problem has been efficiently solved [32], but the problem remains hard for supersingular elliptic curves.Path finding in the supersingular isogeny graph is equivalent to endomorphism ring computation, which was first heuristically proven in [22] and then rigorously (assuming GRH) proven in [58].The key recovery of CSIDH was reduced to endomorphism ring computations in [57].
To study the hardness of the path finding problem it is natural to add some data to the elliptic curves and study how this data interacts with the graph structure.One way to do this is to add the information of an orientation to the elliptic curve vertices.Informally, an orientation on an elliptic curve E is an embedding of an imaginary quadratic order O into the endomorphism ring of E which cannot be extended to a superorder of O.The resulting isogeny graph admits an abelian group action, which is used in cryptographic protocols such as CSIDH [10], Scallop [26], OSIDH [12], and SETA [18].The group action is crucial for defining the Uber isogeny problem, whose hardness underlies all isogeny-based schemes.One might suspect that being given the information of an orientation could weaken the difficulty of the path finding problem, but this depends heavily on the given orientation and does not typically weaken the hardness of the path finding problem [1,57].A natural question to consider would be how to find an orientation on a curve, given that one exists.This is the O-Orienting Problem.It is also natural to consider the decisional version of this problem: given a supersingular elliptic curve E and a quadratic order O, can one decide whether E is orientable by O? Solving either of these problems would completely break OSIDH [12,15,40].Interestingly, they are not even efficiently solved on the quaternion side.In this work, we give reductions between the search and decision variants of these problems, and provide algorithms for the quaternion variant of these problems.
1.1.1.Reduction from Search to Decision O-Orienting Problem.When the discriminant of O is smaller than the characteristic p of the base field, we prove a subexponential reduction from the computational to the decisional version of the O-Orienting Problem.In particular, we provide an explicit algorithm (Algorithm 4.3) to find an O-orientation of an orientable elliptic curve in subexponential time and space when given access to an oracle deciding whether any elliptic curve is O-orientable.In Section 4.2 we provide an in-depth analysis and proof of the complexity of the algorithm (Theorem 4.13).This proves that such an oracle gives non-trivial information since finding an orientation automatically yields a non-scalar endomorphism and the best known algorithms to find a non-scalar endomorphism on a supersingular elliptic curve are exponential [20,23].
Before treating the general case, we prove a polynomial reduction when O is the maximal order of Q( √ −d) and d is the product of small distinct primes in Section 3.This allows us to illustrate the spirit of the more general algorithm in a less complicated setting.We provide an explicit algorithm for this case (Algorithm 3.1) and prove in Theorem 3.10 that this algorithm runs in polynomial time.
1.1.2.Quaternion Order Embedding Problem.In Section 5, we consider the Quaternion Order Embedding Problem (Problem 5.1) which is the quaternion analogue of the O-Orienting Problem.That is, given a maximal quaternion order O ⊂ B p,∞ and a quadratic order O which embeds into O, find an embedding ι : O ֒→ O that cannot be extended to a superorder of O.In Section 5.1 we present a general algorithm to solve the problem of finding embeddings using a factorization oracle.We provide a complexity analysis based on several heuristics in Section 5.2.In Section 5.5 we show that finding embeddings which cannot be extended (i.e., orientations), only adds a small factor to the running time.We prove efficiency for the curve with j-invariant 1728, and describe a practical method of removing the dependence on the factorization oracle.When the discriminant disc(O) is small, our algorithm improves the state of the art being efficient up to disc(O) = O(p).We provide an implementation in Sagemath [53] which, for small discriminant orders, is fast for cryptographically sized p.
Code is available at: https://github.com/jtcc2/finding-orientations Acknowledgements.This project began at KU Leuven Isogeny Days in 2022, and the authors extend their gratitude to the organizers.

Preliminaries
We provide a concise summary of the necessary background and the state of the art algorithms which we use in this paper.
2.1.Supersingular elliptic curves and quaternion algebras.Let p be a prime.An elliptic curve over F p is called supersingular if any one of the following equivalent conditions holds : (1) End(E) is isomorphic to a maximal order in a quaternion algebra (2) E[p r ] = 0 E for all r ≥ 1 (3) j(E) ∈ F p 2 and the multiplication-by-p map [p] is purely inseparable (4) The dual to the p r -power Frobenius is purely inseparable for all r ≥ 1.
See [52] for additional properties and proofs of equivalence.We use the endomorphism ring heavily in what follows, so we describe here the necessary definitions and properties of quaternion objects.For more generality and more detail, we encourage the reader to see [54].
A (definite) quaternion algebra A is a noncommutative algebra which has rank 4 over Q, and can be specified by generators i, j such that: An order O in A is a Z-submodule of A of rank 4 which is also a subring.An order is said to be maximal if it is not properly contained in any other order.For any lattice I in A, we define its left order O L (I) := {α ∈ A : αI ⊆ I}.
The right order O R (I) is defined analogously.A lattice I of A is said to be invertible if there exists a lattice I ′ such that For every order O of A we can define a left class set of equivalence classes of invertible ideals: Two invertible left-O ideals are equivalent in the left class set of O if they differ by a unit of A. The left class set of invertible ideals is finite.The right class set of invertible ideals is analogously defined and is also finite.
For a fixed prime p, we define the (unique up to isomorphism) quaternion algebra B p,∞ to be the definite quaternion algebra ramified precisely at p and ∞.The endomorphism rings of supersingular elliptic curves over F p are isomorphic to maximal orders in B p,∞ : Theorem 2.1 (Deuring [21]).Fix a maximal order M of the quaternion algebra B p,∞ ramified precisely at p and ∞.There is a bijection between isomorphism classes of supersingular elliptic curves over F p and the left class set of the order M .
Given a supersingular elliptic curve E/F p , one might ask to compute End(E) in different forms: to compute endomorphisms of E which generate End(E), or to compute the isomorphism class of End(E) abstractly in the quaternion algebra B p,∞ .This problem is computationally difficult in all formulations.The information of one endomorphism ω of E reveals an imaginary quadratic order Z[ω] embedded within End(E).In Section 2.2, we provide more background information on such embeddings.

Orientations. Definition 2.2 (Orientation).
Let O be an imaginary quadratic order.An Oorientation of a supersingular elliptic curve E/F p is an embedding ι : O ֒→ End(E) which cannot be extended to a larger order containing O. The pair (E, ι) is called an O-oriented supersingular elliptic curve.Definition 2.2 corresponds to the definition of primitive O-orientation found elsewhere in the literature [1,12,40].We omit the word "primitive" in our definition, as almost all of our O-orientations are primitive.When we want to discuss an embedding O ֒→ End(E) which can be extended to a superorder of O, we highlight this by using the term "imprimitive".
The notion of an orientation as in Definition 2.2 was recently introduced to isogeny-based cryptography by Colò and Kohel [12] and was subsequently studied [1,2,26,40,57].The quaternion counterpart of this notion has a longer history, dating back to Chevalley, Hasse, and Noether and often referred to as the theory of optimal embeddings.Supersingular elliptic curves which admit an O-orientation are called O-orientable.There is an action of the class group Cl(O) on the set of O-oriented supersingular elliptic curves induced by the following action of an invertible O-ideal a: where E a is the codomain of the degree-N (a) isogeny For an imaginary quadratic field K, deciding if K embeds into the quaternion algebra B p,∞ is the simple matter of the splitting behavior of p in K.However, for a particular imaginary quadratic order O and a particular supersingular elliptic curve E, it is generally difficult to decide if E is O-orientable.Naturally we are inclined to study the following problems and the relationship between them: One may also consider the group action variant of the Uber-isogeny problem, originally introduced in [18], although we do not pursue this perspective in this work: Problem 2.6 (O-Uber Isogeny Problem).Given a supersingular elliptic curve E with an O-orientation ι : O ֒→ End(E) and an O-orientable supersingular elliptic curve F , find an ideal a ∈ Cl(O) such that a * E = F .

2.4.
Computing an ℓ-isogeny between two j-invariants.Given two supersingular j-invariants j(E) ∈ F p 2 and j(E ′ ) ∈ F p 2 we explain how to find an ℓ-isogeny φ : By [6,Theorem 2], given Weierstrass equations of E and E ′ , we can find (if it exists) a normalized ℓ-isogeny φ : E −→ E ′ with only Õ(ℓ) arithmetic operations over F p 2 .By normalized, we mean that φ pulls back the invariant differential The existence of such a normalized isogeny φ only depends on the choice of Weierstrass equations for E and E ′ which determine the constant λ := φ * ω ′ /ω.Knowing only j(E) and j(E ′ ), we have multiple choices of Weierstrass equations and we have to pick one so that λ = 1.We fix an equation for E : y 2 = x 3 +Ax+B, then find an equation for E ′ so that λ = 1.Following the method given by [50, Section 7] (referring to ideas introduced in [24, Section 3]), we take E ′ : A , where The derivatives ∂Φ ℓ /∂X and ∂Φ ℓ /∂Y can be precomputed with Õ(ℓ 2 log(p)) operations over F p 2 using the techniques in Section 2.3 (see [43,Remark 5.3.10]).Hence, in total, computing an isogeny φ : E −→ E ′ costs Õ(ℓ 2 log(p)) operations over F p 2 when j(E ′ ) = 0, 1728 and ∂Φ ℓ /∂Y (j(E), j(E ′ )) = 0.The cases j(E ′ ) = 0, 1728 are very unlikely (probability O(1/p) in the supersingular isogeny graph).The latter case we split into two: First, suppose ∂Φ ℓ /∂Y (j(E), j(E ′ )) = 0 and ∂Φ ℓ /∂X(j(E), j(E ′ )) = 0.By symmetry of Φ ℓ (X, Y ), we can fix a Weierstrass equation for E ′ and find a normalized Weierstrass equation for E by Equations 1, 2 and 3 (after swapping E and E ′ ).The remaining case is we prove in Appendix A that this is very unlikely when log(ℓ) ≪ log(p) (which will be the case in our paper).
We can still handle singular cases at a higher cost of Õ(ℓ 7/2 ) operations over F p 2 with a naive algorithm.We enumerate all the cyclic subgroups of order ℓ of E[ℓ] (there are ℓ+1 of them) and use [5] to compute each ℓ-isogeny with O( √ ℓ) operations over the field extension K/F p 2 where E[ℓ] is defined.As will be proved in Lemma 2.11, K has degree O(ℓ) over F p 2 so one arithmetic operation over K is equivalent to at most O(ℓ 2 ) operations over F p 2 .Since singular cases are very unlikely when log(ℓ) ≪ log(p), we may assume throughout this paper that computing ℓ-isogenies between j-invariants costs Õ(ℓ 2 log(p)) operations over F p 2 on average by Lemma A.2.

2.5.
Efficiently representing an isogeny of any degree with Kani's lemma.Let ϕ : E −→ E ′ be an isogeny of degree d between supersingular elliptic curves E, E ′ /F p 2 .In general, we can represent ϕ with data of size O(d).We either have direct formulas to evaluate ϕ (given by rational fractions) or equivalently, generators of the kernel (defined over an F p 2 -extension of degree O(d)) from which we can derive these formulas by [56].In this case, evaluating ϕ on a point takes linear time in d.We can do much better when d is smooth by representing ϕ as a product of small degree isogenies.This is an efficient representation, in the sense of the following definition.Definition 2.7.[16, Definition 1.1.1]An efficient representation of an isogeny ϕ : E −→ E ′ defined over a finite field F q is given by a couple (D, A ) where: (i): D is some data of size polynomial in log(deg(ϕ)) and log(q) determining the isogeny ϕ in a unique way.(ii): A is a universal algorithm independent of ϕ returning ϕ(P ) as input D and P ∈ E(F q k ) in polynomial time in k log(q) and log(deg(ϕ)).
We can also efficiently represent ϕ when d is not smooth.Provided we can evaluate ϕ on some torsion points, we know that we can "embed" ϕ in a smooth degree higher dimensional isogeny F .Knowing F , we can evaluate ϕ everywhere in polynomial time.This provides an efficient representation of ϕ.This idea was first introduced in the attacks against SIDH [9,37,44] and then reused for several other applications [16,45,46].Definition 2.8 (d-isogeny in higher dimension).Let α : (A, λ A ) −→ (B, λ B ) be an isogeny between principally polarized abelian varieties (PPAV).We denote by α the isogeny where α is the dual isogeny of α.
We say that α is a d-isogeny We use the following result due to Kani [31,Theorem 2.3].A concise expression of this result may be found in [44,Lemma 3.6].Lemma 2.9 (Kani).Consider a commutative diagram of isogenies between PPAV: where ϕ and ϕ ′ are a-isogenies and ψ and ψ ′ are b-isogenies.
Then, the isogeny F : A × B ′ −→ B × A ′ given in matrix notation by is a d-isogeny with d := a + b, for the product polarizations.If a and b are coprime, the kernel of F is Let N > d be a powersmooth integer coprime with d.We can always write N = d + a 2 1 + a 2 2 + a 2 3 + a 2 4 for some a 1 , a 2 , a 3 , a 4 ∈ Z, by Legendre's four square theorem.Let α ∈ End(E 4 ) be the isogeny written in matrix form as follows: and α ′ be its analogue in End(E ′ ).Let Φ := Diag(ϕ, ϕ, ϕ, ϕ) : E 4 −→ E ′4 .Then, Φ is a d-isogeny, α and α ′ are (N − d)-isogenies and we have a commutative diagram: that yields an 8-dimensional N -isogeny: with kernel: (5) ker(F ) := {( α(P ), Φ(P )) since N and d are coprime.By the above formula, we can compute ker(F ) if we can evaluate ϕ on generators of E[N ].We can then compute F as a product of small degree isogenies (as in dimension 1) and evaluate ϕ efficiently everywhere as a component of F .
Lemma 2.10.Let the setup be that of the previous paragraph, and suppose N = s i=1 q ei i , where q 1 , • • • , q s are distinct primes.Then, we can decompose F as: and where F i is a q ei i -isogeny for all i ∈ {1, • • • , s}.Let K := ker(F ).Moreover, ] and the primes q j = q i for j > i, so we must have Each of the F i in Lemma 2.10 can be decomposed into a chain of q i -isogenies.As suggested in [16, § 5.5], those isogeny chains can be computed in quasi-linear time using the optimal strategy introduced for SIDH [29, § 4.2.2].Each q i -isogeny of the chain can be computed in time O(q 8 i ) with the theta model [36].We summarize this computation in Algorithm 2.1.
Algorithm 2.1: EfficientRep returning an efficient representation of an isogeny.
Data: An isogeny ϕ : E −→ E ′ between supersingular elliptic curves and a smoothness bound D > 0. Result: An 8-dimensional isogeny F of D-powersmooth degree representing F .For all j ∈ {1, 2} and k ∈ {1, 2, 3, 4}, let P i,j,k ∈ E 4 [q ei i ] be the tuple with P i,j in position k and 0 elsewhere; Lemma 2.11.Let E/F p 2 be a supersingular elliptic curve and n ∈ Z >0 .Then E[n] is defined over an extension of degree at most 6φ(n) of F p 2 , where φ is Euler's totient function.
Hence, E[n] is defined over F p 2δ provided that 12 and φ(n) divide δ.If n has an odd prime factor or n = 2 k with k ≥ 2, then φ(n) is even so 12 | 6φ(n) so δ = 6φ(n) satisfy the desired conditions.If n = 2, then E[n] is formed by 0 and the (x, 0) where x is the root of a cubic polynomial equation over F p 2 (Weierstrass equation of E).Hence, E [2] is defined over an extension of degree at most 3 of F p 2 .Proposition 2.12.Finding the a i on line 3 costs O(log(N )/ log log(N )) = O(log(d)) with Pollack and Treviño's algorithm [42].
To generate such a basis, we first sample a random point P ′ ∈ E(F p 2δ ) and compute P i,1 := [M/q ei i ]P ′ , where M := #E(F p 2δ ) until P has order q ei i .By the same method, we sample elementary arithmetic operations by Cipolla-Lehmer's algorithm.E being supersingular, point counting on E to compute M costs O(1) and the scalar multiplication by M/q ei i costs O(log(M )) = O(D log(p)) arithmetic operations over F p 2δ (costing O(D 2 M (p)) each).Testing that (P i,1 , P i,2 ) is a basis costs O(D) elliptic curve additions so O(D) arithmetic operations over F p 2δ (we compute the [k]P i,1 and [l]P i,2 for 1 ≤ k, l ≤ q ei i −1 and conclude that we have a basis if these two sets are disjoint).Only O(1) samplings of P i,1 and P i,2 are necessary before we find a basis.Hence, the overall complexity of the basis computation is Line 7 costs two evaluations of ϕ and line 9 costs eight scalar multiplications by the a i , costing O(log(d)) operations over F p 2δ each.Hence, the total cost of the loop of lines 5-10 is elementary arithmetic operations and 2s = O(log(d)) evaluations of ϕ.
Finally, computing each ) arithmetic operations over F p 2 and computing the basis ) arithmetic operations over an extension of degree O(D) of F p 2 .The total cost of the loop of lines 12 -17 is 2.6.Smoothness test and factorization with the ECM method.In this section, we explain how to test if an integer N is B-smooth and find its factorization if it is the case.A naive method would be to use trial division, but it is not optimal when B is subexponential (which will be the case in the paper).An alternate method would be to factor N with the General Number Field Sieve (GNFS) [8] and test if its prime factors are ≤ B. However, GNFS underperforms with smooth integers so we propose a faster method relying on the elliptic-curve factorization method (ECM) due to Lenstra [34].Finding a prime factor of N with this method takes time L ℓ (1/2, √ 2), where ℓ is the smallest prime divisor of N and with the usual notation where o(1) is for x → ∞.Hence, to test the B-smoothness of N , we simply apply ECM to find a factor k | N after expected time L B (1/2, √ 2).If the running time exceeds what it should be, it means that N is not B-smooth and we stop.Otherwise, we continue and try to factor k and N/k recursively until we have either completely factored N or concluded it is not B-smooth.Algorithm 2.2 follows.
If r is the number of prime divisors of N (θ) (with multiplicity), then r = O(log(|∆ O |)) and Algorithm 2.2 can terminate with at most r calls to ECM so in time rL

Reduction of O-orienting problem for special discriminants
We begin with a nice assumption about the discriminant of our imaginary quadratic order O.The key ideas from this special case provide a foundation for the general cases we consider in Section 4.
Algorithm 2.2: SmoothFact determining if an integer is smooth and returning its factorization.
Data: An integer N ∈ Z >0 and a smoothness bound B > 0. Result: Return ⊥; We use an oracle which solves Problem 2.3 to find an endomorphism ϕ of E to which we map α, thus determining an embedding O ֒→ End(E) either by mapping α = ω to ϕ or (1 + α)/2 = ω to (1 + ϕ)/2.We use the fact that the primes ℓ i are ramified in K.
We walk the O-oriented ℓ i -isogeny volcanoes in order to obtain the endomorphism ϕ on E which is the image of the generator ω under an embedding ι : O ֒→ End(E).The ideals l i above ℓ i in O determine horizontal degree-ℓ i isogenies between O-oriented curves, beginning and ending with E. To see this, we need the following fact about horizontal isogenies of O-oriented elliptic curves: , where ϕ i is an isogeny of degree ℓ i for all i ∈ {1, • • • , r}.For each i, let Oℓ i = (l i ) 2 .The ideals l i determine the horizontal isogenies of O-oriented curves: Lemma 3.2.In the setting described above, all of the isogenies ϕ i in the decomposition of ϕ are horizontal.
Proof.Since N (α) = r i=1 ℓ i , we have Oα = r i=1 l i , l i being the unique prime ideal of O lying above ℓ i for all i ∈ {1, • • • , r}.Hence, the ϕ i intervening in the decomposition of ϕ = ι(α) are horizontal isogenies given by the action of l i .

Now, we describe the steps to obtain an endomorphism
which will be the image of α.Let E 0 := E.
For i = 0, we find the unique isogeny ϕ 1 : E 0 −→ E 1 which corresponds to the action of [l 1 ] on (E 0 , ι) by computing each of the ℓ 1 + 1 outgoing isogenies and querying our oracle to find the one whose codomain E 1 is in fact orientable by O.We continue this process to compute each ϕ i by using the oracle to find the correct degree-ℓ i isogeny to another O-orientable curve.At the last step, we compute the degree-ℓ r isogeny from E r−1 −→ E r and then post-compose with an isomorphism E r ∼ = E 0 : We let ϕ r denote this composition.See Algorithm 3.1 for the algorithmic description of this process.The prime ideal 2 splits in K so there are two horizontal 2-isogenies and one descending 2-isogeny with domain E 0 .However, all these isogenies have the same codomain E 1 (up to isomorphism) with j-invariant j(E 1 ) = 3.So E 1 is both O Koriented and (Z + 2O K )-oriented.

Motivated by this example a question arises: If ψ
i O-orientable, how do we know that ψ i is the unique horizontal isogeny ϕ i given by the action of l i on (E i−1 , ι)?In fact, ϕ i and ψ i could be distinct horizontal isogenies for distinct primitive orientations

3). Or ψ i could even be descending and E ′
i O-oriented by a different orientation than the one induced by ψ i , (E i , (ψ i ) * (ι)).To exclude those cases, we assume p > |∆ O | max 1≤i≤r ℓ i and prove that there is precisely one (primitive) Oorientation ι on E i−1 , which ensures that there is only one isogeny ϕ i corresponding to the action of [l i ] on (E i−1 , ι).We also prove that codomains of descending isogenies are not O-orientable.These are consequences of [30, Theorem 2'], that we recall below.
so it is either the identity or the complex conjugation.The result follows.
O < p 2 by hypothesis.Contradiction.Remark 3.6.Corollary 3.5 holds for any imaginary quadratic order O, not only the special form we consider in this section.
given by the action of l 1 is uniquely determined, and it is the only ℓ 1 -isogeny with O-oriented codomain.In this case, ϕ 1 can be distinguished from other ℓ 1 -isogenies by an oracle query.Similarly for each further i ∈ {2, ..., r}, the isogeny ϕ i : ) by computing each of the ℓ i + 1 isogenies and querying the oracle to find the one whose codomain E i is orientable by O.In particular, the isogeny ϕ r : Possibly post-composing with this isomorphism, we have an endomorphism ϕ = ϕ r • • • • • ϕ 1 ∈ End(E) associated to the action of the ideal r i=1 l i = αO.It follows that ϕ = τ • ι(α) for some automorphism τ ∈ Aut(E).We may post-compose ϕ by τ ∈ Aut(E) until the result has trace zero, as α.The trace can be computed in polynomially many isogeny evaluations using Schoof's algorithm [50, Section 5].
Remark 3.7 (Isomorphisms).Assuming we are working with elliptic curves in Weierstrass form, all isomorphism formulae are known.To find an isomorphism β : E r −→ E 0 , we check the codomain formula for each isomorphism from E r until E 0 is found.
There additional automorphisms in the two special cases of j = 1728 and j = 0 [2, Figure 3.1, Section 6].At each step ϕ i : E i−1 → E i where j(E i ) = 0 or 1728, we must decide whether or not to post-compose with these automorphisms.The automorphisms [±1] will not affect the resulting trace, but we must check one nontrivial automorphism for j = 1728 and two for j = 0.This can be done after the algorithm is completed, as the oracle calls will remain unaffected.
The additional running time of choosing isomorphisms can be bounded by a constant, so does not contribute to the overall complexity.Proof.We justified above that this algorithm terminates and is correct.For all i ∈ {1, ..., r}, this algorithm computes the ℓ i + 1 curves which are ℓ i -isogenous to E i−1 , which costs Õ(ℓ 2 i log(p)) operations over F p 2 by Section 2.3.It calls the oracle IsOrientable O ℓ i + 1 times and computes one ℓ i -isogeny between j(E i−1 ) and j(E i ), which costs on average Õ(ℓ 2 i log(p)) operations over F p 2 by Section 2. The idea is similar the case of special discriminant considered in Section 3. We compute an endomorphism corresponding to a generator of O as a chain of horizontal isogenies of small degrees.However, two difficulties arise.First, the canonical generator ω := (s + √ ∆ O )/2 with s := ∆ O mod 2 of O is not smooth in general.We have to find another smooth generator θ of O. Second, if we denote ϕ := ι(θ) and decompose ϕ := ϕ r • • • • • ϕ 1 as a product of horizontal isogenies of degrees ℓ i , • • • , ℓ r respectively, we may not be able to find the ϕ i simply by using the oracle IsOrientable O as in Section 3. We are no longer guaranteed that ℓ i | ∆ O , so there may be 1 + (∆ O /ℓ i ) = 2 horizontal isogenies of degree ℓ i from a O-oriented elliptic curve.To search for ϕ, starting at root E we fill a binary tree whose nodes are of an endomorphism To optimize the tree search, we propose a meet-in-the middle strategy where two half-depth such trees are computed starting at E instead of a single one: ( We look for θ of the form θ := a + ω with a ∈ Z to be determined.There is no better known method to find a than sampling a randomly and testing whether , the optimal known way to test the Bsmoothness of N (a + ω) is the method introduced in Section 2.6 using ECM with time complexity L B (1/2, √ 2).Algorithm 4.1 presenting the search for θ = a + ω follows.
We assume p > B|∆ O | so any O-orientable curve admits a unique O-orientation up to conjugation by Corollary 3.5(i).Hence, every node of T can be represented by j-invariant (the root and E i,2 are the only two O-orientable curves that are ℓ iisogenous to E i , given by the action of ideals l i , l i above ℓ i .As in Section 3, to find E i,1 and E i,2 we compute the codomain j-invariants of all degree-ℓ i isogenies E i −→ E ′ and apply the decision oracle to see which are O-orientable.Determining such j-invariants can be done using modular polynomials in Õ(ℓ 2 i log(p)) operations over F p 2 , as in Section 2.   Since N (θ) = deg(ϕ) is coprime with ∆ O , ℓ 1 does not divide the conductor of O so ϕ 1 is horizontal or descending by Proposition 3.1(i).Suppose ϕ 1 is descending.Since ϕ is horizontal, there are as many descending as ascending ℓ 1 -isogenies in the decomposition of ϕ.Reordering the ϕ i if necessary, we can assume ℓ 1 = ℓ 2 and ϕ 2 is ascending.Since there is only one ascending isogeny, ϕ 2 must be the dual of ϕ 1 , up to post composition with an isomorphism.Thus, ϕ 2 • ϕ 1 factors through [ℓ 1 ] and so does ϕ.Then ϕ/[ℓ 1 ] ∈ End(E) and ι is not a primitive orientation, contradicting our original assumption.
By the same argument, the ϕ i are also horizontal for i ≥ 2.
Since the ϕ i are horizontal, we may use O-oriented isogeny trees to find these isogenies.Let s := ⌊r/2⌋, T 1 the O-oriented (ℓ 1 , • • • , ℓ s )-isogeny tree starting at E 0 := E and T 2 the O-oriented (ℓ s+1 , • • • , ℓ r )-isogeny tree starting at E. Assume we have found a common leaf E s in T 1 and T 2 .The branch of T 1 of leaf E s is a chain of horizontal ℓ i -isogenies ψ i : Let (E 0 , ι) be an O-oriented supersingular elliptic curve and ψ ∈ End(E 0 ) a horizontal endomorphism of degree coprime to p.Then, there exists α ∈ O such that ψ = ι(α).Proof.
However, the factors ψ i of ψ have subexponential degree so they do not provide an efficient representation of ψ (enabling to evaluate ψ in polynomial time for instance).We apply EfficientRep Algorithm 2.1 to get an efficient representation of ι(ω) or ι(ω) = ±ψ − [a].The search to decision reduction Algorithm 4.3 follows.
For efficiency, only j-invariants are stored in the trees and not the ℓ i -isogenies relating them so we use the method of Section 2.4 to recover them in time Õ(ℓ 2 i log(p)).4.2.Complexity analysis.

Complexity of the smooth norm search (Algorithm 4.1).
To estimate the complexity of Algorithm 4.1, we need to determine the probability that N (a + ω) is ns-(B, r m , ∆ O )-smooth.We have proven results on the distribution of B-smooth integers among random integers but not for random values of quadratic integer polynomials.For that reason, we introduce the following heuristic assumption.
This heuristic assumption is supported by an estimate on B-smooth values of polynomials very similar to random integers.Such an estimate has been proved in [38,Theorem 1.1] under a dual hypothesis on the number of prime values of polynomials when B is in a very tight range.It has been conjectured [27, Equation 1.20] that this result holds for broader values of B. Lemma 4.8.Let Ψ r (x, y, d) denote the number of (y, r, d)-smooth integers ≤ x: Proof.The proof follows from [17, § 2].We have the following inequalities (following from set inclusions): To conclude, we compute the cardinality of for k, n ∈ Z >0 and apply it to the last set in the inequalities above.The set S(k, n) is in bijection with the subsets of k elements in {1, • • • , n + k}, via the maps: Lemma 4.9.Let ψ(x, y) be the number of y-smooth numbers ≤ x.Assume that log(y) ≪ log(x) and log(y) ≫ log log(x).Then log ψ(x, y) x ∼ − log(x) log log(x) log(y) .

Complexity of the search to decision reduction algorithm (algorithm 4.3).
where: • T F S is the execution time of FindSmoothGen (Algorithm 4.1), given by Proposition 4.11: • T tree is the execution time of TreeFill (Algorithm 4.2), given by Proposition 4.12: with s = r m /2 + O(1).• T iso is the time taken in lines 9 and 10 of Algorithm 4.3 to recover the chain of ℓ i -isogenies ψ i : . By Section 2.4, recovering an ℓ i -isogeny from the j-invariants of its domain and codomain costs O(ℓ 2 i polylog(ℓ i ) log(p)) operations over F p 2 .Hence, we have We use Schoof's algorithm [50,Section 5].Namely, we look for primes It follows that: But by Proposition 4.11, we have r m = log(∆)/ log(Bε) with log(1 − ε) ≪ log log(B).We can impose that ε −→ 0, so that log(1 − ε) ≪ log log(B) and that log(ε) ≪ log(B), so that r m ∼ log(∆)/ log(B).Heuristically, the quantity T (∆, B, r m , p) is minimal when the arguments of the two exponentials are close, i.e. when log(∆) log log(∆) log(B) ≃ 2 log(B), the other terms being negligible.Hence, we choose and .
Hence, we can set r m := ⌈ 2 log(∆)/ log log(∆)⌉+1, so that log(ε The space complexity is dominated by the trees T 1 and T 2 , so Algorithm 4.3 uses bits of memory by Proposition 4.12. bits of memory, M (p) being the time complexity of multiplication over F p .

O-orienting problem for quaternion orders
Isogeny problems can often be translated to quaternion problems via the Deuring correspondence, and in many cases the quaternion problems are easier to solve.In this section we consider the quaternion analogue of the O-Orienting Problem stated as follows: Problem 5.1 (Quaternion Order Embedding Problem).Given a maximal quaternion order O ⊂ B p,∞ and an imaginary quadratic order O where an O-orientation of O exists, find the orientation.
Similarly to the curve setting, we define an O-orientation of O to be an embedding ι : O ֒→ O which cannot be extended to a superorder of O, also known as an optimal embedding [54, Chapter 30].
In this section we present a general algorithm, and analyse its complexity, noting special cases.For complexity analysis we assume an efficient factorization oracle exists, however we provide a practical alternative for running the algorithm without such an oracle.For embedding small discriminant quadratic orders O, our algorithm improves the state of the art being efficient up to disc(O) = O(p).
Before moving on to the actual algorithms we give a brief technical overview of the main idea.First we compute a short prime norm N (≈ √ p) connecting ideal between a quaternion order O ′ isomorphic to O and a special extremal order.Our goal is to compute an element of prescribed trace and norm in O ′ and then one can easily construct an element with said trace and norm in O as well.For simplicity assume that the prescribed trace is 0. The trace 0 part of O ′ is a rank 3 lattice and one can compute the Hermite Normal Form (HNF) of this lattice.This means that one has a basis of the form e 11 i + e 12 j + e 13 k, e 22 j + e 23 k, e 33 k and even though e ij are not likely to be integers, their denominator is a divisor of 2N .When looking for an element of trace 0 and norm smaller than p the coefficients of this element with respect to this HNF basis will have a very specific structure.Namely the coefficient of e 11 i + e 12 j + e 13 k will be smaller than p in absolute value and thus can be easily determined by looking at the norm modulo p. Then one only has to work out the two other coefficients which is equivalent to solving a binary quadratic form where the quadratic part is positive definite.This can then essentially be reduced to Cornacchia's algorithm [49].We can extend this to filter out imprimitive solutions.5.1.Finding General Embeddings.First we present an algorithm for finding embeddings, and in the next section we use this to define orientations.Suppose we are given a maximal quaternion order O ⊂ B p,∞ in terms of a Z-basis, and an imaginary quadratic order O = Z[ω], by generator ω of reduced trace t and reduced norm d.
We start with an observation: suppose an embedding ι : Z[ω] ֒→ O exists and let α = ι(ω).Since ω 2 − tω + d = 0 we must also have α 2 − tα+ d = 0. Hence α also has trace t and norm d.Finding any element α of norm d and trace t is enough to define an embedding ι, solving Problem 2.5.This is the approach we take in Algorithm 5.1, finding α ∈ O of a given norm and trace.We make the assumption p = 2 and conventionally use 1, i, j, k as a basis of B p,∞ with i 2 = −q and j 2 = −p.If p ≡ 3 (mod 4) we take q = 1.If p ≡ 5 (mod 8) we take q = 2.If p ≡ 1 (mod 8) we take q to be a prime q ≡ 3 (mod 4) such that p is not a quadratic residue modulo q.While p ≡ 3 mod 4 is the most relevant for isogeny based cryptography, we consider general p.We fix a maximal order O 0 in the following way: Proposition 5.2.[41, Proposition 5.2] The following definitions give a maximal order in B p,∞ for any p = 2: where c is an integer such that q divides c 2 p + 1 where q and c exist by [22, Proposition 1].
We address arbitrary trace in Remark 5.5 and Algorithm 5.1 has no restrictions on trace.However, for simplicity we first describe the algorithm under the assumption that the trace of ω is zero: Step 1: Compute HNF: Put the basis of O into column-style Hermite normal form (HNF).We denote the basis vectors e 0 , e with coefficients e mn ∈ Q.For example see the orders in Proposition 5.2 above.We know the basis is full rank, so e nn = 0 for n = 0, 1, 2, 3, and we prove some additional properties: Lemma 5.3.Given a maximal order O ⊂ B p,∞ with a basis in the above form, the following properties hold (up to finding another basis in right form): (1) e mn ≥ 0 for all n, m

√ p
We return to this in Section 5.3, but for now, by passing to the isomorphic order "closest" to O 0 , we assume that N is of size O( √ p) Step 2: Fix trace: To find a trace zero element α of norm d, we may write an arbitrary element in the following form: Note that since we are working in Hermite Normal Form only e 0 contributes to the trace of α so we set α 0 = 0 to get Tr(α) = 0.
For the condition on the norm, consider the case p ≡ 3 mod 4 for simplicity, however note that this generalizes for any prime p = 2. Then we have the rational ternary quadratic form: (α 1 e 11 ) 2 + p(α 1 e 12 + α 2 e 22 ) 2 + p(α 1 e 13 + α 2 e 23 + α 3 e 33 ) 2 = nrd(α) = d Step 3: Find α 1 mod p: Since α 1 controls the coefficient of i it is the only term without a factor of p. Hence working modulo p removes terms containing α 2 and α 3 , and we can find α 1 ≡ r mod p.
Fix the least positive residue class representative r = r + , as we can execute the remainder of the algorithm a second time on r − if necessary.Then substitute α 1 = r + kp giving a rational ternary quadratic form in k, α 2 and α 3 .
Step 4: A binary quadratic form: As defined in Step 1, we may multiply by the denominator N 2 to obtain integral coefficients.Rearranging we have: where for each iteration, we compute v using the above equation, and with k fixed are left with the integral binary quadratic form v = N 2 (γ 2 1 + γ 2 2 ).
Step 5: Cornacchia's Algorithm: Writing the above form as β 2 1 + β 2 2 = v we solve for integral pairs (β 1 , β 2 ) using Cornacchia's algorithm.For a valid solution we can write it in the form: Note for t odd, this is not optimal as the scaling increases d by a factor of 4, and hence the number of iterations of k by a factor of 2, which can double the running time.Instead we can avoid this by incorporating additional constant terms for the non-zero trace.These details are included in Algorithm 5.1, which we use for our implementation.
The complete algorithm for arbitrary trace is summarised in Algorithm 5.1.Additionally we describe a few further generalisations and improvements: Remark 5.6.Algorithm 5.1 ...
• results in an embedding, but this does not necessarily define an Z[ω]-orientation.This is discussed in Section 5.5.• can be adapted to work with any prime p = 2, not specifically p ≡ 3 mod 4.
For general, B p,∞ = −q,−p Q , q appears in the equations for r, v and the maximum k, and you have to solve β 2 1 + qβ 2 2 = v instead of the sum of two squares.Cornacchia's still works since for B p,∞ , q and p are always coprime.
• can be adapted to non-maximal orders.The value N gains the conductor of the order as a factor.
results and a probabilistic worst case result.We start by attempting to give a worst case running time.Note there are a three reasons why Cornacchia's algorithm may not be efficient at finding all solutions to β 2 1 + qβ 2 2 = v: (1) It requires a factorization of v. To this end, we assume we have an efficient factorization oracle such as Shor's algorithm.See Section 5.3 later on for a practical alternative to using a factorization oracle.(2) Cornacchia's algorithm typically only refers to finding primitive solutions where gcd(β 1 , β 2 ) = 1.To also find imprimitive solutions we must run Cornacchia's on β 2 1 + qβ 2 2 = v/g 2 for every square g 2 | v and scale up the solutions (gβ 1 , gβ 2 ).The number of squares dividing v can be subexponential in v.However, we can say the probability of this for random v is very small, in fact asymptotically there is π 2 6 ∼ 61% chance v is square-free.(3) While just solving for primitive solutions, we must iterate over all the solutions Cornacchia gives.Internally Cornacchia must iterate over all solutions x to the equation x 2 ≡ −q mod v, where the number of solutions can be exponential in v if v has a large number of distinct prime factors.For example, experimentally with p ≡ 3 mod 4 and d ∼ p we get some integers v ∼ p where if v has lots of distinct prime factors, there can be as many as v 0.15 ∼ p 0.15 solutions which is exponential.We resolve this issue by bounding the number of factors of v by the following probability estimate known as the fundamental theorem of probabilistic number theory: • polylog(X) where X is the total size of the inputs, and T (P ) is a value large enough such that the asymptotic probability a random number has less than T (P ) perfect square divisors is larger than P .We define F as the inverse cumulative distribution function of the standard normal distribution where a sample is less than F (P ) with probability P .For example, for P = 0.95 we have F ( P +1 2 ) < 2 and T ( P +1 2 ) ∼ 4. Proof.Steps 1-3 of the algorithm are efficient as polynomial time algorithms exist for computing Hermite normal form [28] [13], and fixing α 0 and solving α 1 modulo p is efficient.In Step 4 a worst case input will result in iterating k over it's full range of values which is O( 1 pe11 d − (α 0 e 00 ) 2 ), where the trace is fixed through α 0 = t 2e00 so (α 0 e 00 ) Then it follows that the number of square roots found in Cornacchia's algorithm is ) ), so we can bound the running time of Cornacchia by O(log(v) (F ( P +1 2 )+1) ) • polylog(v), and clearly for each v we have v ≤ N 2 d < O(pd) and hence polylog(v) = polylog(X).The final consideration is for finding imprimitive solutions using Cornacchia's algorithm which requires repeating for every square dividing v.By definition this is at most T ( P +1 2 ) repetitions with probability P +1 2 .The probability both this condition and v having the correct number of factors is at least P +1 2 + P +1 2 − 1 = P .Now we give a result for the average case running time: Theorem 5.9.Making the following assumptions, regarding iterating over k: • Each v k is distributed like random integers and hence the expected number of distinct prime factors is log log(v k ) by Lemma 5.7, and there is a high probability it only has a few square divisors.• Additionally, the probability each v k is the sum of two squares is independent and at least the probability a random integer less than (N d) 2 is the sum of two squares.• The first solution to Cornacchia's algorithm has β 1 , β 2 uniformly distributed modulo e 22 N and e 33 N respectively.
Then given an efficient factorization oracle, in the case p ≡ 3 mod 4, the average case running time of Algorithm 5.
where X is the total size of all inputs.
We use the following Lemma: Lemma 5.10 (Landau 1908).The number of integers representable as the sum of two squares from from 0 to n ∈ N is the limit C n √ log n as n → ∞, where C ≈ 0.764 is the Landau-Ramanujan constant.Hence for sufficiently large n, the number of integers representable is greater than 1 2 n √ log n .(In fact, experimentally this appears true for all n ≥ 0).
Proof of Theorem 5.9.It's clear the running time is the product of the number of iterations over k, and the running time of Cornacchia's, because all other operations are polynomial time.In the case p ≡ 3 mod 4 we are solving the sum of two squares, hence by Landau, and using the first assumption, we expect (for sufficiently large d) less than 2 log((N d) 2 ) iterations until we find a k where Cornacchia's gives at least one solution.
Now recall that finding one solution to Cornacchia's algorithm is not necessarily enough, since we need to satisfy the conditions α 2 , α 3 ∈ Z.This amounts to checking: Therefore, by the second assumption we expect to have an integral solution after e 22 N × e 33 N solutions from Cornacchias.Noting that e 22 e 33 ≤ N/2 from Prop 5.3, that's N 3 /2 solutions.In total we expect O(N 3 log(N d)) iterations of k.This is bounded above by the maximum number of iterations from Theorem 5.8.Finally, Cornacchia's algorithm uses the efficient factorization oracle to factorize each v k and on average v k is expected to have log log(v k ) distinct prime factors by the first assumption, hence internally Cornacchia's computes at most 2 log log(v k ) = log(v k ) square roots which is efficient.Then to find imprimitive solutions, we only repeat Cornacchia's a constant number of times as the expected number of squares dividing v is very small.Overall this takes time polylog in each v k ≤ N 2 d = O(pd), so this term can be incorperated into polylog(X).
From this we observe the following: Corollary 5.11 (Efficient for orders close to O 0 ).Given an efficient factorization oracle, consider the algorithm applied to the order O 0 .Here we have N = 2, hence the algorithm is efficient; the average case running time is polylog(X).For orders close to O 0 , such as a curve an l-isogeny from the curve with j-invariant 1728, we gain a factor of l in N , hence for small l the algorithm is still efficient.However with each step from O 0 , N gains a factor of the degree of the isogeny, so it gets exponentially harder the further you walk, until we reach the point N ∼ √ p.
For completeness, we now consider the case of arbitrary primes p = 2. Then the quaternion algebra containing order O is B p,∞ = −q,−p Q where q is either 2 or a prime q ≡ 3 mod 4 with Legendre symbol q p = −1.By the same argument as Theorem 5.8, with high probability P the worst case running time is within which is the same as before except the additional factor of 1 √ q appears requiring more iterations over k.Then from Lemma 5.3 part (2), for a different value of K, we get N = 2qN (I) = O(q √ p).Applying this along with Lemma 5.4 gives: Typically q is treated as a constant, so asymptotically the complexity is the same, however, q is actually unbounded; you can construct a prime such that the minimum value for q is larger than a given threshold.Hence we treat q as a variable in our analysis.5.3.Rerandomization and Small Discriminant.Consider the special case of small discriminant orders Z[ω] in Algorithm 5.1.Previously, the best known efficient algorithm, as stated in [57] was simply to look for small vectors in O.This works for | disc(Z[ω])| < 2 √ p − 1 and is stated below.Note that an alternative approach that heuristically works for | disc(Z[ω])| < p 0.8 is outlined in [3] N 2 then the algorithm consists of computing HNF and fast polynomial time arithmetic, and at most M executions of Cornacchia's algorithm.By Cornacchia's algorithm, here we mean mean finding all solutions to x 2 + qy 2 = v, not just primitive ones.
Proof.Naively we use Cornacchia's algorithm once for every k we check.Iterating over k happens twice, once using r + and once using r − , therefore each time we want at most ⌊M/2⌋ iterations.As in proof of Theorem 5.8, and taking into account general p, we can consider the maximum number of iterations over k for each r and bound it by M/2: which, since e 01 ≥ 0 by Prop 5.3, is certainly true if and noting e 11 ≥ 1 N , we get the condition: From this, we obtain a result about disc(Z[ω]) since we can translate generator ω to either − disc(Z The general idea of our rerandomized version is then the common technique of only running Cornacchia on "good" instances.However, if the discriminant is small, then the embedding is with large probability unique, hence we might end up discarding the correct solution.Therefore, we need to rerandomize until all O(1) Cornacchia instances are good, before the algorithm can be sure that no embedding exists.
We define a "good" instance to be x 2 + qy 2 = v k where v k can be factorized in polynomial time, has O(log log(N 2 d)) distinct prime factors, and O(log(N 2 d)) square divisors.The set of prime numbers satisfies these conditions, and with the heuristic that the events of each v k being prime are independent and follow from the density of primes, we expect at most some constant multiple of log(N 2 d) C iterations until C of the Cornacchia instances are primes.With C = O(1) instances from the small discriminant condition, this is efficient.Now we discuss how we rerandomization the order.Let O 0 ⊆ B p,∞ be a maximal order with negligible denominator K (e.g. the "standard" maximal order from Proposition 5.2).As has been pointed out, any maximal order O ⊆ B p,∞ will have denominator bounded by KN , where N is the norm of the connecting ideal from O 0 to O. Hence we can consider any equivalent ideal J = Iγ of small norm, and instead solve the problem in the isomorphic order O R (J) = γ −1 Oγ, before transferring the solution back to O.
As we rerandomize, we might need to try many distinct J ∼ I of small norm.Heuristically, for random orders O, we can expect there to be an abundance of small, equivalent ideals (i.e. with N (J) in O( √ p)).The problem is that this heuristic fails completely if there exists some equivalent ideal I ′ with nrd(I ′ ) ≪ √ p, i.e. in the case that O is too "close" to O 0 .We can fix this by considering other maximal orders O ′ 0 with negligible denominator in other representations of B p,∞ .Specifically, we can generate representations B i = (−q i , −p | Q), where we take q i to be the smallest primes satisfying q i ≡ 3 (mod 4), −q i p = −1.
In each of these representations, regardless of congruence conditions on p, we can take a standard choice of maximal order O 0,i with denominator 2q i as [54, Exercise 15.5].Heuristically, these choices of maximal orders of the different quaternion algebra representations are "independent" in the sense that there is no reason that these should be close to each other, so it is enough to try a small, fixed number of such orders, as (heuristically), the probability that O is close to all O 0,i is negligible.Finally, the explicit isomorphisms B i ∼ = B j are also easy to find and compute, using [25,Lemma 10].The whole algorithm is summarized in Algorithm 5.2.
This gives the following corollary: Corollary 5.17 Now we know primitive embeddings come from primitive solutions, we can determine if an embedding is primitive, and if not extend it to its superorder, very fast using a gcd computation: and to compute a function applying the change of basis transformation taking coefficients of e i s onto coefficients of f i s.This is polynomial time as a variant of the Hermite normal form algorithm, and can be seen as precomputation.Then given a solution α, we change the basis to obtain: Since O is a ring we have 1 ∈ O, and every norm is integral so it cannot contain a rational number less than one.Hence f 0 = f 00 = 1.For α to be a primitive solution there should be no a, b ∈ Z with b > 1 such that α − a = bτ where τ ∈ O.
Equivalently, for any a, when expressing α−a in terms of f i s, the coefficients should not all be divisible by any b > 1.Note that we have: Consider the rational solution x 1 = x 2 = x 3 ∈ Q then as it is rational we can write f (x 1 , x 2 , x 3 ) = wx 2 1 = λ for some w ∈ Q, so |x 1 | = | λ/w|.Since the norm form is positive definite, this means if one variable were to increase, another must decrease in absolute value.Therefore min(|γ 1 |, |γ 2 |, |γ 3 |) ≤ ⌊ λ/w⌋.Now reconsider our integral solution.Suppose the solution is not primitive so gcd(|γ 1 |, |γ 2 |, |γ 3 |) = 1, then there is a prime number ≥ 2 that divides all three numbers.This prime factor must be in the set S = {2, 3, 5, ..., ⌊ λ/w⌋} ∩ {Primes p} as it must divide the smallest of these three numbers.
Heuristically, we assume that γ 1 , γ 2 , γ 3 are distributed like random numbers in the sense that some q ∈ S divides one of them with uniformly random probability 1/q.And assume independence of the probabilities of different factors q 1 , q 2 ∈ S occurring.Then the probability 2 divides all three numbers is 1/2 3 , the probability 3 divides them is 1/3 3 , and the probability any q ∈ S divides them is 1/q 3 .Combined, the probability a number in S divides all three is: P[(γ 1 , γ 2 , γ 3 ) imprimitive] = primes q∈S 1 q 3 ≤ all primes q∈N 1 q 3 = P (3) ≤ 0.175 where P (3) is the prime zeta function at 3. Hence the probability a random solution is primitive is over 80% so if we find 5 independent solutions we would expect at least one to be primitive.Therefore assuming the heuristics above, if we modify algorithm 5.1 to ignore imprimitive solutions, the average running time should only increase by at most a factor of 5.
Note that this argument makes some strong assumptions.Experimentally for some parameters we see the probability the first solution found is primitive is around 80%, but on other parameter choices it is considerably lower.With all parameters we tested, we found the probability is always over 50%, which still suggests the average running time is only worsened by a small factor, but it gives reason to doubt these assumptions.In particular consider independence.Existence of embeddings come with symmetry hence we may find two solutions where there is only a change of sign in the defining formulae.This means the probability q divides a coefficient of one solution might have a strong dependence on whether q divides the coefficient of the second solution.We leave a more complete analysis of the probability of finding primitive embeddings to future research.

Date:
August 23, 2023.This research was funded in part by the UK Engineering and Physical Sciences Research Council (EPSRC) (grant number EP/V011324/1.), by the Agence Nationale de la Recherche under grant ANR MELODIA (ANR-20-CE40-0013), and the France 2030 program under grant agreement No. ANR-22-PETQ-0008 PQ-TLS.This research is also supported by the Hungarian Ministry of Innovation and Technology NRDI Office within the framework of the Quantum Information National Laboratory Program, the János Bolyai Research Scholarship of the Hungarian Academy of Sciences.

Problem 2 . 3 (
Decision O-Orienting Problem).Given an elliptic curve E and an imaginary quadratic order O, determine if E is orientable by O.Problem 2.4 (O-Orienting Problem).Given an elliptic curve E oriented by an imaginary quadratic order O, find the orientation.We explore the following quaternion variant of Problem 2.4 in Section 5.Problem 2.5 (Quaternion Order Embedding Problem).Given a maximal quaternion order O and an imaginary quadratic order O which embeds into O, find the embedding.
Algorithm 2.1 terminates and is correct.It requires O(log(d)) evaluations of the input d-isogeny ϕ on points defined over an extension of degree O(D) of F p 2 and O(D 3 log 3 (p) log(d) + D 10 M (p) log 2 (d)) elementary (bitwise) arithmetic operations, where M (p) is the complexity of the multiplication over F p .Proof.Correctness has been justified by Equation 5 and Lemma 2.10.Termination is clear.Now we compute the complexity.If N is only slightly bigger than d, then s = O(log(d)) and finding a suitable N on line 2 of the algorithm takes O(log(d)) arithmetic operations.

Theorem 3 . 4 .
[30, Theorem 2']  Let O ⊂ B p,∞ be a maximal order in the quaternion algebra ramifying at p and ∞.Let j i : O i ֒−→ O (i ∈ {1, 2}) be two primitive embeddings of orders in the same imaginary quadratic field K

Example 3 . 8 .
Let p = 83 and O = Z[ √ −21], the ring of integers of K = Q( √ −21).We see p is inert in O so K embeds into the quaternion algebra End 0 (E) for any supersingular E over F p .Now, let E/F 2 p be the O-oriented curve y 2 = x 3 + x, we find the orientation by finding an endomorphism ω with N (ω) = 21 = 3 • 7 and T r(ω) = 0. From E we pick a 3-isogeny to y 2 = x 3 + 32x + 38 √ −1, this is also Ooriented.Then we pick a horizontal 7-isogeny which has codomain y 2 = x 3 + 26x.
4. The number of isomorphisms β : E r −→ E 0 is O(1).Using [50, Section 5], we compute the trace of β • ϕ r • • • • • ϕ 1 on line 18 of Algorithm 3.1 in polynomial time in log(p), r = O(log(p)) and max 1≤i≤r ℓ i .The total cost is polynomial in log(p) and max 1≤i≤r ℓ i .4. Solving the O-orienting problem with a decision oracle 4.1.Description of the algorithms.Let O be an imaginary quadratic order with general discriminant ∆ O .Given access to an oracle IsOrientable O for the Decision O-Orienting Problem 2.3, we solve the O-Orienting Problem 2.4 finding a an O-orientation ι : O ֒→ End(E) of any given supersingular elliptic curve E/F p 2 if it exists.

Algorithm 3 . 1 :
Algorithm to solve the O-Orienting Problem 2.4 with an oracle for the Decision O-Orienting Problem 2.3, special discriminant.Data: A supersingular elliptic curve E 0 /F p 2 ; the maximal order O of Q( √ −d), where d := r i=1 ℓ i is a product of small distinct primes, where p > |∆ O | max 1≤i≤r ℓ i ; an oracle IsOrientable O for the Decision O-Orienting Problem 2.3.Result: If E 0 is O-orientable, an efficient representation (as defined in 2.7)

1 )Definition 4 . 2 .
Find a generator θ of O of B-smooth norm N (θ) := r i=1 ℓ i .(2) Starting at E, compute O-oriented (ℓ 1 , • • • , ℓ s )-isogeny tree T 1 and an Ooriented (ℓ s+1 , • • • , ℓ r )-isogeny tree T 2 (with s ≃ r/2).(3) Find a matching leaf in T 1 and T 2 .(4)Extract the corresponding endomorphism ϕ ∈ End(E).(5)Infer from ϕ = ι(θ) an efficient representation of the canonical generator ϕ 0 := ι(ω) (in the sense of Definition 2.7).We explain each step in detail in the following paragraphs.4.1.1.Finding a smooth norm generator.Let O be an imaginary quadratic order and ω be a generator.We want to find another generator θ of O with smooth norm N (θ) = r i=1 ℓ i .The computation of ϕ = ϕ r • • • • • ϕ 1 associated to θ is exponential in the ℓ i and r, so we require the ℓ i and 2 r to be subexponential in log(|∆ O |).For technical reasons (see Lemma 4.3), N (θ) should also be non-square and coprime to ∆ O .In summary, we look for a generator θ of ns-(B, r m , ∆ O )-smooth norm, in the sense of Definition 4.2, with B and 2 rm subexponential in log(|∆ O |).An integer N ∈ N is (B, r m , d)-smooth when its decomposition into prime factors N = r i=1 ℓ i satisfies r ≤ r m , ℓ i ≤ B, and ℓ i ∤ d for all i ∈ {1, • • • , r}.We say that N is ns-(B, r m , d)-smooth when it is (B, r m , d)-smooth and not a square.

Algorithm 4 . 1 :
FindSmoothGen finding a smooth generator of an imaginary quadratic order O. Data: The discriminant ∆ O of an imaginary quadratic order O and smoothness parameters B > 0 and r m ∈ Z >0 .Result: A generator θ of O having ns-(B, r m , ∆ O )-smooth norm N (θ) and its prime factors ℓ 1 , • • • , ℓ r (with multiplicity).

Remark 4 . 6 .
The cases of ∆ O = −3, −4 are excluded from this lemma because in those cases, we have a very simple way to find the orientation:If ∆ O = −3, then O = Z[ζ 3 ], with ζ 3 := (1 + √ −3)/2 so any elliptic curve E that is O-oriented contains an automorphism of order 3.By [52, Theorem III.10.1], we must have j(E) = 0, so E is given by the Weierstrass equation y 2 = x 3 + 1 (up to isomorphism), and ζ 3 corresponds to the automorphism (x, y) ∈ E −→ (ξ 3 x, y) ∈ E, where ξ 3 is a primitive third root of unity in F p 2 .

2 ( 4 ) 2 , and e 11 = 1 Ke22e33Lemma 5 . 4 .
For all e mn , denominators divide K • N (I) where K = 2, 4 or 2q depending on whether p ≡ 3 mod 4, or ≡ 5 mod 8 or ≡ 1 mod 8 respectively, and where I := N OO 0 , N := [O : O ∩ O 0 ] is the connecting ideal from O 0 (3) e 00 = 1 e 22 e 33 ≤ N (I) (5) e 01 = 0 or e 01 = 1/(2Ke 22 e 33 ) where K is defined in (2) Proof.(1) Requirement of HNF.(2) As defined, I is the connecting ideal between O 0 and O.I is contained in both O 0 and O and N (I) • O = I Ī ⊆ O 0 [54, Proposition 16.6.15].Therefore the largest denominator of all e mn s is at most N (I) times the largest denominator of O 0 as given in Prop 5.2.(3) The trace of any element must be integral hence 2e 00 ∈ Z.We must also have 1 ∈ O hence e 00 | 1 so either e 00 = 1 2 or 1 and Tr(O) = Z or 2Z respectively.The (non-reduced) discriminant of any maximal order in B p,∞ is p 2 , so by definition p 2 = det(Tr(e m e n )) ∈ Tr(O), but p is odd, so p 2 ∈ 2Z so we must have e 00 = 1 2 .(4) As O and O 0 are maximal, they both have the same discriminant.Hence the change of basis matrix must have determinant 1 [54, Lemma 15.2.5], which means e nn = f nn = 1 2•K , where (f ) n is the basis of O 0 specified in Proposition 5.2.Then using (3) we have e 11 = 1/(Ke 22 e 33 ).The result follows from (2). (5) 1 ∈ O so there is some n ∈ Z such that 1 e00 e 01 − ne 11 = 0. From above e 00 = 1 so 2e 01 = n Ke22e33 .But to be in HNF we must have already reduced e 01 as much as possible hence n = 0 or 1.In general, we can replace the order O by an isomorphic order O ′ , having denominator bounded by N := K •N (I ′ ) = O( √ p), where I ′ is a connecting (O 0 , O)-ideal equivalent to I, and where K is defined in (2) of Prop 5.3.Such an ideal always exists by the following lemma (taken from [16, Lemma 5.2.2])Let O ⊂ B p,∞ be a maximal order with connecting ideal I = I(O 0 , O), then there exists an equivalent ideal J ∼ I with N (J) ≤ 2 √ 2 π

β 2 = 33 Finally, we must check α 2 , α 3 ∈
N γ 2 = N α 1 e 13 + N α 2 e 23 + N α 3 e 33 and solve for α 2 and α 3 α 2 = β 1 − N α 1 e 12 N e 22 α 3 = β 2 − N α 1 e 13 − N α 2 e 23 N e Z.If this is true we have a valid solution α = α 1 e 1 + α 2 e 2 + α 3 e 3 .If not we continue trying the next solution to Cornacchia's, or move on to the next iteration of k in Step 4. If no solutions are found it means Z[ω] does not embed into O.Remark 5.5 (Arbitrary trace t).Suppose the element we are searching for does not have trace zero.We can always reduce the problem to finding an element of trace zero.Suppose t ∈ 2Z, then since O is a ring we have 1 ∈ O so α−t/2 ∈ O has trace zero and norm d − t 2 /4 ∈ Z.We can search for this trace zero element then translate back to find α.Similarly if t is odd we have trace zero element 2α − t ∈ O of norm 4d − t 2 , once found we translate back, divide by 2 and check α ∈ O in Step 5 of the algorithm.If not we continue searching.

Lemma 5 . 7 (Theorem 5 . 8 .
Erdős-Kac theorem).For a natural number n, the number of distinct prime factors of n follows the standard normal distribution with mean log log n and standard deviation √ log log n as n → ∞.This gives us the following result: Let 0.5 ≤ P < 1.Assuming the heuristic that v is distributed like random integers and hence the number of distinct prime factors follows Lemma 5.7, and given an efficient factorization oracle, the running time of Algorithm 5.1 is within α 1 e 12 ≡ 0 mod e 22 N β 2 − N α 0 e 03 − N α 1 e 13 − N α 2 e 23 ≡ 0 mod e 33 N

.
For arbitrary p = 2, given a maximal order O ⊆ B p,∞ , and a quadratic order O with | disc(O)| in O(p), Algorithm 5.2 heuristically computes an For any element α ∈ O, write α for its class in the lattice O/Z.Consider the discriminant form ∆ : O/Z −→ Q α −→ Tr(α) − 4 nrd(α), an integral quadratic form of rank 3. It does not depend on the choice of a representative α of the class α, and we also write ∆(α).Let d be an integer.We say that a solution α ∈ O of ∆(α) = d is primitive if α is a primitive element of the lattice O/Z, i.e., it is not of the form α = b β for some element β ∈ O/Z and integer b > 1.We now show primitive solutions correspond directly to primitive orientations.Lemma 5.18.An element α ∈ O is a primitive solution of ∆(α) = d if and only if Z[α] ⊆ O is a primitive embedding.Proof.Suppose α is an imprimitive solution of ∆(α) = d, i.e., there are integers a and b > 1 such that (α − a)/b ∈O.Then, Z[α] Z[(α − a)/b] ⊆ O, hence Z[α] ⊆ O is not primitive.Conversely, suppose Z[α] ⊆ O is not primitive, so there exists β ∈ (Q[α] ∩ O) \ Z[α].There exists integers a and b > 1 such that α = a + bβ.In particular, α = b β, so α is not a primitive solution.

Lemma 5 . 19 .
Given a maximal order O with basis e 0 , e 1 , e 2 , e 3 and an element α ∈ O of trace t = Tr(ω) and norm d = nrd(ω), there is a polynomial time algorithm on order O which:• determines whether embedding ι defined by extending ω → α is a primitive or imprimitive embedding of Z[ω],• and if it's imprimitive outputs (a, b, α ′ ) defining a superorder Z[ ω−a b ] ⊃ Z[ω] which ι can be extended to, through the map ω−a b → α ′ .Proof.First convert the upper diagonal basis e 0 , e 1 , e 2 , e 3 of O into an lower diagonal basis f0 , f 1 , f 2 , f 3 , O = f 0 , f 1 , f 2 , f 3 Z = f 00 , f 10 + f 11 i, f 20 i + f 21 i + f 22 j, f 30 + f 31 i + f 32 j + f 33 k Z a 2 , a 3 , a 4 ∈ Z such that a 2 7Compute ϕ(P i,1 ) and ϕ(P i,2 ); Since log(y) ≫ log log(x), we have y ≫ log 2 (x), so that Complexity of the tree filling algorithm (Algorithm 4.2).Proposition 4.12.With inputs B > 0, an imaginary quadratic order O with |∆ O |B < p, primes ℓ 1 , • • • , ℓ s ≤ B splitting in O and an oracle IsOrientable O for Problem 2.3 running in constant time, Algorithm 4.2 runs in time O 2 s B 2 polylog(B) log(p)M (p) , where M (p) is the time complexity of the multiplication in F p .It also uses O(2 s log(p)) bits of memory.Proof.Filling-in tree T in Algorithm 4.2 costs for all 1 ≤ i ≤ s, 2 i−1 calls to IsOrientable O and the computation of 2 i−1 sets of j-invariants ℓ i -isogenous to the same elliptic curve.Each call to IsOrientable O costs O(1) and each j-invariants computation costs O(ℓ 2 i polylog(ℓ i ) log(p)) operations over F p 2 by Section 2.3.Arithmetic operations over F p 2 cost O(M (p)).Hence, the total cost of filling tree T is We already have proved the termination of Algorithm 4.3 when B∆ < p.This is a consequence of Lemma 4.3, Lemma 4.4 and Heuristic 4.7 (which prove that TreeFill and FindSmoothGen terminate).On the whole, the total time complexity of Algorithm 4.3 is 1 , e 2 , e 3 .Then we can write O as: O = e 0 , e 1 , e 2 , e 3 Z = e 00 + e 01 i + e 02 j + e 03 k, e 11 i + e 12 j + e 13 k, e 22 j + e 23 k, e 33 k Z 2 = t 2 4 .And by Prop 5.3 we have 1 e11 ≤ N .Then for each iteration over k, Cornacchia's algorithm is used in Step 5. To be efficient at finding primitive solutions we have to bound the number of distinct prime factors of v, by . Under certain heuristics, we can rerandomize Algorithm 5.1 by considering isomorphic orders, potentially in different representations of B p,∞ , to avoid factoring and bound the denominator N < O( √ p).The result of this is a corollary (Corollary 5.17) which gives a heuristic polynomial time algorithm for solving Problem 2.5 for disc(Z[ω]) in O(p), or deciding that no solution exists.The first step is to bound the number of values to try Cornacchia on in Algorithm 5.1: Lemma 5.16.Fix a positive integer M .As in the context of Algorithm 5.1, take d = nrd(ω), N as the common denominator of the HNF basis, and k the variable we iterate over.Suppose we have Recalling that we can bound N , the denominator of O, to O( √ p), we see that Lemma 5.16 says that when disc(Z[ω]) is in O(p), and assuming q in O(1), the only potentially expensive part in Algorithm 5.1 is Cornacchia on a constant number of instances.