Craig Costello
email: craigco@microsoft.com

Huseyin Hisil
email: huseyin.hisil@yasar.edu.tr

Benjamin Smith
email: smith@lix.polytechnique.fr

Faster Compact Diffie-Hellman: Endomorphisms on the x-line

Keywords: Elliptic curve cryptography, scalar multiplication, twist-secure, side channel attacks, endomorphism, Kummer variety, addition chains, Montgomery curve

may come

Introduction

In this paper, we discuss the design and implementation of state-of-the-art Elliptic Curve Diffie-Hellman key exchange (ECDH) primitives for security level of approximately 128 bits. The major priorities for our implementation are 1. Compactness: We target x-coordinate-only systems. These have the advantage of shorter keys, simple and fast algorithms, and (in a well-designed system) the ability to use arbitrary x-values, not just legitimate x-coordinates of points on a curve (the "illegitimate" values are simply x-coordinates on the quadratic twist). For x-coordinate ECDH, the elliptic curve exists only to supply formulae for scalar multiplications, and a hard elliptic curve discrete logarithm problem (ECDLP) to underwrite a hard computational Diffie-Hellman problem (CDHP) on x-coordinates. The users should not have to verify whether given values correspond to points on a curve, nor should they have to compute any quantity that cannot be derived simply from x-coordinates alone. In particular, neither a user nor an algorithm should have to distinguish between the curve and its quadratic twist-and the curve must be chosen to be twist-secure. 2. Fast, constant-time execution: Every Diffie-Hellman key exchange is essentially comprised of four scalar multiplications,5 so optimizing scalar multiplication P → [m]P for varying P and m is a very high priority. At the same time, a minimum requirement for protecting against side-channel timing attacks is that every scalar multiplication P → [m]P must be computed in constant time (and ideally with the same execution pattern), regardless of the values of m and P .

In our concrete implementation, we target a security level of approximately 128 bits (comparable to Curve25519 [START_REF] Daniel | Curve25519: New Diffie-Hellman speed records[END_REF], secp256r1 [START_REF]Standards for Efficient Cryptography 2 (SEC 2)[END_REF], or brainpoolP256t1 [START_REF] Ecc Brainpool | ECC Brainpool standard curves and curve generation[END_REF]). The reference system with respect to our desired properties is Bernstein's Curve25519, which is based on an efficient, uniform differential addition chain applied to a well-chosen pair of curve and twist. Curve25519 and its twist are presented as Montgomery models. These models not only provide highly efficient group operations, but they are optimized for x-coordinate-only operations; and crucially, they do not distinguish between the curve and its twist. Essentially, well-chosen Montgomery curves offer compactness straight out of the box.

Having chosen Montgomery curves as our platform, we must implement a fast, uniform, and constant-time scalar multiplication on x-coordinates of Montgomery models. To turbocharge our scalar multiplication, we apply a combination of efficiently computable pseudo-endomorphisms and two-dimensional differential addition chains. The use of efficient endomorphisms follows in the tradition of [START_REF] Koblitz | CM-curves with good cryptographic properties[END_REF][START_REF] Solinas | An improved algorithm for arithmetic on a family of elliptic curves[END_REF][START_REF] Gallant | Faster point multiplication on elliptic curves with efficient endomorphisms[END_REF][START_REF] Galbraith | Endomorphisms for faster elliptic curve cryptography on a large class of curves[END_REF], but to the best of our knowledge, this work represents the first use of endomorphism scalar decompositions in the pure x-coordinate setting (that is, without additional input to the addition chain).

Our implementation is built on a curve-twist pair (E, E ′ ) equipped with efficiently computable endomorphisms (ψ, ψ ′ ). The family of Q-curve reductions in [START_REF] Smith | Families of fast elliptic curves from Q-curves[END_REF] offer a combination of fast endomorphisms and compatiblity with fast underlying field arithmetic. Crucially (and unlike earlier endomorphism constructions such as [START_REF] Gallant | Faster point multiplication on elliptic curves with efficient endomorphisms[END_REF][START_REF] Galbraith | Endomorphisms for faster elliptic curve cryptography on a large class of curves[END_REF]), they also offer the possibility of twist-secure group orders over fast fields. One of these curves, with almost-prime order over a 254-bit field, forms the foundation of our construction (see §2). Any other curve from the same family over the same field could be used with only very minor modifications to the formulae below and the source code for our implementations; we explain our specific curve choice in Appendix B. The endomorphisms ψ and ψ ′ induce efficient pseudo-endomorphisms ψ x and ψ ′

x on the x-line; we explain their construction and use in §3. The key idea of this work is to replace single scalar multiplications (m, x(P )) → x([m]P ) with multiscalar multiexponentiations ((a, b), x(P )) -→ x([a]P ⊕ [b]ψ(P )) or x([a]P ⊕ [b]ψ ′ (P )) , where (a, b) is either a short multiscalar decomposition of a random full-length scalar m (that is, such that [m]P = [a]P ⊕ [b]ψ(P ) or [a]P ⊕ [b]ψ ′ (P )), or a random short multiscalar. The choice of ψ or ψ ′ formally depends on whether P is on E or E ′ , but there is no difference between ψ and ψ ′ on the level of x-coordinates: they are implemented using exactly the same formulae. Given that every element of the base field is the x-coordinate of a point on E or E ′ , we may view the transformation above as acting purely on field elements, and not curve points.

From a practical point of view, the two crucial differences compared with conventional ECDH over a 254-bit field are 1. The use of 128-bit multiscalars (a, b) in Z 2 in place of the 254-bit scalar m in Z. We treat the geometry of multiscalars, the distribution of their corresponding scalar values, and the derivation of constant-bitlength scalar decompositions in §4.

The use of two-dimensional differential addition chains to compute x([a]P ⊕[b]ψ(P ))

given only (a, b) and x(P ). We detail this process in §5.

We have implemented three different two-dimensional differential addition chains: one due to Montgomery [START_REF] Peter | Evaluating recurrences of form Xm+n = f (Xm, Xn, Xm-n) via Lucas chains[END_REF] via Stam [START_REF] Stam | Speeding up subgroup cryptosystems[END_REF], one due to Bernstein [START_REF] Daniel | Differential addition chains[END_REF], and one due to Azarderakhsh and Karabina [START_REF] Azarderakhsh | A new double point multiplication algorithm and its application to binary elliptic curves with endomorphisms[END_REF]. We provide implementation details and timings for scalar multiplications based on each of our chains in §6. Each offers a different combination of speed, uniformity, and constant-time execution. The differential nature of these chains is essential in the xcoordinate setting, which prevents the effective use of the vector chains traditionally used in the endomorphism literature (such as [START_REF] Ernst | Addition chains of vectors[END_REF]).

We will be putting detailed Magma code online. A complete mixed-assembly-and-C code is publicly available in eBATS [START_REF]eBACS: ECRYPT Benchmarking of Cryptographic Systems[END_REF] format at http://hhisil.yasar.edu.tr/files/hisil20131026compact.tar.gz

The curve

We begin by defining our curve-twist pair (E, E ′ ). We work over the field F p 2 = F p (i), where

p := 2 127 -1 and i 2 = -1 .
We have chosen this Mersenne prime for its compatibility with a range of fast techniques for modular arithmetic, including Montgomery-and NIST-style approaches. We build efficient F p 2 -arithmetic on top of the fast F p -arithmetic described in [START_REF] Bos | Fast cryptography in genus 2[END_REF]. Appendix A provides a complete description of our arithmetic routines.

In what follows, it will be convenient to define the constants

u := 1466100457131508421 , v := 1 2 (p -1) = 2 126 -1 , and w := 1 4 (p + 1) = 2 125 .
The curve E and its twist E ′ . We define E to be the elliptic curve over F p 2 with affine Montgomery model E :

y 2 = x(x 2 + Ax + 1) ,
where

A = A 0 + A 1 • i with A 0 = 45116554344555875085017627593321485421 , A 1 = 2415910908 .
The element 12/A is not a square in F p 2 , so the curve over F p 2 defined by

E ′ : (12/A)y 2 = x(x 2 + Ax + 1)
is a model of the quadratic twist of E. The twisting

F p 4 -isomorphism δ : E → E ′ is defined by δ : (x, y) -→ (x, (A/12) 1/2 y). The map δ 1 : (x, y) → (x W , y W ) = ( 12 A x + 4, 12 2 A 2 y) defines an F p 2 -isomorphism between E ′ and the Weierstrass model E 2,-1,s : y 2 W = x 3 W + 2(9(1 + si) -24)x W -8(9(1 + si) -16) of [32, Theorem 1] with s = i(1 -8/A 2 ) = 86878915556079486902897638486322141403
, so E is a Montgomery model of the quadratic twist of E 2,-1,s (in the notation of [32, §5] we have

E ∼ = E ′ 2,-1,s and E ′ ∼ = E 2,-1,s ). All of these curves have j-invariant j(E) = j(E ′ ) = j(E 2,-1,s ) = 2 8 (A 2 -3) 3 A 2 -4 = 2 6 (5 -3si) 3 (1 -si) (1 + s 2 ) 2 .
Group orders, structures, and generators. Using the SEA algorithm [START_REF] Schoof | Counting points on elliptic curves over finite fields[END_REF], we find that

#E(F p 2 ) = 4N and #E ′ (F p 2 ) = 8N ′
where

N = v 2 + 2u 2 and N ′ = 2w 2 -u 2
are 252-bit and 251-bit primes, respectively. Looking closer, we see that

E(F p 2 ) ∼ = Z/2Z × Z/2Z × Z/N Z and E ′ (F p 2 ) ∼ = Z/2Z × Z/4Z × Z/N ′ Z .
Recall that every element of F p 2 is either the x-coordinate of two points in E(F p 2 ), the xcoordinate of two points in E ′ (F p 2 ), or the x-coordinate of one point of order two in both E(F p 2 ) and E ′ (F p 2 ). The x-coordinates of the points of exact order 2 in E(F p 2 ) (and in E ′ (F p 2 )) are 0 and -

1 2 A ± 1 2 √ A 2 -4; the points of exact order 4 in E ′ (F p 2 ) have x-coordinates ±1.
Either of the points with x-coordinate 2 will serve as a generator for the cryptographic subgroup

E(F p 2 )[N ]; either of the points with x-coordinate 2 -i generate E ′ (F p 2 )[N ′ ].
Curve points, x-coordinates, and random bitstrings. Being Montgomery curves, both E and E ′ are compatible with the Elligator 2 construction [6, §5]. For our curves, [6, Theorem 5] defines efficiently invertible injective maps F p 2 → E(F p 2 ) and F p 2 → E ′ (F p 2 ). This allows points on E and/or E ′ to be encoded in such a way that they are indistinguishable from uniformly random strings of bitlength 254. Since we work with x-coordinates only in this article, a square root is saved when computing the injection (see [6, §5.5] for more details).

The ECDLP on E and E ′ . Suppose we want to solve an instance of the DLP in E(F p 2 ) or E ′ (F p 2 ). Applying the Pohlig-Hellman-Silver reduction [START_REF] Pohlig | An improved algorithm for computing logarithms over GF(p) and its cryptographic significance. Information Theory[END_REF], we almost instantly reduce to the case of solving a DLP instance in either E(F p 2 )[N ] or E ′ (F p 2 )[N ′ ]. The best known approach to solving such a DLP instance is Pollard's rho algorithm [START_REF] Pollard | Monte carlo methods for index computation (mod p)[END_REF], which (properly implemented) can solve DLP instances in

E(F p 2 )[N ] (resp. E ′ (F p 2 )[N ′ ]) in around 1 2 √ πN ∼ 2 125.8 (resp. 1 2 √ πN ′ ∼ 2 125.
3 ) group operations on average [START_REF] Daniel J Bernstein | On the correct use of the negation map in the Pollard rho method[END_REF]. One might expect that working over F p 2 would imply a √ 2-factor speedup in the rho method by using Frobenius classes; but this seems not to be the case, since neither E nor E ′ is a subfield curve (cf. [36, §6]).

The embedding degrees of E and E ′ with respect to N and N ′ are 1 50 (N -1) and 1 2 (N ′ -1), respectively, so ECDLP instances in E(F p 2 )[N ] and E(F p 2 )[N ′ ] are not vulnerable to the Menezes-Okamoto-Vanstone [START_REF] Menezes | Reducing elliptic curve logarithms to logarithms in a finite field[END_REF] or Frey-Rück [START_REF] Frey | The Tate pairing and the discrete logarithm applied to elliptic curve cryptosystems[END_REF] attacks. The trace of E is p 2 +1-4N = ±1, so neither E nor E ′ are amenable to the Smart-Satoh-Araki-Semaev attack [START_REF] Smart | The discrete logarithm problem on elliptic curves of trace one[END_REF][START_REF] Satoh | Fermat quotients and the polynomial time discrete log algorithm for anomalous elliptic curves[END_REF][START_REF] Semaev | Evaluation of discrete logarithms in a group of p-torsion points of an elliptic curve in characteristic p[END_REF].

While our curves are defined over a quadratic extension field, this does not seem to reduce the expected difficulty of the ECDLP when compared to elliptic curves over similar-sized prime fields. Taking the Weil restriction of E (or E ′ ) to F p as in the Gaudry-Hess-Smart attack [START_REF] Gaudry | Constructive and destructive facets of Weil descent on elliptic curves[END_REF], for example, produces a simple abelian surface over F p ; and the best known attacks on DLP instances on simple abelian surfaces over F p offer no advantage over simply attacking the ECDLP on the original curve [START_REF] Smart | How secure are elliptic curves over composite extension fields[END_REF][START_REF] Gaudry | Index calculus for abelian varieties of small dimension and the elliptic curve discrete logarithm problem[END_REF]. (See [14, §9] for further discussion.)

Superficially, E is what we would normally call twist-secure (in the sense of Bernstein [3] and Fouque-Réal-Lercier-Valette [START_REF] Fouque | Fault attack on elliptic curve Montgomery ladder implementation[END_REF]), since its twist E ′ has a similar security level. Indeed, E (and the whole class of curves from which it was drawn) was designed with this notion of twist-security in mind. However, twist-security is more subtle in the context of endomorphismbased scalar decompositions; we will return to this subject in §4 below.

The endomorphism ring. Let π E denote the Frobenius endomorphism of E. The curve E is ordinary (its trace t E is prime to p), so its endomorphism ring is an order in the quadratic field K := Q(π E ). (The endomorphism ring of an ordinary curve and its twist are always isomorphic, so what holds below for E also holds for E ′ .) We will see below that E has an endomorphism ψ such that

ψ 2 = -[2]π E . The discriminant of Z[ψ] is the fundamental discriminant D K = -8 • 5 • 397 • 10528961 • 6898209116497 • 1150304667927101 of K, so Z[ψ] is the maximal order in K; hence, End(E) = Z[ψ].
The safecurves specification [START_REF] Daniel | SafeCurves: choosing safe curves for elliptic-curve cryptography[END_REF] suggests that the discriminant of the CM field should have at least 100 bits; our E easily meets this requirement, since D K has 130 bits. We note that well-chosen GLS curves can also have large CM field discriminants, but GLV curves have tiny CM field discriminants by construction: for example, the endomorphism ring of the curve secp256k1 [START_REF]Standards for Efficient Cryptography 2 (SEC 2)[END_REF] (at the heart of the Bitcoin system) has discriminant -3.

Brainpool [START_REF] Ecc Brainpool | ECC Brainpool standard curves and curve generation[END_REF] requires the ideal class number of K to be larger than 10 6 ; this property is never satisfied by GLV curves, which have tiny class numbers (typically ≤ 2) by construction. But E easily meets this requirement: the class number of End(E) is

h(End(E)) = h(D K ) = 2 7 • 31 • 37517 • 146099 • 505117 ∼ 10 19 .
3 The endomorphism on E and pseudo-endomorphisms on the x-line Theorem 1 of [START_REF] Smith | Families of fast elliptic curves from Q-curves[END_REF] defines an efficient endomorphism

ψ 2,-1,s : (x W , y W ) -→ -x p W 2 - 9(1 -si) x p W -4 , y p W √ -2 -1 2 + 9(1 -si) (x p W -4) 2
of degree 2p on the Weierstrass model E 2,-1,s , with kernel (4, 0) . To avoid an ambiguity in the sign of the endomorphism, we must fix a choice of √ -2 in F p 2 . We choose the "small" root:

√ -2 := 2 64 • i . (1)
Applying the isomorphisms δ and δ 1 , we define efficient F p 2 -endomorphisms

ψ := (δ 1 δ) -1 ψ 2,-1,s δ 1 δ and ψ ′ := δψδ -1 = δ -1 1 ψ 2,-1,s δ 1
of degree 2p on E and E ′ , respectively, each with kernel (0, 0) . More explicitly: if we let

n(x) := A p A x 2 + Ax + 1 , d(x) := -2x , m(x) := n ′ (x)d(x) -n(x)d ′ (x), r(x) := A p A (x 2 -1) , and s(x) := n(x) p d(x) p ,
then ψ and ψ ′ are defined (using the same value of √ -2 fixed in Eq. ( 1)) by

ψ : (x, y) → s(x) , - 12 v A v y p √ -2 m(x) p d(x) 2p and ψ ′ : (x, y) → s(x) , - √ -2 12 2v A 2v y p r(x) p d(x) 2p .
The action of the endomorphisms on curve points. Theorem 1 of [START_REF] Smith | Families of fast elliptic curves from Q-curves[END_REF] tells us that

ψ 2 = -[2]π E and (ψ ′ ) 2 = [2]π E ′ , (2)
where π E and π E ′ are the p 2 -power Frobenius endomorphisms of E and E ′ , respectively, and

P (ψ) = P (ψ ′ ) = 0 , where P (T ) = T 2 -4uT + 2p . (3)
If we restrict to the cryptographic subgroup E(F p 2 )[N ], then ψ must act as multiplication by an integer eigenvalue λ, one of the two roots of P (T ) modulo N . Similarly,

ψ ′ acts on E ′ (F p 2 )[N ′ ]
as multiplication by one of the roots λ ′ of P (T ) modulo N ′ . The correct eigenvalues are

λ ≡ - v u (mod N ) and λ ′ ≡ - 2w u (mod N ′ ) .
Equation ( 2) implies that λ 2 ≡ -2 (mod N ) and λ ′2 ≡ 2 (mod N ′ ). (Note that choosing the other square root of -2 in Eq. ( 1) negates ψ, ψ ′ , λ, λ ′ , and u.)

To complete our picture of the action of ψ on E(F p 2 ) and ψ ′ on E ′ (F p 2 ), we describe its action on the points of order 2 and 4 listed above. Choosing square roots appropriately:

(0, 0) -→ 0 under ψ and ψ ′ , -1 2 A ± 1 2 √ A 2 -4, 0 -→ (0, 0) under ψ and ψ ′ , 1, ± 1 2 A(A + 2)/3 -→ -1 2 A -1 2 √ A 2 -4, 0 under ψ ′ , -1, ± 1 2 -A(A + 2)/3 -→ -1 2 A + 1 2 √ A 2 -4, 0 under ψ ′ . (4)
Pseudo-endomorphisms on the x-line. One advantage of the Montgomery model is that it allows a particularly efficient arithmetic using only the x-coordinate. Technically speaking, this corresponds to viewing the x-line P 1 as the Kummer variety of E: that is, P 1 ∼ = E/ ±1 . The x-line is not a group: if P and Q are points on E, then x(P ) and x(Q) determine the pair {x(P ⊕ Q), x(P ⊖ Q)}, but not the individual elements x(P ⊕ Q) and x(P ⊖ Q). However, the x-line inherits part of the endomorphism structure of E: every endomorphism φ of E induces a pseudo-endomorphism 6 φ x : x → φ x (x) of P 1 , which determines φ up to sign; and if φ 1 and φ 2 are two endomorphisms of E, then

(φ 1 ) x (φ 2 ) x = (φ 2 ) x (φ 1 ) x = (φ 1 φ 2 ) x = (φ 2 φ 1 ) x .
In addition 7 to pseudo-doubling (DBL), pseudo-addition (ADD), and combined pseudodoubling and pseudo-addition (DBLADD) on P 1 , we need expressions for both ψ x and (ψ ± 1) x to initialise the addition chains in Section 5. Moving to projective coordinates, write x = X/Z and y = Y /Z; then the negation map on E is [-1] : (X : Y : Z) → (X : -Y : Z), and the double cover

E → E/ [±1] ∼ = P 1 is (X : Y : Z) → (X : Z). The pseudo-doubling on P 1 is [2] x ((X : Z)) = (X + Z) 2 (X -Z) 2 : (4XZ) (X -Z) 2 + A+2 4 • 4XZ . (5)
Our endomorphism ψ induces the pseudo-endomorphism

ψ x ((X : Z)) = A p (X -Z) 2 -A+2 2 (-2XZ) p : A(-2XZ) p .
Composing ψ x with itself and comparing with Eq. ( 5), we confirm that ψ x ψ x = - [START_REF] Babai | On Lovász' lattice reduction and the nearest lattice point problem[END_REF] x (π E ) x .

6 "Pseudo-endomorphisms" are, strictly speaking, endomorphisms of the x-line P 1 . We use the term pseudoendomorphism here to avoid confusion with endomorphisms of elliptic curves, and also to fit with the use of terms like "pseudo-addition" for operations on the x-line. 7 Montgomery's formulae for these operations are available in explicit form on the EFD [START_REF] Daniel | Explicit-formulas database[END_REF].

Proposition 1. With the notation above, and with the value of √ -2 chosen as in Eq. ( 1),

(ψ ± 1) x (x) = (ψ ′ ± 1) x (x) = 2s 2 nd 4p -x(xn) p m 2p A p-1 ∓ 2(1/A) (p-1)/2 m p s(xn) (p+1)/2 d 2p A p-1 √ -2 2s(x -s) 2 d 4p A p-1 . (6)
Proof. If P and Q are points on a Montgomery curve By 2 = x(x 2 + Ax + 1), then

x(P ± Q) = B (x(P )y(Q) ∓ x(Q)y(P )) 2 x(P )x(Q) (x(P ) -x(Q)) 2 .
Taking P = (x, y) to be a generic point on E (where B = 1), setting Q = ψ(P ), and using y 2 = -A p 2A dn to eliminate y yields the expression for (ψ ± 1) x above. The same process for E ′ (with B = 12 A ), eliminating y with 12 A y 2 = -A p 2A dn, yields the same expression for (ψ ′ ±1) x . ⊓ ⊔ Deriving explicit formulae to compute the expression(s) in ( 6) is straightforward, so we omit them for space considerations (but see our code online). If P ∈ E, then on input of x(P ), the combined computation of the three projective elements (X λ-1 : Z λ-1 ), (X λ : Z λ ), (X λ+1 : Z λ+1 ), which respectively correspond to the three affine elements x([λ-1]P ), x([λ]P ), x([λ + 1]P ), incurs 15 multiplications, 129 squarings and 10 additions in F p 2 . The bottleneck of this computation is raising dn to the power of (p + 1)/2 = 2 126 , which incurs 126 squarings. We note that squarings are significantly faster than multiplications in F p 2 (see Appendix A).

Scalar decompositions

We want to evaluate scalar multiplications [m]P as [a]P ⊕ [b]ψ(P ), where m ≡ a + bλ (mod N ) and the multiscalar (a, b) has a significantly shorter bitlength8 than m. For our applications we impose two extra requirements on multiscalars (a, b), so as to add a measure of side-channel resistance:

1. both a and b must be positive, to avoid branching and to simplify our algorithms; and 2. the multiscalar (a, b) must have constant bitlength (that is, bitlength independent of m as m varies over Z), so that multiexponentiation can run in constant time.

In some protocols-notably Diffie-Hellman-we are not interested in the particular values of our random scalars, as long as those values remain secret. In this case, rather than starting with m in Z/N Z (or Z/N ′ Z) and finding a short, positive, constant-bitlength decomposition of m, it would be easier to randomly sample some short, positive, constant-bitlength multiscalar (a, b) from scratch. The sample space must be chosen to ensure that the corresponding distribution of values a + bλ in Z/N Z does not make the discrete logarithm problem of finding a + bλ appreciably easier than if we started with a random m.

The lattices of zero decompositions. The problems of finding good decompositions and sampling good multiscalars are best addressed using the geometric structure of the spaces of decompositions for E and E ′ . The multiscalars (a, b) such that a + bλ ≡ 0 (mod N ) or a + bλ ′ ≡ 0 (mod N ′ ) form lattices L = (N, 0), (-λ, 1)

and

L ′ = (N ′ , 0), (-λ ′ , 1) , respectively, with a + bλ ≡ c + dλ (mod N ) if and only if (a, b) -(c, d) is in L (similarly, a + bλ ′ ≡ c + dλ ′ (mod N ′ ) if and only if (a, b) -(c, d) is in L ′ ).
L = e 1 := (v, u) , e 2 := (-2u, v)
and

L ′ = e ′ 1 := (u, w) , e ′ 2 := (2u -2w, 2w -u) .
2 ] is • ∞ -reduced if b 1 ∞ ≤ b 2 ∞ ≤ b 1 -b 2 ∞ ≤ b 1 + b 2 ∞ . These conditions are satisfied for [b 1 , b 2 ] = [e 2 , -e 1 ] and [e ′ 1 , e ′ 2 ], so ( e 2 ∞ , e 1 ∞ ) (resp. ( e ′ 1 ∞ , e ′ 2
A := (a, b) ∈ R 2 : 0 ≤ vb -ua < N and 0 ≤ 2ub + va < N , A ′ := (a, b) ∈ R 2 : 0 ≤ ub -wa < N ′ and 0 ≤ (2u -2w)b -(2w -u)a < N ′ .
Every integer m has precisely one decomposition for

E(F p 2 )[N ] (resp. E ′ (F p 2 )[N ′ ]) in any translate of A by L (resp. A ′ by L ′ ).
9 That is, reduced with respect to Kaib's generalized Gauss reduction algorithm [START_REF] Kaib | The Gauss lattice basis reduction algorithm succeeds with any norm[END_REF] for the infinity norm. 10 For the Euclidean norm, the bases [e1, e2] and [e ′ 1 , 2e

) ∞ ≤ 1 2 ( e 1 ∞ + e 2 ∞
). This decomposition is approximately the shortest possible, in the sense that the true shortest decomposition is at most ±e 1 ± e 2 away.

Observe that e 1 ∞ = e 2 ∞ = 2 126 -1, so (ã, b) has bitlength at most 126.

However, ã or b may be negative (violating the positivity requirement), or have fewer than 126 bits (violating the constant bitlength requirement). Indeed, m → (ã, b) maps Z onto (A - 1 2 (e 1 + e 2 )) ∩ Z 2 . This region of the (a, b)-plane, "centred" on (0, 0), contains multiscalars of every bitlength between 0 and 126-and the majority of them have at least one negative component. We can achieve positivity and constant bitlength by adding a carefully chosen offset vector from L, translating (A - 1 2 (e 1 + e 2 )) ∩ Z 2 into a region of the (a, b)-plane where every multiscalar is positive and has the same bitlength. Adding 3e 1 or 3e 2 ensures that the first or second component always has precisely 128 bits, respectively; but adding 3(e 1 + e 2 ) gives us a constant bitlength of 128 bits in both. Theorem 3 makes this all completely explicit. Twist-security with endomorphisms. We saw in §2 that DLPs on E and its twist E ′ have essentially the same difficulty, while Proposition 4 shows that the real DLP instances presented to an adversary by 127-bit multiscalar multiplications are not biased into a significantly more attackable range. But there is an additional subtlety when we consider the fault attacks considered in [START_REF] Daniel | Curve25519: New Diffie-Hellman speed records[END_REF] and [START_REF] Fouque | Fault attack on elliptic curve Montgomery ladder implementation[END_REF]: If we try to compute [m]P for P on E, but an adversary sneaks in a point P ′ on the twist E ′ instead, then in the classical context the adversary can derive m after solving the discrete logarithm [m mod

N ′ ]P ′ in E ′ (F p 2 ).

Two-dimensional differential addition chains

Addition chains are used to compute scalar multiplications using a sequence of group operations (or pseudo-group operations). A one-dimensional addition chain computes [m]P for a given integer m and point P ; a two-dimensional addition chain computes [a]P ⊕ [b]Q for a given multiscalar (a, b) and points P and Q. In a differential addition chain, the computation of any ADD P ⊕ Q is always preceded (at some earlier stage in the chain) by the computation of its associated difference P ⊖ Q. The simplest differential addition chain is the original one-dimensional "Montgomery ladder" [START_REF] Peter | Speeding the Pollard and elliptic curve methods of factorization[END_REF], which computes scalar multiplications [m]P for a single exponent m and point P . Every ADD in the Montgomery ladder is in the form [i]P ⊕ [i + 1]P , so every associated difference is equal to P . Several two-dimensional differential addition chains have been proposed, targeting multi-exponentiations in elliptic curves and other primitives; we suggest [START_REF] Daniel | Differential addition chains[END_REF] and [START_REF] Stam | Speeding up subgroup cryptosystems[END_REF] for overviews.

In any two-dimensional differential chain computing [a]P ⊕ [b]Q for general P and Q, the input consists of the multiscalar (a, b) and the three points P , Q, and P ⊖ Q. The initial difference P ⊖Q (or equivalently, the initial sum P ⊕Q) is essential to kickstart the chain on P and Q, since otherwise (by definition) P ⊕ Q cannot appear in the chain. As we noted in §1, computing this initial difference is an inconvenient obstruction to pure x-coordinate only multiexponentiations on general input: the pseudo-group operations ADD, DBL, and DBLADD can all be made to work on x-coordinates (the ADD and DBLADD operations make use of the associated differences that are available in a differential chain), but in general it is impossible to compute the initial difference x(P ⊖ Q) in terms of x(P ) and x(Q).

For our application, we want to compute x([a]P ⊕ [b]ψ(P )) given inputs (a, b) and x(P ). Crucially, we can compute x(P ⊖ ψ(P )) as (ψ -1) x (x(P )) using Proposition 1; this allows us to compute x([a]P ⊕ [b]ψ(P )) using two-dimensional differential addition chains with input (a, b), x(P ), ψ x (x(P )), and (ψ -1) x (x(P )).

We have implemented one one-dimensional differential addition chain (LADDER) and three two-dimensional differential addition chains (PRAC, AK, and DJB). We briefly describe each chain, with its relative benefits and drawbacks, below.

(Montgomery) LADDER chains. We implemented the full-length one-dimensional Montgomery ladder, as a reference to assess the speedup that our techniques offer over conventional scalar multiplication (It is also used as a subroutine within PRAC). The ladder can be made constant-time by adding a suitable multiple of N to the input scalar.

PRAC chains. Montgomery [START_REF] Peter | Evaluating recurrences of form Xm+n = f (Xm, Xn, Xm-n) via Lucas chains[END_REF] proposed a number of algorithms for generating differential addition chains that are often much shorter than his eponymous ladder. His one-dimensional "PRAC" routine contains an easily-implemented two-dimensional subroutine, which computes the double-exponentiation [a]P ⊕ [b]Q very efficiently. The downside for our purposes is that the chain is not uniform: different inputs (a, b) give rise to different execution patterns, rendering the routine vulnerable to a number of side-channel attacks. Our implementation of this chain follows Algorithm 3.25 of [START_REF] Stam | Speeding up subgroup cryptosystems[END_REF] AK chains. Azarderakhsh and Karabina [START_REF] Azarderakhsh | A new double point multiplication algorithm and its application to binary elliptic curves with endomorphisms[END_REF] recently constructed a two-dimensional differential addition chain which offers some middle ground in the trade-off between uniform execution and efficiency. While it is less efficient than PRAC, their chain has the advantage that all but one of the iterations consist of a single DBLADD; this uniformity may be enough to thwart some simple side-channel attacks. The single iteration which does not use a DBLADD requires a separate DBL and ADD, and this slightly slower step can appear at different stages of the algorithm. The location of this longer step could leak some information to a side-channel adversary under some circumstances, but we can protect against this by replacing all of the DBLADDs with separate DBL and ADDs, incurring a very minor performance penalty. A more serious drawback for this chain is its variable length: the total number of iterations depends on the input multiscalar. This destroys any hope of achieving a runtime that is independent of the input. Nevertheless, depending on the physical threat model, this chain may still be a suitable alternative. Our implementation of this chain follows Algorithm 1 in [START_REF] Azarderakhsh | A new double point multiplication algorithm and its application to binary elliptic curves with endomorphisms[END_REF].

DJB chains. Bernstein gives the fastest known two-dimensional differential chain that is both fixed length and uniform [4, §4]. This chain is slightly slower than the PRAC and AK chains, but it offers stronger resistance against many side-channel attacks. 13 If the multiscalar (a, b) has bitlength ℓ, then this chain requires precisely ℓ-1 iterations, each of which computes one ADD and one DBLADD. In our case, Theorem 3 allows us to fix the number of iterations at 127. The execution pattern of the multiexponentiation is therefore independent of the input, and will run in constant time. It takes some work to organise the description in [START_REF] Daniel | Differential addition chains[END_REF] into a concrete algorithm; we give an algorithm specific to our chosen curve in Appendix C.

Operation counts. Table 1 profiles the number of high-level operations required by each of our addition chain implementations on E. We used the decomposition in Theorem 3 to guarantee positive constant-bitlength multiscalars. In situations where side-channel resistance is not a priority, and the AK or PRAC chain is preferable, variable-length decompositions could be used: these would give lower operation counts and slightly faster average timings.

Table 1. Pseudo-group operation counts per scalar multiplication on the x-line of E for the 2-dimensional DJB, AK and PRAC chains (using endomorphism decompositions) and the 1dimensional LADDER. The counts for LADDER and DJB are exact; those for PRAC and AK are averages, with corresponding standard deviations, over 10 6 random trials (random scalars and points). In addition to the operations listed here, each chain requires a final F p 2 -inversion to convert the result into affine form. The LADDER and DJB chains offer some slightly faster high-level operations. In these chains, the "difference elements" fed into the ADDs are fixed; if these points are affine, then this saves one F p 2 -multiplication for every ADD. In LADDER, the difference is always the affine x(P ), so these savings come for free. In DJB, the difference is always one of the four values x(P ), ψ x (x(P )), or (ψ ± 1) x (x(P )), so a shared inversion is used to convert ψ x (x(P )) and (ψ ± 1) x (x(P )) from projective to affine coordinates. While this costs one F p 2 -inversion and six-F p 2 multiplications, it saves 253 F p 2 -inversions inside the loop.

Timings

Table 2 lists cycle counts for our implementations run on an Intel Core i7-3520M (Ivy Bridge) processor at 2893.484 MHz with hyper-threading turned off, over-clocking ("turboboost") disabled, and all-but-one of the cores switched off in BIOS. The implementations were compiled with gcc 4.6.3 and tested on a 64-bit Linux environment. All cycle counts were obtained using the SUPERCOP toolkit [START_REF]eBACS: ECRYPT Benchmarking of Cryptographic Systems[END_REF].

The most meaningful performance comparison we can draw is with Bernstein's Curve25519 software. Like our software, Curve25519 works entirely on the x-line, from start to finish; using the uniform one-dimensional Montgomery ladder, it runs in constant time. Thus, fair performance comparisons can only be made between his implementation and the two of ours that are also both uniform and constant-time: LADDER and DJB. Benchmarked on our hardware with all settings as above, Curve25519 scalar multiplications ran in 182,000 cycles on average. Looking at Table 2, we see that using the one-dimensional LADDER on the Table 2. Performance timings for four different implementations of compact, x-coordinateonly scalar multiplications targeting the 128-bit security level. Timings are given for the one-dimensional Montgomery ladder, as well as the two-dimensional (DJB, AK and PRAC) chains that benefit from the application of an endomorphism and subsequent short scalar decompositions. x-line of E gives a factor 1.20 speed up over Curve25519, while combining an endomorphism with the two-dimensional DJB chain on the x-line of E gives a factor 1.26 speed up over Curve25519.

While there are several other implementations targeting the 128-bit security level that give faster performance numbers than ours, we reiterate that our aim was to push the boundary in the arena of implementations that are purely x-coordinate-only.

Hamburg [START_REF] Hamburg | Fast and compact elliptic-curve cryptography[END_REF] has also documented a fast software implementation employing xcoordinate-only Montgomery arithmetic. However, it is difficult to compare Hamburg's software with ours: his is not available to be benchmarked, and his figures were obtained on the Sandy Bridge architecture (they were also manually scaled back to compensate for turboboost being enabled). Nevertheless, Hamburg's own comparison with Curve25519 suggests that a fair comparison between our constant-time implementations and his would be close.

A Efficient arithmetic in F p and F p 2

We access lower level integer arithmetic for efficient addition, subtraction, multiplication and squaring operations in F p and F p 2 where p = 2 127 -1, see §2. At this level, elements of F p are represented by integer values in the usual way, with the exception that the representation of 0 is not unique: 0 is allowed to be represented in "semi-reduced" form by the integers 0 and p. Semi-reduced values can be used in any chain of operations without causing an exception, since all of our algorithms are designed to accept inputs and produce outputs in the interval [0, p]. The implementor should reduce each output into the range [0, p) at the very end of the target computation, in order to satisfy unique representation field elements. This type of arithmetic has already been exploited in earlier works, such as [START_REF] Bos | Fast cryptography in genus 2[END_REF], but a thorough exposition has not yet appeared.

We will be frequently referring back to the divisibility lemma of integers.

Lemma 5. Let u, v ∈ Z with v > 0. Then there exist unique q, r ∈ Z such that u = r + qv and 0 ≤ r < v. In particular, q = ⌊u/v⌋ and r = u -⌊u/v⌋v where ⌊.⌋ is the floor function.

In what follows, the "mod 2 128 " and "mod 2 256 " operators, are included (even though they are often unnecessary) to reinforce the fact that all arithmetic operations are being performed on an unsigned integer arithmetic circuit over a 128-bit data type. We let k i denote the i th significant bit of an integer k and use (k i , . . . , k j ) to denote the integer formed by the bit-string that starts with k i , continues with bits in order of increasing significance, and ends with k j (with 0 ≤ i ≤ j ≤ 127). This is a Big-endian type representation. Although it is possible to provide much shorter arguments for sections A.1-5, we prefer to keep the notes in longer format in order to assist easier verification.

It should be noted that all of the techniques in this section avoid branching. This is highly desirable for an efficient implementation, especially on an architecture with pipelining capability.

A.1 Semi-reduced addition modulo p

The operation (a + b) mod p is replaced by Algorithm 1. -Line-3: A semi-reduced output is given by f := (d+e) mod 2 128 , observing that 0 ≤ f ≤ p.

Input: a, b ∈ Z such that 0 ≤ a, b ≤ p. Output: f ∈ Z such that f ≡ (a + b) (mod p) and 0 ≤ f ≤ p.
A.2 Semi-reduced subtraction modulo p

The operation (a -b) mod p is replaced by Algorithm 2.

Input: a, b ∈ Z such that 0 ≤ a, b ≤ p.

Output: f ∈ Z such that f ≡ (a -b) (mod p) and 0 ≤ f ≤ p. Output: f ∈ Z such that f ≡ (ab) (mod p) and 0 ≤ f ≤ p. -Line-2: Use Lemma 5 to write c = d + 2 127 (e + 2 127 f ) for integers 0 ≤ d, e < 2 127 and f . There are two cases to investigate:

• Case 1: Assume that a b + âb < (p + 1) 2 . Then 0 ≤ c < (p + 1) -Line-3: Set g := (e + f ) mod 2 128 where 0 ≤ g ≤ p.

-Line-4: Noting that d + 2 127 (e + 2 127 f ) ≡ d + 2 127 g ≡ d + g (mod p), that 0 ≤ d, g ≤ p, and that 0 ≤ d + g ≤ 2p, a semi-reduced output is obtained by Algorithm 1 applied on the operands d and g.

A.5 Lazy semi-reduction modulo p following a double-word subtraction

The lazy reduction (ab -âb ) mod p is replaced by Algorithm 5.

Input: a, â, b, b ∈ Z such that 0 ≤ a, â, b, b ≤ p.
Output: h ∈ Z such that h ≡ (ab -âb ) (mod p) and 0 ≤ h ≤ p. -Line- in F p . This means that, for any fixed value of A 1 in F p , there are at most four corresponding possibilities for A 0 ∈ F p , which can be found as the roots of a quartic equation over F p (and vice versa). Keeping in mind that we need to search over many values of A to find a twistsecure pair (E, E ′ ), the upshot of ( 7) is that we cannot expect to find a cryptographically suitable curve with both A 0 and A 1 being small. The F p 2 -arithmetic described in Appendix A places no preference on which of these two coefficients is small, so we flipped a coin and conducted our search choosing A 1 to be small 14 : we restricted our search to A 1 values with integer representations less than 2 32 , so that A 1 would only occupy one word on 32bit and 64-bit platforms. Our search also prioritised A 1 whose integer representations had low signed Hamming-weight, in the hope that multiplication by A 1 might possibly be faster when computed via sequence of additions and shifts. When the integer representation of A 1 had a signed Hamming-weight of 1,2 or 3, we did not find any curve-twist pairs with optimal cofactors; however, when further loosening the search to include those A 1 with a corresponding signed Hamming-weight of 4, we found 10 such pairs. Among these, there are 3 whose A 1 had an integer representation of precisely 32 bits -the curve in Section 2 corresponds to the smallest such A 1 . Although the low signed Hamming-weight of A 1 did not end up being useful in our implementation, its small size did give rise to a minor but noticeable speedup. The takeaway message is that the construction in [32, §5] is flexible enough to find a vast number of twist-secure curves over any quadratic extension field, for which all of the discussion in this paper can either be directly applied or easily modified (regardless of how the parameter search is designed). Such curve-twist pairs can be readily found in a verifiably random manner, following, for instance, the method described in [11, §5].

For example: let H be a hash function, and let π(i) be the string consisting of the first i digits of π (without the decimal), i.e. π(1) = "3", π(2) = "31", π(3) = "314", and so on. Starting with i = 1, we conducted some other searches by taking s := H(π(i)) mod p until we found the first curve-twist pair E and E ′ with optimal cofactors 4 and 8. With H = SHA -1, the first pair was found at i = 19244; and, with H = SHA -256, the first pair was found at i = 41004. Our discussion and software can easily be modified to handle such cases.

C Bernstein's uniform two-dimensional differential addition chain

Algorithm 6 is a concrete adaptation of Bernstein's addition chain [4, §4] to our curve E, following the multiscalar decomposition described in Section 4. We use the usual formulae (see [START_REF] Daniel | Explicit-formulas database[END_REF]) for pseudo-doubling, pseudo-addition, and for the combination of the two, writing their inputs and outputs as follows. For pseudo-doubling, we write

x([2]R) = DBL(x(R)) ; for pseudo-addition, we write x(T ⊕ U ) = ADD(x(T ), x(U ), x(T ⊖ U )) ; and for their combined computation, we write

x([2]R), x(R ⊕ S) = DBLADD(x(R), x(S), x(R ⊖ S)) . 14 We also included the restriction that the integer representation of A1 be congruent to 2 modulo 4, since it is actually the constant (A + 2)/4 which is used in Montgomery's formulae [22, p. 261].

The main iterations in the chain compute a DBLADD alongside a standalone ADD, so we denote combined pseudo-doubling and pseudo-addition by

x([2]R), x(R ⊕ S), x(T ⊕ U ) = DBLDBLADD(x(R), x(S), x(R ⊖ S), x(T ), x(U ), x(T ⊖ U )) .

The chain is determined in its entirety using only bit operations before any arithmetic is done on the x-line (the symbols ⊕ and ⊖ denote bit operations in Lines 3-9 of Algorithm 2, but curve operations everywhere else).

Observe that the associated differences in pseudo-additions are always one of the four affine input points x(P ), x(Q), x(P ⊖ Q), or x(P ⊕ Q). On the other hand, the three running values T 0 = (X 0 : Z 0 ), T 1 = (X 1 : Z 1 ) and T 2 = (X 2 : Z 2 ) are projective. Thus, the final step (which chooses one of the three running values to output) will involve an F p 2 -inversion to output T ind f inal in affine form.

Proof.= 1 4 ẽ′ 2

 42 The proof of [32, Proposition 2] constructs sublattices ẽ1 := -2(v, u), ẽ2 := -2(2u, v) of index 4 in L and ẽ′ 1 := 2(2w, -u), ẽ′ 2 := 4(u, w) of index 8 in L ′ . We easily verify that e 1 = -1 2 ẽ2 and e 2 = -1 2 ẽ1 are both in L; since ẽ1 , ẽ2 has index 4 in e 1 , e 2 , we must have L = e 1 , e 2 . Similarly, both e ′ 1 and e ′ 2 = 1 2 (ẽ ′ 2 -ẽ′ 1 ) are in L ′ , and thus form a basis for L ′ . According to [19, Definition 3], an ordered lattice basis [b 1 , b

  ∞ )) are the successive minima of L (resp. L ′ ) by [19, Theorem 5]. 10 ⊓ ⊔ In view of Proposition 2, the fundamental parallelograms of L and L ′ are the regions of the (a, b)-plane defined by

Theorem 3 . 5 √

 35 Given an integer m, let (a, b) be the multiscalar defined by a := m + (3 -⌊α⌉) v -2 (3 -⌊β⌉) u and b := (3 -⌊α⌉) u + (3 -⌊β⌉) v , where α and β are the rational numbers α := (v/N )m and β := -(u/N )m . Then 2 127 < a, b < 2 128 , and m ≡ a + bλ (mod N ). In particular, (a, b) is a positive decomposition of m, of bitlength exactly 128, for any m. Proof. We have m ≡ a + bλ (mod N ) because (a, b) = (ã, b) + 3(e 1 + e 2 ) ≡ (m, 0) (mod L), where (ã, b) is the translate of (m, 0) by the Babai roundoff ⌊α⌉e 1 + ⌊β⌉e 2 described above. Now (ã, b) lies in A -1 2 (e 1 + e 2 ), so (a, b) lies in A + 5 2 (e 1 , e 2 ); our claim on the bitlength of (a, b) follows because the four "corners" of this domain all have 128-bit components. ⊓ ⊔ Random multiscalars. As we remarked above, in a pure Diffie-Hellman implementation it is more convenient to simply sample random multiscalars than to decompose randomly sampled scalars. Proposition 4 shows that random multiscalars of at most 127 bits correspond to reasonably well-distributed values in Z/N Z and in Z/N ′ Z, in the sense that none of the values occur more than one more or one fewer times than the average, and the exceptional values are in O( √ N ). Such multiscalars can be trivially turned into constant-bitlength positive 128-bit multiscalars-compatible with our implementation-by (for example) completing a pair of 127-bit strings with a 1 in the 128-th bit position of each component. Proposition 4. Let B = [0, p] 2 ; we identify B with the set of all pairs of strings of 127 bits. 1. The map B → Z/N Z defined by (a, b) → a+bλ (mod N ) is 4-to-1, except for 4(p-6u+4) ≈ 4 √ 2N values in Z/N Z with 5 preimages in B, and 8(u 2 -3u + 2) ≈ 1 N values in Z/N Z with only 3 preimages in B.

2 . 7 √

 27 The map B → Z/N ′ Z defined by (a, b) → a + bλ ′ (mod N ′ ) is 8-to-1, except for 8u 2 ≈ 2 N ′ values with 9 preimages in B.Proof (Sketch). For (1): the map (a, b) → a + bλ (mod N ) defines a bijection between each translate of A ∩ Z 2 by L and Z/N Z. Hence, every m in Z/N Z has a unique preimage (a 0 , b 0 ) in A ∩ Z 2 , so it suffices to count ((a 0 , b 0 ) + L) ∩ B for each (a 0 , b 0 ) in A ∩ Z 2 . Cover Z 2 with translates of A by L; the only points in Z 2 that are on the boundaries of tiles are the points in L. Dissecting B along the edges of translates of A and reassembling the pieces, we see that 8v -24u + 20 < 4p multiscalars in B occur with multiplicity five, 8u 2 -24u + 16 < p/9 with multiplicity three, and every other multiscalar occurs with multiplicity four. There are therefore 4N + (8v -24u + 20) -(8u 2 -24u + 16) = (p + 1) 2 preimages in total, as expected. The proof of (2) is similar to (1), but counting ((a, b) + L ′ ) ∩ B as (a, b) ranges over A ′ . ⊓ ⊔

  But in the endomorphism context, we compute [m]P as [a]P ⊕ [b]ψ(P ), and the attacker sees [a + bλ ′ ]P ′ , which is not [m mod N ′ ]P ′ (or even [a + bλ mod N ′ ]P ′ ); we should ensure that the values (a + bλ ′ mod N ′ ) are not concentrated in a small subset of Z/N ′ Z when (a, b) is a decomposition for E(F p 2 )[N ]. This can be achieved by a similar argument to that of Proposition 4: the map Z/N Z → Z/N ′ Z defined by m → (a, b) → a + bλ ′ (mod N ′ ) is a good approximation of a 2-to-1 mapping.

  11 : given a multiscalar (a, b) and points P , Q, and P -Q, this algorithm computes d = gcd(a, b) and R = [ a d ]P ⊕ [ b d ]Q. To finish computing [a]P ⊕ [b]Q, we write d = 2 i e with i ≥ q and e odd, then compute S = [2 i ]R with i consecutive DBLs, before finally computing [e]S with a one-dimensional LADDER chain 12 .

  chain

1 cAlgorithm 1 :• Case 1 :• Case 2 :

 1112 := (a + b) mod 2 128 ; 2 d := (c0, c1, . . . , c126), e := (c127); 3 f := (d + e) mod 2 128 ; 4 return f ; Semi-reduced addition modulo p -Line-1: Notice that 0 ≤ c = a + b ≤ 2p < 2 128 . -Line-2: Use Lemma 5 to write c = d + 2 127 e for integers 0 ≤ d < 2 127 and e. There are two cases to investigate: Assume that a + b ≤ p. The bounds on c and d imply that 0/2 127 ≤ c/2 127 = (d + 2 127 e)/2 127 = d/2 127 + 2 127 e/2 127 = e ≤ p/2 127 , so e = 0. Thus a + b ≡ d + 2 127 e ≡ d + 2 127 • 0 ≡ d + 0 ≡ d + e (mod p). Assume that a + b > p. Then p < c ≤ 2p. The bounds on c and d imply that (p + 1)/2 127 ≤ e ≤ 2p/2 127 , so e = 1. The bounds on c also imply that p -2 127 < c -2 127 ≤ 2p -2 127 and we have d = c -2 127 e = c -2 127 , so 0 ≤ d < p. Thus a + b ≡ d + 2 127 e ≡ d + 2 127 • 1 ≡ d + 1 ≡ d + e (mod p).

1Algorithm 2 :-Line- 1 :

 21 c := (a -b) mod 2 128 ; 2 d := (c0, c1, . . . , c126), e := (c127); 3 f := (d -e) mod 2 128 ; 4 return f ; Semi-reduced subtraction modulo p Notice that 0 ≤ c < 2 128 . -Line-2: Use Lemma 5 to write c = d + 2 127 e for integers 0 ≤ d < 2 127 and e. There are two cases to investigate: • Case 1: Assume that a ≥ b. Then 0 ≤ c = a -b ≤ p. The bounds on c and d imply that 0/2 127 ≤ c/2 127 = (d + 2 127 e)/2 127 = e ≤ p/2 127 , so e = 0. Thus a -b ≡ d + 2 127 e ≡ d -e (mod p). • Case 2: Assume that a < b. Then c = 2 128 +a-b and -p ≤ a-b < 0. So, 2 127 < c < 2 128 . The bounds on c and d imply that (2 127 + 1)/2 127 ≤ e ≤ (2 128 -1)/2 127 , so e = 1. The bounds on c also imply that 2 127 -2 127 < c -2 127 < 2 128 -2 127 , and we have d = c -2 127 e = c -2 127 . So, 0 < d ≤ p and d ≥ e. Thus a -b ≡ (2 128 + a -b) -2 128 ≡ c -2 128 ≡ d + 2 127 e -2 128 ≡ d -e (mod p). Line-3: A semi-reduced output is given by f := (d-e) mod 2 128 , observing that 0 ≤ f ≤ p. A.3 Semi-reduced multiplication modulo p The operation (ab) mod p is replaced by Algorithm 3. -Line-1: Notice that 0 ≤ c = ab ≤ p 2 < 2 256 . -Line-2: Use Lemma 5 to write c = d + 2 127 e for integers 0 ≤ d < 2 127 and e. The bounds on c and d imply that 0/2 127 ≤ c/2 127 = (d + 2 127 e)/2 127 = e ≤ p 2 /2 127 , so 0 ≤ e < p. -Line-3: Noting that ab ≡ d + 2 127 e ≡ d + (2 127 -1)e + e ≡ d + pe + e ≡ d + e (mod p), that 0 ≤ d, e ≤ p, and that 0 ≤ d + e ≤ 2p, a semi-reduced output is obtained by Algorithm 1 applied on the operands d and e. Input: a, b ∈ Z such that 0 ≤ a, b ≤ p.

1 cAlgorithm 3 :

 13 := (ab) mod 2 256 ; 2 d := (c0, c1, . . . , c126), e := (c127, c128, . . . , c253); 3 f := semi-add(d, e); 4 return f ; Semi-reduced multiplication modulo p A.4 Lazy semi-reduction modulo p following a double-word addition The lazy reduction (a b + âb) mod p is replaced by Algorithm 4. Input: a, â, b, b ∈ Z such that 0 ≤ a, â, b, b ≤ p.Output: h ∈ Z such that h ≡ (a b + âb) (mod p) and 0 ≤ h ≤ p.

1Algorithm 4 :

 4 c := (a b + âb) mod 2 256 ; 2 d := (c0, c1, . . . , c126), e := (c127, c128, . . . , c253), f := (c254); 3 g := (e + f ) mod 2 128 ; 4 h := semi-add(d, g); 5 return h; Lazy semi-reduction modulo p following a double-word addition -Line-1: Notice that 0 ≤ c = a b + âb ≤ 2p 2 < 2 256 .

1 cAlgorithm 5 :

 15 := (ab -âb ) mod 2 256 ; 2 d := (c0, c1, . . . , c126), e := (c127, c128, . . . , c253), f := (c254), g := (c255); 3 h := (e -f ) mod 2 128 ; 4 j := semi-add(d, g); 5 return j; Lazy semi-reduction modulo p following a double-word subtraction -Line-1: Notice that 0 ≤ c < 2 256 .

  Let a, â, b, b ∈ Z and 0 ≤ a, â, b, b ≤ p. We use the obvious method which computes (a + âi) + (b + bi) as ((a + b) mod p) + ((â + b) mod p)i. Both modular additions are replaced by Algorithm 1. Analogous comments apply for the case of subtraction which uses Algorithm 2. A.7 Multiplication in F p 2 Let a, â, b, b ∈ Z and 0 ≤ a, â, b, b ≤ p. On the target architecture, we experienced the best performance for computing (a+âi)(b+ bi) by coupling a Karatsuba-based operation scheduling

  The sets of decompositions of m for E(F p )[N ] and E(F p 2 )[N ′ ] therefore form lattice cosets

	(m, 0) + L	and	(m, 0) + L

′

, respectively; we can find short decompositions of m for E(F p )[N ] (resp. E(F p 2 )[N ′ ]) by subtracting vectors near (m, 0) in L (resp. L ′ ) from (m, 0). To find these vectors, we need

• ∞ -reduced 9 bases for L and L ′ .

Proposition 2 (Definition of e 1 , e 2 , e ′ 1 , e ′ 2 ). Up to order and sign, the shortest possible bases for L and L ′ (with respect to • ∞ ) are given by

  Short, constant-bitlength scalar decompositions. Returning to the problem of finding short decompositions of m: let (α, β) be the (unique) solution in Q 2 to the system αe 1 +βe 2 = (m, 0). Since e 1 , e 2 is reduced, the closest vector to (m, 0) in L is one of the four vectors ⌊α⌋e 1 + ⌊β⌋e 2 , ⌊α⌋e 1 + ⌈β⌉e 2 , ⌈α⌉e 1 + ⌊β⌋e 2 , or ⌈α⌉e 1 + ⌈β⌉e 2 by[START_REF] Kaib | The Gauss lattice basis reduction algorithm succeeds with any norm[END_REF] Theorem 19]. Following Babai[START_REF] Babai | On Lovász' lattice reduction and the nearest lattice point problem[END_REF], we subtract ⌊α⌉e 1 + ⌊β⌉e 2 from (m, 0) to get a decomposition (ã, b) of m; by the triangle inequality, (ã, b

′ 1e ′ 2 ] are • 2-reduced, but [e ′ 1 , e ′ 2 ]

is not.

  ≤ 2p 2 . The bounds on c, d, and e imply that (p + 1) 2 /(2 127 ) 2 ≤ f ≤ 2p 2 /(2 127 ) 2 , so f = 1. The bounds on c also imply that (p + 1) 2 -(2 127 ) 2 ≤ c -(2 127 ) 2 ≤ 2p 2 -(2 127) 2 , and we have d + 2 127 e = c -(2 127 ) 2 f = c -(2 127 ) 2 . So, 0 ≤ d + 2 127 e ≤ ((p -1) 2 -2). The bounds on d + 2 127 e imply that 0/2 127 ≤ (d + 2 127 e)/2 127 ≤ ((p -1) 2 -2)/2 127 , so 0 ≤ e < (p -2). Thus a b + âb ≡ d + 2 127 (e + 2 127 f ) ≡ d + 2 127 (e + 2 127 • 1) ≡ d + 2 127 (e + 1) ≡ d + 2 127 (e + f ) (mod p) and 0 ≤ e + f < p.

2

. The bounds on c, d, and e imply that 0/(2 127

) 2 ≤ c/(2 127 ) 2 = (d + 2 127 e + (2 127 ) 2 f )/(2 127 ) 2 = f ≤ ((p + 1) 2 -1)/(

2

127 ) 2 , so f = 0. Thus a b + âb ≡ d + 2 127 (e + 2 127 f ) ≡ d + 2 127 (e + 2 127 • 0) ≡ d + 2 127 (e + 0) ≡ d + 2 127 (e + f ) (mod p) and 0 ≤ e + f < p. • Case 2: Assume that a b + âb ≥ (p + 1) 2 . Then (p + 1) 2 ≤ c

  2: Use Lemma 5 to write c = d + 2 127 (e + 2 127 (f + 2g)) for integers 0 ≤ d, e < 2 127 , 0 ≤ f < 2, and g. There are two cases to investigate: • Case 1: Assume that ab ≥ âb . Then 0 ≤ c = ab -âb ≤ p 2 . The bounds on c, d, e and f imply that 0/(2 127 ) 2 ≤ c/(2 127 ) 2 = (d + 2 127 e + (2 127 ) 2 (f + 2g))/(2 127 ) 2 = f +2g ≤ p 2 /(2 127 ) 2 ; that is f +2g = 0. So, f = g = 0. Thus d+2 127 (e+2 127 (f +2g)) ≡ d + 2 127 (e + 2 127 • 0) ≡ d + 2 127 (e -0) ≡ d + 2 127 (e -f ) (mod p). Noting that d + 2 127 (e + 2 127 (f + 2g)) ≡ d + 2 127 h ≡ d + h (mod p), that 0 ≤ d, h ≤ p, and that 0 ≤ d + h ≤ 2p, a semi-reduced output is obtained by Algorithm 1 applied on the operands d and h.

	A.6 Addition and subtraction in F p 2

• Case 2: Assume that ab < âb . Then c = 2 256 + ab -âb and -p 2 ≤ ab -âb < 0. So, 2 256 -p 2 ≤ c < 2 256 . As in the previous case, the bounds on c, d, e and f imply that

(2 256 -p 2 )/(2 127 ) 2 ≤ f + 2g ≤ (2 256 -1)/(2 127 ) 2 , so f + 2g = 3 and f = g = 1.

The bounds on c also imply that

2 256 -p 2 -3(2 127 ) 2 = 2 128 -1 ≤ c -3(2 127 ) 2 < 2 256 -3(2 127

) 2 and we also have d + 2 127 e = c -(2 127 ) 2 (f + 2g) = c -3(2 127 ) 2 . So, 2 128 -1 ≤ d + 2 127 e. The bounds on d + 2 127 e imply that (2 128 -1)/2 127 ≤ (d + 2 127 e)/2 127 < (2 256 -3(2 127 ) 2 )/2 127 = 2 127 , so 1 ≤ e < 2 127 and e ≥ f . Thus ab -âb ≡ (2 256 + ab -âb ) -2 256 = c -2 256 ≡ c -4 ≡ d + 2 127 (e + 2 127 (f + 2g)) -4 ≡ d + 2 127 (e + 2 127 (1 + 2 • 1)) -4 ≡ d + 2 127 (e -1) ≡ d + 2 127 (e -f ) (mod p).

-Line-3: Set h := (e -f ) mod 2 128 where 0 ≤ h ≤ p.

-Line-4:

We do not count the cost of authenticating keys, etc., here. Two of the scalar multiplications can be computed in advance in the case of static Diffie-Hellman; in this fixed-base scenario (where the point P is constant but m varies) we could profit from extensive precomputations. For simplicity, in this work we concentrate on the dynamic case (where P and m are variable).

The bitlength of a scalar m is ⌈log 2 |m|⌉; the bitlength of a multiscalar (a, b) is ⌈log 2 (a, b) ∞⌉.

We implemented the binary version of Montgomery's two-dimensional PRAC chain, neglecting the ternary steps in [23,Table 4] (see also[START_REF] Stam | Speeding up subgroup cryptosystems[END_REF] Table 3.1]). Including these ternary steps could be significantly faster than our implementation, though it would require fast explicit formulae for tripling on Montgomery curves.

In practice d is very small, so there is little benefit in using a more complicated chain for this final step.

It would be interesting to implement our techniques with Bernstein's non-uniform two-dimensional extendedgcd differential addition chain[START_REF] Daniel | Differential addition chains[END_REF]. This chain can outperform two-dimensional PRAC, though it "takes more time to compute and is not easy to analyse".

Acknowledgements The authors would like to thank Joppe W. Bos for independently benchmarking our code on his computer. The second author acknowledges that the notes of Appendix A grew from discussions with Joppe W. Bos on an earlier work [10].

with two lazy reductions. This computes the product as (ab -âb ) mod p + (a + â)(b + b) -ab -âb mod p i.

The routine starts with two integer additions t 0 := a + â and t 1 := b + b satisfying 0 ≤ t 0 , t 1 < (2 128 -1). The routine continues with the 3 integer multiplications t 2 := t 0 t 1 , t 3 := ab and t 4 := âb satisfying 0 ≤ t 2 ≤ (2 128 -2) 2 < 2 256 and 0 ≤ t 3 , t 4 ≤ (2 127 -1) 2 < 2 254 . Since t 2 > t 3 and (t 2 -t 3 ) > t 4 , the integer value t 5 := (t 2 -t 3 ) -t 4 is positive and satisfies both 0 ≤ t 5 ≤ 2p 2 < 2 255 and t 5 = a b + âb. The reduction of t 5 is performed as in Algorithm 4. The reduction of t 6 := (t 3 -t 4 ) mod 2 256 is performed as in Algorithm 5.

A.8 Squaring in F p 2

Let a, â, b, b ∈ Z and 0 ≤ a, â, b, b ≤ p. On the target architecture, we experienced that a lazy semi-reduction strategy gives the same timings as the (non-lazy) semi-reduction strategy for computing (a + âi) 2 = (a -â)(a + â) + 2aâ i.

A.9 Other operations in F p 2 Many other F p 2 operations can be efficiently performed by F p arithmetic only. For instance, negation can be performed as -a = (0 -a) + (0 -â)i, p-th powering as a p = a + (0 -â)i, and inversion as a -1 = a(a 2 + â2 ) p-2 + (0 -â(a 2 + â2 ) p-2 )i -our F p 2 -inversion implementation incurs 128 F p -squarings, 12 F p -multiplications and 2 F p -additions/subtractions.

B How was this curve chosen?

The particular curve-twist pair described in Section 2 was chosen from the family of degree-2 Q-curve reductions described in [START_REF] Smith | Families of fast elliptic curves from Q-curves[END_REF]Proposition 3]. This family (over F p 2 ) is parameterised by a free parameter s in F p , and covers at least p -3 different isomorphism classes for any given p. Every curve in the family is equipped with an efficient endomorphism, and the arithmetic properties of the family make it possible to find twist-secure curves in the family.

We chose our curve E from this family as follows. First, we fixed p = 2 127 -1; a Mersenne prime, this p facilitates very fast modular arithmetic. Next, we chose a tiny nonsquare to define

); this makes F p 2 -arithmetic slightly easier, and our formulae much simpler. We then needed to find a twist-secure curve E in the family with a rational transformation into Montgomery form. This means that E had to have order divisible by 4, so optimal group order for E was therefore #E = 4N , with N prime; it follows (from p 2 ≡ 1 mod 4) that the optimal group order for E ′ is then #E ′ = 8N ′ , with N ′ is prime. Our curve search therefore discarded parameter values corresponding to pairs that did not have these optimal cofactors.

To optimise performance, we searched for parameter values s in F p that gave rise to Montgomery representations whose constant A in F p 2 had "small" coefficients: that is, for A = A 0 + A 1 i with A 0 and A 1 in F p , we wanted the integer representation of A 0 and A 1 to be small. However, the construction in [32, §5] requires the constant A to satisfy 8/A 2 = 1 + si, from which it plainly follows that

Input: a, b ∈ Z + (both 128 bits -see Theorem 3), and x(P ), x(Q), x(Q ⊖ P ), x(Q ⊕ P ) (four affine elements on the x-line, where Q = ψ(P ) on E) Output: x([a]P ⊕ [b]ψ(P ))

1 initialization: (a)2 = (a127, . . . , a0) ∈ {0, 1} 128 , (b)2 = (b127, . . . , b0) ∈ {0, 1} 128 . 2 z0, z1, z2, z3 ← ().

/* z's start as empty bit-sequences */

* zj = (zj,0, . . . , zj,126) ∈ {0, 1} 127 , j = 0, . . . 3 */ 18 case [0, 0, 0, 0] : T1, T0, T2 ← DBLADDADD(T1, T0, x(Q ⊕ P ), T2, T1, x(Q)). case [0, 0, 0, 1] : T1, T0, T2 ← DBLADDADD(T1, T0, x(Q ⊕ P ), T2, T1, x(P )). Algorithm 6: Bernstein's uniform 2-D chain, tailored to the curve in Section 2.