Ping Ngai Chung

Craig Costello
email: craigco@microsoft.com

Benjamin Smith
email: smith@lix.polytechnique.fr

Craig Costello

Fast, uniform, and compact scalar multiplication for elliptic curves and genus 2 Jacobians with applications to signature schemes

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

In terms of per-bit security, elliptic curves and Jacobians of genus 2 curves appear to be roughly equivalent. However, when it comes to efficient and sidechannel-aware implementations, we see a curious divergence. Full genus 2 Jacobian arithmetic is relatively slow compared with elliptic curves, and hard to implement in a uniform and constant-time fashion. On the other hand, we have an extremely fast and uniform scalar pseudomultiplication for genus 2 Kummer surfaces, which often outperforms its elliptic equivalent; but this comes at the cost of identifying elements with their inverses, so Kummer surfaces are widely believed to be suitable only for Diffie-Hellman protocols, where no individual full group operations are required.

Thus, genus 2 beats elliptic curve performance for Diffie-Hellman (where we can use Kummer surfaces), but is beaten in almost every other implementation scenario (where we are confined to Jacobians). In this article we show how to exploit Kummer surface arithmetic to carry out full Jacobian scalar multiplications, bringing the speed and side-channel security of Kummers to implementations of other discrete-log-based cryptographic protocols. In particular, our results show that many existing competitive Diffie-Hellman implementations based on pseudomultiplication can be wrapped with simple and efficient preand post-computations to yield competitive full scalar multiplication algorithms, ready for use in more general cryptosystems, including signature schemes.

Scalar multiplication. To make things precise, let G be a subgroup of an elliptic curve or a genus 2 Jacobian (with ⊕ denoting the group law, and ⊖ its inverse). We are interested in computing one-dimensional scalar multiplications (m, P ) -→ [m]P := P ⊕ • • • ⊕ P m times for m ∈ Z ≥0 , P ∈ G , which is the operation at the heart of all discrete logarithm and Diffie-Hellman problem-based cryptosystems. We are also interested in two-dimensional multiscalar multiplications

((m, n), (P, Q)) -→ [m]P ⊕ [n]Q for m, n ∈ Z ≥0 , P, Q ∈ G ;
these appear explicitly in many cryptographic protocols, including signature verification, but they are also a key ingredient in endomorphism-accelerated one-dimensional scalar multiplication techniques such as GLV [START_REF] Gallant | Faster point multiplication on elliptic curves with efficient endomorphisms[END_REF] and its descendants.

If the scalar m is secret, then [m]P must be computed in a uniform and constant-time way to protect against even the most elementary side-channel attacks. This means that the execution path of the algorithm must be independent of the scalar m (we may assume that the bitlength of m is fixed).

Uniform, constant-time algorithms for elliptic curve scalar multiplication are well-known (and even widely-used). In contrast, if G is a subgroup of a genus 2 Jacobian, then this requirement has represented an insurmountable obstacle until now: the usual Cantor arithmetic [START_REF] Cantor | Computing in the Jacobian of a hyperelliptic curve[END_REF] and its derivatives [START_REF] Hisil | Jacobian coordinates on genus 2 curves[END_REF] have so many incompatible special cases that it has appeared impossible to implement it in a uniform way without abandoning all hope of competitive efficiency (see §6).

Kummer surfaces and x-lines. The situation changes dramatically when we pass to the quotient G/ ±1 , identifying elements with their inverses. Let

x : G -→ G/ ±1 be the quotient map, so x(P ) = x(⊖P ) for all P in G. In the elliptic case, x is projection onto the x-coordinate (whence our notation); in genus 2, x is the map from the Jacobian to its Kummer surface. We emphasize that G/ ±1 is not a group, and at first glance this prevents us instantiating group-based protocols in G/ ±1 . Nevertheless, scalar multiplication in G induces a welldefined pseudomultiplication (m, x(P )) -→ x([m]P ) in G/ ±1 , because [m](±P ) = ±([m]P ). This suffices for implementing protocols like Diffie-Hellman key exchange which only involve scalar multiplication, and not individual group operations (which we lose in passing from G to G/ ±1 ).

Pseudomultiplication algorithms are typically faster and simpler than full scalar multiplication algorithms-the x-only Montgomery ladder being a case in point for elliptic curves-and these algorithms have therefore become the basis of the fastest and safest Diffie-Hellman implementations (such as Bernstein's Curve25519 software [START_REF] Daniel | Curve25519: New Diffie-Hellman speed records[END_REF] and its heirs). We also see excellent Diffie-Hellman implementations based on Kummer surfaces [START_REF] Daniel | Kummer strikes back: New DH speed records[END_REF].

But it is widely believed that pseudomultiplication cannot be used for general cryptosystems, because G/ ±1 is not a group: since the "sign" of the output of a pseudomultiplication is ambiguous, it cannot be used as the input to individual group operations. This belief has so far disqualified Kummers as candidates for implementing most common signature schemes, as well as encryption schemes such as ElGamal [START_REF] Elgamal | A public key cryptosystem and a signature scheme based on discrete logarithms[END_REF]. As a result, we tend not to see highly competitive genus 2 implementations of signatures and public-key encryption schemes, because we are hamstrung by relatively slow and non-uniform Jacobian arithmetic.

Recovering group elements after pseudomultiplication. The folklore that Kummer surfaces cannot be used in true group-based cryptosystems may appear mathematically true-but it is algorithmically false.

Indeed, it has long been known that x-only arithmetic on elliptic curves can be used for full scalar multiplication: López and Dahab [START_REF] López | Fast multiplication on elliptic curves over GF(2 m ) without precomputation[END_REF] (followed by Okeya and Sakurai [START_REF] Okeya | Efficient elliptic curve cryptosystems from a scalar multiplication algorithm with recovery of the y-coordinate on a Montgomeryform elliptic curve[END_REF] and Brier and Joye [START_REF] Brier | Weierstraß elliptic curves and side-channel attacks[END_REF]) showed that the auxiliary values computed by the x-only Montgomery ladder can be used to recover the missing y-coordinate, and hence to compute full scalar multiplications on elliptic curves.

In this work we extend this technique from elliptic curves to genus 2, and from one-to two-dimensional scalar multiplication. In the abstract, our algorithms follow the same common pattern:

1. First, Project the inputs (and possibly some auxiliary elements) from G to G/ ±1 using the x map; 2. then pseudomultiply in G/ ±1 (that is, compute x([m]P ) or x([m]P ⊕[n]Q) using a differential addition chain); 3. finally, efficiently Recover the correct preimage

[m]P (or [m]P ⊕ [n]Q) in G
from the outputs of the pseudomultiplication.

The one-dimensional version of this pattern in §3 uses the well-known Montgomery ladder [START_REF] Peter | Speeding the Pollard and elliptic curve methods of factorization[END_REF] for its pseudomultiplication; the two-dimensional version in §4 is based on Bernstein's binary differential addition chain [START_REF] Daniel | Differential addition chains[END_REF]. In §5 we apply our algorithms to various models of elliptic curves. We recover Okeya-Sakurai and Brier-Joye multiplication, along with new uniform, compact two-dimensional scalar multiplication algorithms.

Moving to genus 2, the Jacobian point recovery method we present in §6 solves the problem of uniform genus 2 arithmetic (at least for scalar multiplication): rather than wrestling with the special cases of Cantor's algorithm, we can descend to the faster, more uniform Kummer surface, pseudomultiply there, and then recover the right Jacobian point afterwards. Our methods work on the most general form of the Kummer, and are easy to adapt to more specialized models.

In §7 we specialize to the Gaudry model for Kummer surfaces, which have the fastest known arithmetic. The result is exceptionally fast one-and twodimensional scalar multiplication algorithms for the genus 2 Jacobians that admit these Kummer surfaces. Finally, in §8 we give a concrete example of how our algorithms can be applied to create fast instances of Schnorr signature schemes.

Remark 1. Damien Robert has pointed out to us that his recent preprint with David Lubicz [START_REF] Lubicz | Arithmetic on abelian and kummer varieties[END_REF] uses similar techniques to improve the efficiency of their explicit arithmetic of general abelian varieties in arbitrarily high dimension based on theta functions. See Remark 6 below for further details.

Notation and conventions Throughout, F q denotes a finite field of characteristic > 3. As usual, we use M, S, I and a to respectively denote the costs of one field multiplication, squaring, inversion, and addition in F q . We also use m c to denote the cost of a field multiplication by a fixed constant c, which is often significantly different to M, i.e., when c is a small curve constant.

In our algorithms, we use the notation (x 1 , . . . , x n ) ← (y 1 , . . . , y n ) to denote parallel assignment. In sequential terms, (x, y) ← (z, w) is equivalent to x ← z; y ← w (and to y ← w; x ← z), while (x, y) ← (y, x), which swaps x and y, is equivalent to t ← x; x ← y; y ← t, where t is a temporary variable.

Key subroutines

In this paper we work with various models of elliptic curves, Jacobians, and Kummer surfaces, but viewed from a high level the algorithms are essentially the same; we can therefore make some substantial simplifications by presenting them as template algorithms acting on an abstract abelian group G, and its quotient G/ ± 1, by black-box subroutines (whose implementation details we will provide later). We therefore assume the existence of six algorithms acting on elements of G and G/ ±1 : 1. Project : G → G/ ±1 implements the mapping x : G → G/ ±1 ; that is, Project(P ) = x(P ) for all P ∈ G . This is trivial in the elliptic context, where it amounts to dropping one of the coordinates, and only slightly less straightforward in genus 2.

2. xDBL : G/ ±1 → G/ ±1 implements pseudo-doubling in G/ ±1 : that is, xDBL(x(P )) = x([2]P ) for all P ∈ G .

Efficient formulae for xDBL are known in each of our contexts. 3. xADD : (G/ ±1 ) 3 → G/ ±1 implements the standard differential addition:

xADD(x(P ), x(Q), x(P ⊖ Q)) = x(P ⊕ Q) for all P = Q ∈ G .
As with xDBL, efficient formulae for xADD are known in each of our contexts. 4. xDBLADD : (G/ ±1 ) 3 → (G/ ±1 ) 2 implements a simultaneous doubling and differential addition:

xDBLADD(x(P ), x(Q), x(Q ⊖ P )) = (x([2]P ), x(P ⊕ Q)) for all P = Q ∈ G .
While xDBLADD may be naïvely defined by computing an xDBL and an xADD separately-that is, as

xDBLADD(x(P ), x(Q), x(Q ⊖ P )) = (xDBL(x(P )), xADD(x(P ), x(Q), x(Q ⊖ P )))
-sometimes in practice it can be implemented more efficiently by exploiting shared intermediate operands, so we treat it as a distinct operation. 5. ADD : G × G → G computes the group law in G:

ADD(P, Q) = P ⊕ Q for all P, Q ∈ G .
ADD is only used once in the two-dimensional algorithm; it is not used at all in the one-dimensional algorithm. Minimizing ADDs is a deliberate strategy: in practice, it is often a relatively costly and potentially non-constant-time operation compared with xADD (especially in genus 2). 6. Recover : G × (G/ ±1 ) 2 → G/ ±1 computes preimages under x:

Recover(P, x(Q), x(Q ⊕ P )) = Q for all P, Q ∈ G \ G[2] .
This operation was introduced for binary elliptic curves by López and Dahab [START_REF] López | Fast multiplication on elliptic curves over GF(2 m ) without precomputation[END_REF], for Montgomery models of elliptic curves by Okeya and Sakurai [START_REF] Okeya | Efficient elliptic curve cryptosystems from a scalar multiplication algorithm with recovery of the y-coordinate on a Montgomeryform elliptic curve[END_REF], and for short Weierstrass models of elliptic curves by Brier and Joye [START_REF] Brier | Weierstraß elliptic curves and side-channel attacks[END_REF]. We extend these techniques to genus 2 in §6 and §7.

Uniform one-dimensional scalar multiplication

Algorithm 1 is a template for uniform one-dimensional scalar multiplication; it is suitable for use anywhere in curve-based cryptosystems where the calculation (m, P ) → [m]P is required, but especially those where m is secret. Algorithm 1 applies the Project-pseudomultiply-Recover pattern to lift the x-only Montgomery ladder for pseudomultiplication in G/ ±1 to a full scalar multiplication routine for G, generalizing and abstracting the methods of [START_REF] López | Fast multiplication on elliptic curves over GF(2 m ) without precomputation[END_REF], [START_REF] Okeya | Efficient elliptic curve cryptosystems from a scalar multiplication algorithm with recovery of the y-coordinate on a Montgomeryform elliptic curve[END_REF], and [START_REF] Brier | Weierstraß elliptic curves and side-channel attacks[END_REF]. Proof. Line 1 of Algorithm 1 calls Project to map P into G/ ±1 . Lines 2-9 are just the standard Montgomery ladder [START_REF] Peter | Speeding the Pollard and elliptic curve methods of factorization[END_REF]. At the end of each of the β iterations of the loop (each of which calls xDBLADD once), we have

Input : A positive integer m = β-1 i=0 mi2 i , with m β-1 = 0,
(t 1 , t 2 ) = (x([⌊m/2 i ⌋]P ), x([⌊m/2 i ⌋ + 1]P )) ;
so at the end of the loop, at Line 10, (t 1 , t 2 ) = (x([m]P ), x([m]P ⊕ P )). The Recover(P, t 1 , t 2 ) in Line 10 therefore yields [m]P .

⊓ ⊔

In its abstract form, Algorithm 1 is uniform and constant-time with respect to fixed-length m. In practice, the implementation of xDBLADD must also be uniform and constant-time. If uniform, constant-time behaviour is required with respect to P , then the implementation of Project must also be uniform and constant-time.

Uniform two-dimensional scalar multiplication

Algorithm 3 is a template for uniform two-dimensional scalar multiplication. It is intended for use in cryptographic routines that require computing [m]P ⊕[b]Q, especially when at least one of m and n are secret, but it is also useful for implementing endomorphism-accelerated scalar multiplications, which compute [m]P as [m 0 ]P ⊕ [m 1 ]φ(P ). Algorithm 3 applies the Project-pseudomultiply-Recover pattern to a pseudomultiplication based on Bernstein's binary differential addition chain [START_REF] Daniel | Differential addition chains[END_REF]. Our algorithm is similar to its x-only counterpart, which was used (without any proof of correctness) in [START_REF] Costello | Faster compact Diffie-Hellman: Endomorphisms on the x-line[END_REF].

Bernstein's binary differential addition chain

Bernstein defined his binary differential addition chain in [2, §4] as follows. First, set C 0 (0, 0) = C 1 (0, 0) := ((0, 0), (1, 0), (0, 1), (1, -1)) ; then the chain C D (A, B) is defined recursively by

C D (A, B) := C d (a, b) || (O, E, M )
where || denotes concatenation, O, E, and M are defined by

O := (A + (A + 1 mod 2), B + (B + 1 mod 2)) , (1)
E := (A + (A + 0 mod 2), B + (B + 0 mod 2)) , (2)
M := (A + (A + D mod 2), B + (B + D + 1 mod 2)) , (3)
and a, b, and d are defined by

a := ⌊A/2⌋ , b := ⌊B/2⌋ , d :=          D if a ≡ A and b ≡ B (mod 2) , 0 if a ≡ A and b ≡ B (mod 2) , 1 if a ≡ A and b ≡ B (mod 2) , 1 -D if a ≡ A and b ≡ B (mod 2) .
Observe that O contains two odd integers, E two even integers, and M is "mixed", with one even and one odd integer. The differences M -O, M -E, and O -E depend only on D and the parities of A and B, as shown in Table 1.

Table 1. The differences between M , O, and E as functions of D and A, B (mod 2).

A (mod 2) B (mod 2) O -E M -O M -E 0 0 (1, 1) (D -1, -D) (D, 1 -D) 0 1 (1, -1) (D -1, D) (D, D -1) 1 0 (-1, 1) (1 -D, -D) (-D, 1 -D) 1 1 (-1, -1) (1 -D, D) (-D, D -1)
By definition, the triple (O, E, M ) contains three of the four pairs (A, B), (A + 1, B), (A, B + 1), and (A + 1, B + 1). The missing pair is (A + (A + D + 1 mod 2), B + (B + D mod 2)), from which it follows immediately that ,o,e,m,O,E,M ). Then O, E, and M can be expressed in terms of o, e, m, D, and the parities of A, B, a, and b using the relations in Table 2.

C (A mod 2) (A, B) contains (A, B) .
Proof. The result follows-after elementary but lengthy calculations-from the definition of C D (A, B), considerations of parity, and the values in Table 1 applied to o, e, and m, with d derived from the first four columns.

⊓ ⊔ Table 2. Relations between the adjacent triples (o, e, m) and (O, E, M ).

A -B A -a B -b D O difference E M difference (mod 2) (mod 2) (mod 2) of summands of summands 0 0 0 0 o + e ±(1, 1) 2e m + e ±(0, 1) 0 0 0 1 o + e ±(1, 1) 2e m + e ±(1, 0) 0 0 1 0 o + e ±(1, 1) 2m m + e ±(0, 1) 0 0 1 1 o + e ±(1, 1) 2m m + o ±(1, 0) 0 1 0 0 o + e ±(1, 1) 2m m + o ±(0, 1) 0 1 0 1 o + e ±(1, 1) 2m m + e ±(1, 0) 0 1 1 0 o + e ±(1, 1) 2o m + o ±(0, 1) 0 1 1 1 o + e ±(1, 1) 2o m + o ±(1, 0) 1 0 0 0 o + e ±(1, -1) 2e m + e ±(0, 1) 1 0 0 1 o + e ±(1, -1) 2e m + e ±(1, 0) 1 0 1 0 o + e ±(1, -1) 2m m + e ±(0, 1) 1 0 1 1 o + e ±(1, -1) 2m m + o ±(1, 0) 1 1 0 0 o + e ±(1, -1) 2m m + o ±(0, 1) 1 1 0 1 o + e ±(1, -1) 2m m + e ±(1, 0) 1 1 1 0 o + e ±(1, -1) 2o m + o ±(0, 1) 1 1 1 1 o + e ±(1, -1) 2o m + o ±(1, 0)

Encoding the chain

Our aim is to use

C m0 (m, n) to compute x([m]P ⊕ [n]Q), where m = β-1 i=0 m i 2 i and n = β-1 i=0 n i 2 i .
Treating the m i and n i as bits, with ⊕ denoting binary addition (xor) and ⊗ binary multiplication (and), we define the sequence of transition vectors

c i := ((m i ⊕ n i ), (m i ⊕ m i+1 ), (n i ⊕ n i+1 ), d i ) for 0 ≤ i < β -1 ,
where d 0 := m 0 and

d i+1 := ((d i ⊕ 1) ⊗ (m i ⊕ m i+1 )) ⊕ (d i ⊗ (n i ⊕ n i+1 ⊕ 1)) for i ≥ 0 .
The coordinates of c i correspond to the first four columns of Table 2, with

A = ⌊m/2 i ⌋, a = ⌊m/2 (i+1) ⌋, B = ⌊n/2 i ⌋, b = ⌊n/2 (i+1) ⌋, D = d i , and d = d i+1 .
In view of Lemma 2, given m and n, we can iteratively reconstruct the entire chain C m0 (m, n) from d β-1 and the sequence of transition vectors c 0 , . . . , c β-2 . Algorithm 2 computes precisely this data to encode C m0 (m, n). We note that the c i also encode the difference elements required for each differential addition.

Two-dimensional scalar multiplication

If we map each pair (a, b) in C m0 (m, n) to the element x([a]P ⊕ [b]Q) in G/ ±1 , then C m0 (m, n) yields a method of computing x([m]P ⊕[n]Q) using a sequence of Input : Positive β-bit integers m = β-1 i=0 mi2 i and n = β-1 i=0 ni2 i . Output: A sequence C of β -1 transition vectors and one additional bit d β-1 ,
encoding the differential addition chain Cm 0 (m, n).

1 C ← [ ] ; 2 d ← m0 ; 3 for i ← 0 up to β -2 do 4 Append ((mi ⊕ mi+1), (ni ⊕ ni+1), (mi+1 ⊕ ni+1), d) to C ; 5 d ← ((d ⊕ 1) ⊗ (mi ⊕ mi+1)) ⊕ (d ⊗ (ni ⊕ ni+1 ⊕ 1)) ; 6 end 7 return (C, d) ;
Algorithm 2: ChainEncoding: Computes and encodes Bernstein's twodimensional "binary" differential addition chain xADDs (and xDBLs) with the fixed differences x(P ), x(Q), x(P ⊕Q), and x(P ⊖Q). Thus, the sequence of transition vectors in our encoding of C m0 (m, n) gives us a natural iterative algorithm for computing x([m]P ⊕ [n]Q) starting from the fixed differences. Using this as the core of our Project-pseudomultiply-Recover

pattern yields Algorithm 3, which computes [m]P ⊕ [n]Q in G.
Theorem 1. Let P and Q be elements of G, let m and n be positive integers, and let β be the bitlength of max(m, n). Algorithm 3 computes [m]P ⊕ [n]Q using one call to ADD, three calls to Project, β calls to each of xADD and xDBLADD, and one call to Recover.

Proof. Algorithm 3 begins by calling Algorithm 2 to compute the sequence of transition vectors and d β-1 for the given m and n. Lines 2-4 compute the required differences x P := x(P ), x Q := x(Q), x ⊕ := x(P ⊕Q), and x ⊖ := x(P ⊖Q) by computing S := P ⊕ Q with the single call to ADD, then using the three calls to Project to compute x P , x Q , and

x ⊕ = x(S), before x ⊖ = xADD(x P , x Q , x ⊕ ).
Throughout the rest of the algorithm O corresponds to the odd-odd pair, E to the even-even pair, and M to the mixed pair in a segment of C m0 (m, n), where each pair (a, b) is associated with x([a]P ⊕ [b]Q). Lines 5-12 use a single xDBLADD to initialize O, E, and M such that

C d β-1 (m β-1 , n β-1 ) = C 0 (0, 0)||(O, E, M ).
Lines 13-31 iterate over the sequence of transition vectors to compute the last three elements of C m0 (m, n). After each iteration, the triple (O, E, M ) satisfies

C di (⌊m/2 i ⌋, ⌊n/2 i ⌋) = C di+1 (⌊m/2 i+1 ⌋, ⌊n/2 i+1 ⌋) || (O, E, M ) .
Each iteration requires two differential additions and a pseudodoubling, but we observe that the pseudodoubling always applies to one of the arguments of a differential addition, so each of the β -1 iterations requires exactly one xADD and one xDBLADD.

After the loop is completed, at Line 33, Eqs. ( 1), (2), and (3) imply

O = x(R ⊕ ∆ O ) , E = x(R ⊕ ∆ E ) , and M = x(R ⊕ ∆ M ) , where R = [m]P ⊕ [n]Q and ∆ O , ∆ E
, and ∆ M are given in the following table:

Input : Positive β-bit integers m = β-1 i=0 mi2 i and n = β-1 i=0 ni2 i with m β-1 and n β-1 not both zero, and elements

P, Q of G Output: R = [m]P ⊕ [n]Q ((c0, . . . , c β-2 ), d β-1 ) ← ChainEncoding(m,n) ; S ← ADD(P, Q) ; (xP , xQ, x⊕) ← (Project(P ), Project(Q), Project(S)) ; x⊖ ← xADD(xP , xQ, x⊕) ; switch (m β-1 , n β-1 , d β-1 ) do Algorithm 3: Two-dimensional uniform scalar multiplication template (m 0 , n 0 ) ∆ O ∆ E ∆ M if d 0 = 0 ∆ M if d 0 = 1 [m]P ⊕ [n]Q (0, 0) P ⊕ Q 0 G Q P Recover(P ⊕ Q, E, O) (0, 1) P Q 0 G P ⊕ Q Recover(P, M, O) (1, 0) Q P P ⊕ Q 0 G Recover(P, M, E) (1, 1) 0 G P ⊕ Q P P ⊕ Q Recover(P ⊕ Q, O, E)
In each case, we can recover [m]P ⊕[n]Q by applying Recover with the arguments specified by the last column of the corresponding row. This is precisely what is done in Lines 33-38; Line 39 then returns the result, [m]P ⊕ [n]Q.

⊓ ⊔

Like Algorithm 1, Algorithm 3 is uniform and constant-time in its abstract form, at least for fixed-length multiscalars (m, n). In practice, for Algorithm 3 to be uniform and constant-time with respect to m and n, the implementations of xADD and xDBLADD must be uniform and constant-time. For uniform and constant-time behaviour with respect to P and Q, the implementations of ADD and Project must also be uniform and constant-time.

Remark 2. The core of Algorithm 3 is similar to the algorithm in [14, App. C], but with a different (and slightly simpler) encoding of the addition chain. Here, the i-th transition vector is ((

m i ⊕n i ), (m i ⊕m i+1 ), (n i ⊕n i+1 ), d i ); in [14, App. C] it is ((m i ⊕ m i+1 ⊕ n i ⊕ n i+1 ), (m i ⊕ m i+1 ), (m i ⊕ n i ), d i )
, and the order of the vectors is reversed (as is the order of the loop iterations).

Remark 3. Implementers should notice that every adjacent pair of xADD and xDBLADD operations in Algorithm 3 share one operand in their differential additions. Further savings might therefore be made by sharing a few intermediate calculations between the xADD and xDBLADD calls-that is, by implicitly defining an xDBLADDADD operation, as in [14, App. C]. We have not done this here, for two reasons. First, it somewhat obscures the fundamental symmetry of the addition chain. Second, looking ahead at the explicit formulae for xADD and xDBLADD suggests that there are no intermediate calculations that can be shared in the elliptic curve scenarios. However, close inspection of the fast Kummer arithmetic in §7.3 reveals that two adjacent differential additions sharing a common summand can be merged to save 8 field additions. It is unclear to us that whether this potential saving would give a net benefit after the code complexity required to perform the merging, so we have not included this potential saving in our operation counts.

Remark 4. These techniques should readily extend to the higher-dimensional Montgomery-like differential addition chains described by Brown [START_REF] Daniel | Multi-dimensional montgomery ladders for elliptic curves[END_REF]. We do not investigate this here.

Efficient scalar multiplication on elliptic curves

We now pass from the abstract to the concrete. Applying the one-dimensional pattern to elliptic curves, we recover Okeya-Sakurai-style multiplication algorithms; applying the two-dimensional pattern yields something new.

Montgomery models

The most obvious application of our methods is to elliptic curves with Montgomery models

BY 2 Z = X(X 2 + AXZ + Z 2 ) ,
which are defined to optimize the xADD, xDBL, and xDBLADD operations. Indeed, historically, the first appearance of the Project-pseudomultiply-Recover pattern was in the context of one-dimensional scalar multiplication on Montgomery models [START_REF] Okeya | Efficient elliptic curve cryptosystems from a scalar multiplication algorithm with recovery of the y-coordinate on a Montgomeryform elliptic curve[END_REF]. Important examples of Montgomery curves include Curve25519, recently recommended for standardization by the CFRG.

In this context, Project : (x, y) → x is trivially computed. The ADD, xADD, xDBL, and xDBLADD operations are all presented-and thoroughly costed-in the EFD [START_REF] Daniel | Explicit formulas database[END_REF]. The Recover operation was derived by Okeya and Sakurai in [31, §3]. The operation counts for all of these operations are compiled in Table 3. We now apply the Project-pseudomultiply-Recover pattern using the routines above. In the one-dimensional case, we recover the Okeya-Sakurai algorithm [START_REF] Okeya | Efficient elliptic curve cryptosystems from a scalar multiplication algorithm with recovery of the y-coordinate on a Montgomeryform elliptic curve[END_REF].

Theorem 2. Let E/F q be an elliptic curve in Montgomery form.

1. Let P be a point in E(F q ) \ E [START_REF] Daniel | Differential addition chains[END_REF], and let m be a positive β-bit integer. Then Algorithm 1 computes [m]P in

(5β + 10)M + (4β + 2)S + (β + 1)m A + 1m B + (8β + 6)a + 1I .
2. Let P and Q be points in E(F q ) \ E [START_REF] Daniel | Differential addition chains[END_REF], and let m and n be positive β-bit integers. Then Algorithm 3 computes

[m]P ⊕ [n]Q in (8β + 14)M + (6β + 2)S + (β + 1)m A + 1m B + (14β + 13)a + 2I .

Edwards models

As a byproduct of Theorem 2, we obtain new scalar multiplication algorithms for Edwards models [START_REF] Edwards | A normal form for elliptic curves[END_REF]. Indeed, every twisted Edwards model can be transformed into a Montgomery model with a linear change of coordinates, so one option to perform scalar multiplications on Edwards curves is to pass to and from an associated Montgomery model, making use of the operations summarized in Table 3. Important examples of Edwards models in practice include Ed25519 [START_REF] Daniel | High-speed high-security signatures[END_REF] and Goldilocks [START_REF] Hamburg | Ed448-Goldilocks, a new elliptic curve[END_REF].

Another option here is to exploit pseudomultiplication formulae that are native to Edwards curves, as described by Gaudry and Lubicz [20,§6.2]. In this case, for Edwards curves E/F q : x 2 + y 2 = 1 + dx 2 y 2 with d = r 2 and r ∈ F q , we define Project : (x, y) → y and use the formulae in the EFD [START_REF] Daniel | Explicit formulas database[END_REF] to compute the pseudomultiplication via Algorithm 1.

The Recover:

(P, y(Q), y(Q ⊕ P )) → Q = (x(Q), y(Q))
operation is defined by rearranging the Edwards addition law between P and Q, to make x(Q) the subject, i.e.,

x(Q) = y(Q ⊕ P ) -y(Q)y(P ) x(P ) (dy(P )y(Q)y(Q ⊕ P ) -1)
.

Scalar multiplication on short Weierstrass models

Not every elliptic curve has a Montgomery or Edwards model; the most general form for an elliptic curve over F q in characteristic greater than 3 is the short Weierstrass model E :

y 2 = x 3 + Ax + B ⊂ P 2 .
Important examples of curves commonly implemented as Weierstrass models include the NIST [START_REF]Recommended elliptic curves for federal government use[END_REF] and Brainpool curves [START_REF] Brainpool | ECC Brainpool standard curves and curve generation v. 1.0[END_REF].

Our key subroutines are implemented as follows. For brevity, we use affine coordinates here, but the resulting formulae are easily projectivized: see Brier and Joye or the EFD. In this context, Project : P = (x, y) → x is trivially computed. The ADD, xADD, and xDBL, and xDBLADD operations are all specified as efficient straight-line programs in the EFD [START_REF] Daniel | Explicit formulas database[END_REF]. Brier and Joye describe Recover for short Weierstrass models in [START_REF] Brier | Weierstraß elliptic curves and side-channel attacks[END_REF]Prop. 3

]: it maps (P, x(Q), x(P ⊕ Q)) to Q = (x(Q), y(Q)), where y(Q) = 2B + (A + x(P )x(Q))(x(P ) + x(Q)) -x(P ⊕ Q)(x(P ) -x(Q)) 2 2y(P ) .
Table 4 summarizes the operation counts for all of this routines.

Applying Algorithm 1 with these subroutines yields the scalar multiplication algorithm for Weierstrass models in [START_REF] Brier | Weierstraß elliptic curves and side-channel attacks[END_REF]. Applying Algorithm 3 with these subroutines yields a new algorithm for two-dimensional scalar multiplication on Weierstrass curves. Theorem 3. Let E/F q be an elliptic curve in short Weierstrass form. Proof. Take the values from Table 4 in Lemma 1 (for the first part) and Theorem 1 (for the second). ⊓ ⊔

Genus 2 Jacobians and General Kummer Surfaces

We now turn our focus to genus 2 cryptosystems, where G is (a subgroup of) the Jacobian J C of a genus 2 curve

C : y 2 = f (x) = f 6 x 6 + f 5 x 5 + • • • + f 1 x + f 0 over F q .
Elements of J C (F q ) are presented in their standard Mumford representation:

P ∈ J C (F q ) ←→ a(x) = x 2 + a 1 x + a 0 , b(x) = b 1 x + b 0
where a 1 , a 0 , b 1 , and b 0 are in F q and b(x

) 2 ≡ f (x) (mod a(x)) .
The quotient G/ ±1 is (a subset of) a Kummer surface, which is a quartic surface in P 3 . There are two main ways to represent the Kummer: the "general" model K gen

C

(see Cassels and Flynn [11,Ch. 3]) and the "fast" model K fast C algorithmically developed by Chudnovsky and Chudnovsky [START_REF] Chudnovsky | Sequences of numbers generated by addition in formal groups and new primality and factorization tests[END_REF], and introduced in cryptography by Gaudry [START_REF] Gaudry | Fast genus 2 arithmetic based on Theta functions[END_REF].

Broadly speaking, the fast Kummer model corresponds to the Montgomery model for elliptic curves, while the general model corresponds to the Weierstrass model. Indeed, as in the elliptic curve situation, fast Kummers offer significant gains in performance and uniformity (indeed, they are at the heart of a number of record-breaking Diffie-Hellman implementations), but at the price of a lot of rational 2-torsion: hence, not every Kummer can be put in fast form.

We will define the Project and Recover operations for general Kummers; since any Kummer surface can be transformed into a general model over the ground field, these Project and Recover routines can easily be specialized to fast Kummers (or to any other interesting models of Kummer surfaces). We do not lose much in taking this approach compared with deriving a new, specialized Project and Recover for Jacobians with fast Kummers from scratch, because Algorithms 1 and 3 only call Project a few times, and Recover once. Scalar pseudomultiplication on general Kummers is relatively slow, and has held little interest for cryptographers thus far; we briefly discuss their arithmetic and pseudomultiplication in §6.5.

Genus 2 arithmetic and side-channel attacks

The group law on J C (and hence the ADD operation) is typically computed using Cantor's algorithm, specialized to genus 2. This style of arithmetic has a serious drawback in cryptographic contexts: it is highly susceptibile to simple sidechannel attacks. Similar to the textbook addition for short Weierstrass models of elliptic curves, Cantor's algorithm in genus 2 treats several input cases differently, branching off into distinct explicit computations. Explicit formulae that are derived for generic additions fail to compute correctly when one or both of the inputs are special points-that is, points in J C where a(x) = x -α or a(x) = (x -α) 2 . While such special points are sparse enough in J C that random scalar multiplications do not encounter them, they are plentiful enough that attackers could easily mount exceptional procedure attacks [START_REF] Izu | Exceptional procedure attack on elliptic curve cryptosystems[END_REF], forcing legitimate users into special cases and using timing variabilities to recover secret data.

The Project-pseudomultiply-Recover approach detailed in §6.2 addresses this problem. Like x-only (pseudo)scalar multiplication on elliptic curves, the explicit formulae for doublings and differential additions on certain Kummer surfaces associated to J C are well behaved on all inputs, including the images of special points in J C under the Project map. The only time special points might be encountered is at the two end points of the scalar multiplication; but suitable point validation can detect and reject special points as inputs, and special outputs in J C are not the goal of an exceptional procedure attack-adversaries can only hope to trigger exceptional points early on in a scalar multiplication, where the number of possible intermediate points is subexponential.

Point recovery on the general Kummer

The general model of the Kummer surface K gen C associated to J C is defined by

K gen C : K 2 (ξ 1 , ξ 2 , ξ 3 )ξ 2 4 + K 1 (ξ 1 , ξ 2 , ξ 3 )ξ 4 + K 0 (ξ 1 , ξ 2 , ξ 3 ) = 0 , (4)
where

K 2 = ξ 2 2 -4ξ 1 ξ 3 , K 1 = -2(2f 0 ξ 3 1 + f 1 ξ 2 1 ξ 2 + 2f 2 ξ 2 1 ξ 3 + f 3 ξ 1 ξ 2 ξ 3 + 2f 4 ξ 1 ξ 2 3 + f 5 ξ 2 ξ 2 3 + 2f 6 ξ 3 3 ), K 0 = (f 2 1 -4f 0 f 2 )ξ 4 1 -4f 0 f 3 ξ 3 1 ξ 2 -2f 1 f 3 ξ 3 1 ξ 3 -4f 0 f 4 ξ 2 1 ξ 2 2 + 4(f 0 f 5 -f 1 f 4 )ξ 2 1 ξ 2 ξ 3 + (f 2 3 + 2f 1 f 5 -4f 2 f 4 -4f 0 f 6 )ξ 2 1 ξ 2 3 -4f 0 f 5 ξ 1 ξ 3 2 + 4(2f 0 f 6 -f 1 f 5 )ξ 1 ξ 2 2 ξ 3 + 4(f 1 f 6 -f 2 f 5 )ξ 1 ξ 2 ξ 2 3 -2f 3 f 5 ξ 1 ξ 3 3 -4f 0 f 6 ξ 4 2 -4f 1 f 6 ξ 3 2 ξ 3 -4f 2 f 6 ξ 2 2 ξ 2 3 -4f 3 f 6 ξ 2 ξ 3 3 + (f 2 5 -4f 4 f 6 )ξ 4 3 .

Projection from Jacobians to general Kummers

Project implements the map

J C → K gen C described in [11, Eqs. (3.1.3-5)]; it maps generic points x 2 + a 1 x + a 0 , b 1 x + b 0 in J C to (ξ 1 : ξ 2 : ξ 3 : ξ 4 ) in K gen C
, where

ξ 1 = 1 , ξ 2 = -a 1 , ξ 3 = a 0 , and (5)
ξ 4 = b 2 1 + (a 2 1 -a 0 ) f 5 a 1 -f 6 (a 2 1 -a 0 ) + a 1 (f 3 -f 4 a 1 ) -f 2 . ( 6
)
This costs 5M + 2S + 7a (saving 1M if f 6 = 0).

Recovering Jacobian points from general Kummers

We now derive explicit formulae for

Recover : (P, x(Q), x(Q ⊕ P )) -→ Q in genus 2.
We follow the approach of Okeya-Sakurai [START_REF] Okeya | Efficient elliptic curve cryptosystems from a scalar multiplication algorithm with recovery of the y-coordinate on a Montgomeryform elliptic curve[END_REF] and Brier-Joye [START_REF] Brier | Weierstraß elliptic curves and side-channel attacks[END_REF] for elliptic curves, rewriting the equations used for computing x(Q ⊕ P ) in terms of P = (x(P ), y(P )) and Q = (x(Q), y(Q)), making y(Q) the unknown. We may suppose that P and Q are nonzero and not points of order 2. Suppose we are given4

P = x 2 + a 1 (P )x + a 0 (P ), b 1 (P )x + b 0 (P ) ∈ J C (F q ) , x(Q) = (1 : ξ 2 : ξ 3 : ξ 4 ) ∈ K gen C (F q ) , x(Q ⊕ P ) = (1 : ξ ⊕ 2 : ξ ⊕ 3 : ξ ⊕ 4 ) ∈ K gen C (F q ) ;
we want to Recover

Q = x 2 + a 1 (Q)x + a 0 (Q), b 1 (Q)x + b 0 (Q) ∈ J C (F q ) .
We already have a 1 (Q) = -ξ 2 and a 0 (Q) = ξ 3 from Eq. ( 5); it remains to compute b 1 (Q) and b 0 (Q).

Okeya and Sakurai noticed that the formulae for y-coordinate recovery on Montgomery curves are simpler if x(Q ⊖ P ) is also known (see [31, pp. 129-130]); we observed the same simplification in genus 2, where it was evident that computing x(Q ⊖ P ) from x(Q), x(P ) and x(Q ⊕ P ) (using one more differential addition on the corresponding Kummer) yielded a faster overall recovery. We therefore begin our Recover with a call to xADD, to compute

x(Q ⊖ P ) = (1 : ξ ⊖ 2 : ξ ⊖ 3 : ξ ⊖ 4 ) = xADD(Q, P, x(Q ⊕ P )) .
Since P and Q are nonzero, they correspond to unique degree-2 divisors on C (cf. [START_REF] Cassels | Prolegomena to a middlebrow arithmetic of curves of genus 2[END_REF]Ch. 1]):

P ←→ [(u P , v P ) + (u ′ P , v ′ P )] , Q ←→ [(u Q , v Q ) + (u ′ Q , v ′ Q )] .
We do not compute the values of u P , v P , u ′ P , v ′ P , u Q , v Q , u ′ Q , and v ′ Q (which are generally defined over a quadratic extension): here they simply serve as formal devices, to aid our derivation of recovery formulae. Let

G 1 := E 1 + E 2 , G 2 := u ′ P E 1 + u P E 2 , G 3 := E 3 + E 4 , G 4 := u ′ Q E 3 + u Q E 4 ,
where

E 1 = v P (u P -u ′ P )(u P -u Q )(u P -u ′ Q ) , E 2 = v ′ P (u ′ P -u P )(u ′ P -u Q )(u ′ P -u ′ Q ) , E 3 = v Q (u Q -u ′ P )(u Q -u P )(u Q -u ′ Q ) , E 4 = v ′ Q (u ′ Q -u P )(u ′ Q -u ′ P )(u ′ Q -u Q )
.

The G i are functions of P and Q, because they are symmetric with respect to

(u P , v P ) ↔ (u ′ P , v ′ P ) and (u Q , v Q ) ↔ (u ′ Q , v ′ Q )
; below, we will compute them as functions of P , x(Q ⊕ P ), and x(Q ⊖ P ). For notational convenience, we define

Z 1 := -(ξ 2 + a 1 (P )) , Z 2 := ξ 3 -a 0 (P ) , (7)
Z 3 := a 1 (P )ξ 3 + a 0 (P )ξ 2 , Z 4 := ξ 3 Z 2 -ξ 2 Z 3 , (8)
D := Z 2 2 -Z 1 Z 3 , (9)
∆ := -4G 2 2 + 2(ξ ⊕ 2 + ξ ⊖ 2 )G 1 G 2 -2(ξ ⊕ 3 + ξ ⊖ 3 )G 2 1 , (10)
and

T := f 6 -G 2 1 -G 2 3 = f 6 -G 2 1 - 1 D 2 ξ 4 D + f 0 Z 2 1 -f 1 Z 1 Z 2 + f 2 Z 2 2 -f 3 Z 2 Z 3 +f 4 Z 2 3 -f 5 Z 3 Z 4 + f 6 Z 2 4 . (11)
We now use the fact that the cubic polynomial

l(x) = E 1 (x -u ′ P )(x -u Q )(x -u ′ Q ) + E 2 (x -u P )(x -u Q )(x -u ′ Q ) + E 3 (x -u P )(x -u ′ P )(x -u ′ Q ) + E 4 (x -u P )(x -u ′ P )(x -u Q ) = (G 1 x -G 2 )(x 2 + a 1 (Q)x + a 0 (Q)) + (G 3 x -G 4 )(x 2 -ξ 2 x + ξ 3 )
satisfies ℓ(x) ≡ b(x) mod a(x) when a(x), b(x) is the Mumford representation of P , Q or ⊖(P ⊕ Q) (this is just the geometric definition of the group law on J C ; the cubic ℓ(x) is analogous to the line through P , Q, and ⊖(P ⊕ Q) in the classic elliptic curve group law). Together with b(x) 2 ≡ f (x) (mod a(x)), which is satisfied by the Mumford representation a(x), b(x) of every point on J C , this gives the relations

b 1 (Q) b 0 (Q) = ξ 2 Z 1 + Z 2 Z 1 ξ 3 Z 1 Z 2 G 3 G 4 , (12)
G 3 G 4 = T ∆ ξ ⊖ 3 -ξ ⊕ 3 ξ ⊕ 2 -ξ ⊖ 2 ξ ⊕ 2 ξ ⊖ 3 -ξ ⊖ 2 ξ ⊕ 3 ξ ⊕ 3 -ξ ⊖ 3 G 1 G 2 , (13)
and

G 1 G 2 = 1 D Z 2 -Z 1 -a 0 (P )Z 1 a 1 (P )Z 1 -Z 2 b 1 (P ) b 0 (P ) . (14)
Our Recover operation is now defined as follows: first, we use Eqs. ( 7) and (8) to compute the Z i ; then Eq. ( 9) yields D. Then we can use Eq. ( 14) to obtain G 1 and G 2 , which we use to compute ∆ and T using Eqs. ( 10) and [START_REF] Cassels | Prolegomena to a middlebrow arithmetic of curves of genus 2[END_REF]. Equation (13) then yields G 3 and G 4 , which we substitute into Eq. ( 12) to compute b 1 (Q) and b 0 (Q). We have thus recovered the full point

Q = x 2 + a 1 (Q)x + a 0 (Q), b 1 (Q)x + b 0 (Q) in J C .
Altogether, computing the map (P, x(Q), x(Q ⊕ P ), x(Q ⊖ P )) → Q costs 71M + 8S + 8m c + 35a + 1I; adding the cost of computing x(Q ⊖ P ) via an xADD operation, we obtain the full cost of computing Recover(P, x(Q), x(Q ⊕ P )).

Remark 5. It is worth noting that ξ ⊕ 4 and ξ ⊖ 4 do not appear in our recovery procedure; this may be useful in scenarios where it is advantageous to omit their computation or transmission.

Scalar multiplication on Jacobians via general Kummers

The use of general Kummers in cryptography was investigated by Smart and Siksek [START_REF] Smart | A fast Diffie-Hellman protocol in genus 2[END_REF] and Duquesne [START_REF] Sylvain Duquesne | Montgomery scalar multiplication for genus 2 curves[END_REF]. While they represent a natural generalization of x-only arithmetic for elliptic curves, in their full generality they do not offer competitive performance. The xADD and xDBL operations are defined by complicated biquadratic forms in the ξ i (see [START_REF] Cassels | Prolegomena to a middlebrow arithmetic of curves of genus 2[END_REF]), which are hard to evaluate quickly for general curve parameters. While these formulae are completely compatible with our Project-pseudomultiply-Recover pattern, and yield scalar multiplication algorithms that may be useful to number theorists, the cryptographic applications of these algorithms are à priori limited. We will therefore skip any investigation of these scalar multiplication algorithms here, and move on directly to Gaudry's fast Kummer surfaces.

We remark, nevertheless, that the use of these biquadratic forms in conjunction with the Project and Recover defined above yields a solution to the problem of defining uniform and constant-time scalar multiplication algorithms for general genus 2 Jacobians. We leave the eventual optimization of these algorithms as an open problem, along with their application to Jacobians with special endomorphism structure, or Jacobians whose Kummers can be put into intermediate forms between the general and fast models.

Efficient Scalar Multiplication via Fast Kummers

While the performance of general Kummer models is somewhat disappointing, Gaudry [START_REF] Gaudry | Fast genus 2 arithmetic based on Theta functions[END_REF] showed if we allow a certain 2-torsion structure on the Jacobian, then an alternative classical model for the Kummer (investigated algorithmically by Chudnovsky and Chudnovsky [START_REF] Chudnovsky | Sequences of numbers generated by addition in formal groups and new primality and factorization tests[END_REF]) yields a dramatic speedup, competitive with-and regularly outperforming-elliptic curve arithmetic. Jacobians equipped with these "fast" Kummers are ideal candidates for our techniques; here we use the model with squared theta coordinates described in [START_REF] Cosset | Applications of theta functions for hyperelliptic curve cryptography[END_REF]Ch. 4]. We set

Construction of fast Kummers

λ = ac bd , µ = ce df , ν = ae bf ,
and define an associated genus 2 curve C in Rosenhain form:

C : y 2 = f (x) = x(x -1)(x -λ)(x -µ)(x -ν) .
The fast Kummer surface for C is the quartic surface K fast [START_REF] Sylvain Duquesne | Montgomery scalar multiplication for genus 2 curves[END_REF] where

C ⊂ P 3 defined by K fast C : (X 2 + Y 2 + Z 2 + T 2 ) -F (XT + Y Z) -G(XZ + Y T ) -H(XY + ZT ) 2 = EXY ZT
E := 4abcd (ABCD/((ad -bc)(ac -bd)(ab -cd))) 2 , F := (a 2 -b 2 -c 2 + d 2 )/(ad -bc) , G := (a 2 -b 2 + c 2 -d 2 )/(ac -bd) , H := (a 2 + b 2 -c 2 -d 2 )/(ab -cd) .
For efficient arithmetic on K fast C , we precompute the theta constants

x 0 := 1 , y 0 := a/b , z 0 := a/c , t 0 := a/d , x ′ 0 := 1 , y ′ 0 := A/B , z ′ 0 := A/C , t ′ 0 := A/D .
The image of the identity element of

J C in K fast C is (a : b : c : d) = (1/x 0 : 1/y 0 : 1/z 0 : 1/d 0 ) ;
we also observe that (A : B :

C : D) = (1/x ′ 0 : 1/y ′ 0 : 1/z ′ 0 : 1/d ′ 0 ).

Projection from Jacobians to fast Kummers

Project implements the map

J C → K fast C from [13, §5.3], which is defined for generic points in J C by x 2 + a 1 x + a 0 , b 1 x + b 0 → (X : Y : Z : T ), where X = a a 0 (µ -a 0 )(λ + a 1 + ν) -b 2 0 , Y = b a 0 (νλ -a 0 )(1 + a 1 + µ) -b 2 0 , Z = c a 0 (ν -a 0 )(λ + a 1 + µ) -b 2 0 , T = d a 0 (µλ -a 0 )(1 + a 1 + ν) -b 2 0 .
This costs 11M + 1S + 3m c + 12a + 1I, assuming the output point is normalized with X = 1.

For special points x -α, β in J C , Project is defined by first adding a point of order 2 in J C , e.g., adding x -λ, 0 to get x 2 -(α + β)x + αβ, β α-λ x -βλ α-λ , then using the above map to K fast C , where the translation is undone by using the coordinate permutation-and-negation that corresponds to translation by xλ, 0 (see [19, §3.4]). Finally, 0 JC (whose Mumford representation is 1, 0 ) maps to (a : b : c : d).

Basic fast Kummer arithmetic

We recall the formulae and operation counts for xDBL and xADD on fast Kummers from [19, §3.2]. If

x(P ) = (X 1 : Y 1 : Z 1 : T 1 ) , x(Q) = (X 2 : Y 2 : Z 2 : T 2 ) , x(P ⊖ Q) = (X ⊖ : Y ⊖ : Z ⊖ : T ⊖ ) ,
then xADD maps (x(P ), x(Q), x(P ⊖ Q)) to (X ⊕ : Y ⊕ : Z ⊕ : T ⊕ ), where

X ⊕ = (X ′ + Y ′ + Z ′ + Y ′ ) 2 /X ⊖ , Y ⊕ = (X ′ + Y ′ -Z ′ -Y ′ ) 2 /Y ⊖ , Z ⊕ = (X ′ -Y ′ + Z ′ -Y ′ ) 2 /Z ⊖ , T ⊕ = (X ′ -Y ′ -Z ′ + Y ′ ) 2 /T ⊖ ,
where

X ′ = x ′ 0 (X 1 + Y 1 + Z 1 + T 1 )(X 2 + Y 2 + Z 2 + T 2 ) , Y ′ = y ′ 0 (X 1 + Y 1 -Z 1 -T 1 )(X 2 + Y 2 -Z 2 -T 2 ) , Z ′ = z ′ 0 (X 1 -Y 1 + Z 1 -T 1 )(X 2 -Y 2 + Z 2 -T 2 ) , T ′ = t ′ 0 (X 1 -Y 1 -Z 1 + T 1 )(X 2 -Y 2 -Z 2 + T 2 ) .
We can compute xADD using 8 squares, 7 products (3 by constants, ignoring the multiplication by x ′ 0 = 1), and 4 divisions (or, projectively, another 10 products); but in our applications, (X ⊖ : Y ⊖ : Z ⊖ : T ⊖ ) is fixed. so we can precompute 1/X ⊖ , 1/Y ⊖ , 1/Z ⊖ , and 1/T ⊖ , and then the 4 divisions become 4 products.

The xDBL operation is defined by exactly the same formulae on setting (X 2 :

Y 2 : Z 2 : T 2 ) = (X 1 : Y 1 : Z 1 : T 1 ) and (X ⊖ : Y ⊖ : Z ⊖ : T ⊖ ) = (a : b : c : d) (so X ′ = x ′ 0 (X 1 + Y 1 + Z 1 + T 1 ) 2
, and so on). The combined xDBLADD operation is outlined in Gaudry [19, §3.3]: we share the computation of x ′ 0 (X 1 + Y 1 + Z 1 + T 1 ) 2 , etc., between the calculation of xDBL and xADD. The resulting operation costs 16 products and 1 square.

Finally, the function ADD : P, Q → x(Q⊕P ) in J C is computed via the formulae in [START_REF] Hisil | Jacobian coordinates on genus 2 curves[END_REF]Eq. (12)] at a cost of 22M + 2S + 1I + 27a.

Recovering Jacobian points from fast Kummers

The Recover operation for fast Kummers uses the Recover for general Kummers (defined in §6.2) as a subroutine. To move between the fast Kummer K fast C and the general Kummer K gen C , we modify the map from J C to K fast

C

given in [13, §5.3] (which was in turn derived from [START_REF] Paul Van Wamelen | Equations for the Jacobian of a hyperelliptic curve[END_REF]). This gives a linear isomorphism

τ : K gen C -→ K fast C , defined on generic points by τ : (ξ 1 : ξ 2 : ξ 3 : ξ 4 ) -→ (X : Y : Z : T ) = (ξ 1 : ξ 2 : ξ 3 : ξ 4 )M (a, b, c, d) t , where M (a, b, c, d) =     µ(λ + ν) νλ(1 + µ) ν(λ + µ) µλ(1 + ν) -µ -νλ -ν -µλ µ + 1 λ + ν ν + 1 λ + µ -1 -1 -1 -1     .
The inverse map τ With x : J C → K, the operation Recover: P, x(Q), x(Q ⊕ P ) → Q is computed as follows. As was done in §6.2, the first step is to compute x(Q ⊖ P ) from x(P ), x(Q) and x(Q ⊕ P ) via an xADD operation (note that x(P ) was already obtained during the Project operation). We then use τ -1 to map all four elements, x(P ), x(Q), x(Q ⊕ P ) and x(Q ⊖ P ) to the corresponding general Kummer K gen C . We can then use the Recover operation on K gen C to output a normalized point Q in Mumford coordinates using 71M + 4m c + 8S + 34a + 1I. Recall from §6.2 that the fourth Kummer coordinate ξ 4 is only needed for the Kummer point corresponding to Q. Thus, the map τ -1 is only computed in full once (costing 15M + 12a), but for the three other points can omit the fourth coordinate (costing 11M + 9a).

In total, including the initial xADD operation (costing 14M+4S+3m c +12a), the map Recover:

P, x(Q), x(Q ⊕ P ) → Q when x : J C → K fast C costs 133M + 12S + 7m c + 97a + 1I.
Here we count multiplications by the curve constants f i as full multiplications.

Scalar multiplication on Jacobians with fast Kummers

The costs of our key subroutines for fast Kummers are summarized in Table 5. Proof. Take the values from Table 5 in Lemma 1 for the first part and Theorem 1 for the second, using a simultaneous inversion [START_REF] Peter | Speeding the Pollard and elliptic curve methods of factorization[END_REF] to replace the 3I (corresponding to the 3 Project operations) by 6M + 1I. ⊓ ⊔ Remark 6. As we noted in Remark 1, similar techniques appear in Lubicz and Robert [START_REF] Lubicz | Arithmetic on abelian and kummer varieties[END_REF] for general abelian varieties in higher dimension embedded in projective space via theta functions. After suitable precomputations, their "compatible addition" operation views the results of Montgomery-like pseudomultiplication algorithms on a Kummer variety K of dimension g as points on the corresponding abelian variety A embedded in K × K, and therefore chooses the "correct" result of an addition on K. They explain how to map the resulting point into the 4Θ embedding of A (in P 4 g -1 ; for genus 2, this is P 15 ). A major difference with our treatment here is that Robert and Lubicz cannot treat A as a Jacobian (since general abelian varieties of dimension g > 3 are not Jacobians); hence, when specializing to genus 2, there is no connection with any curve C or Jacobian J C (and in particular, the starting and finishing points do not involve the Mumford representation). Kohel [START_REF] Kohel | Arithmetic of split kummer surfaces: Montgomery endomorphism of edwards products[END_REF] explores similar ideas for elliptic curves, leading to a very interesting interpretation of Edwards curve arithmetic.

Applications to signatures

Smart and Siksek [START_REF] Smart | A fast Diffie-Hellman protocol in genus 2[END_REF] first showed that the action of Z on G/ ±1 could be used to instantiate Diffie-Hellman on G/ ±1 rather than G. In [34, §5], they pondered whether protocols like ElGamal encryption could be instantiated on G/ ±1 , and observed that the main obstruction appeared to be that such protocols require an addition in the group-that is, a true group structure and not a mere Z-action.

Our Project-pseudomultiply-Recover technique allows a range of cryptographic protocols beyond Diffie-Hellman key exchange to take advantage of the fast and uniform (multi-)scalar multiplications offered on G/ ±1 , where G can now be any elliptic curve or any genus 2 Jacobian (defined over a large characteristic field).

In this section we illustrate this by showing how Schnorr signatures [START_REF] Schnorr | Efficient signature generation by smart cards[END_REF] can take advantage of the fast arithmetic available when G/ ±1 is the fast Kummer surface K fast C . More specifically, we describe an instantiation of a Schnorr signature scheme where G is the Jacobian of the Gaudry-Schost curve [START_REF] Gaudry | Genus 2 point counting over prime fields[END_REF] used recently to set Diffie-Hellman speed records [START_REF] Daniel | Kummer strikes back: New DH speed records[END_REF] at the 128-bit security level. Following the design of EdDSA signatures [START_REF] Daniel | High-speed high-security signatures[END_REF], we hash both the message M and signer's public key Q alongside the first half of the signature. We emphasise, however, that both of these choices are merely for illustrative purposes, and that any known variant of ElGamal signatures (Schnorr or otherwise) can use G/ ±1 for any suitable choice of G.

The Gaudry-Schost Jacobian

Let q = 2 127 -1, choose an α in F q such that 363α 2 + 833 = 0, define the six constants a = 11 , b = -22 , c = -19 , c = -3 , e = 1 + α , f = 1 -α ; setting λ := (ac)/(bd), µ := (ce)/(df ), and ν := (ae)/(bf ), we let C be the genus 2 curve over F q defined by C : y 2 = x(x -1)(x -λ)(x -µ)(x -ν) .

The Jacobian J C of C has cardinality #J C (F q ) = 2 4 N , where N = 2 250 -0x334D69820C75294D2C27FC9F9A154FF47730B4B840C05BD is a 250-bit prime.

We define the 127-bit encoding of F q to be the little-endian encoding of {0, 1, . . . , 2 127 -2}, and let H : {0, 1} * → {0, 1} 512 be any suitable hash function with a 512-bit output (eg. SHA-512 or SHA3-512).

To ensure that all our scalars are positive and of fixed bitlength, when computing scalar multiplications [m]P on J C we systematically replace m with m ′ := (m mod N ) + 3N ; then m ′ always has exactly 252 bits (since both 3N and 4N have 252 bits). Following the lead of [START_REF] Daniel | High-speed high-security signatures[END_REF], the cofactor 16 is included in the key generation and verification operations to void any potential threat of small subgroup attacks [START_REF] Hoon | A key recovery attack on discrete log-based schemes using a prime order subgroupp[END_REF]. In this case, the scalar 16z is parsed by appending four zeroes to the parsing of z described above, so a multiplication by the 252-bit scalar z is followed by 4 doubling operations. Where applicable, the operation counts stated below include the cost of these 4 additional DBL or xDBL operations.

To compare our Project-pseudomultiply-Recover approach to traditional arithmetic in J C , we take the operation counts from [START_REF] Hisil | Jacobian coordinates on genus 2 curves[END_REF]Table 2], where ADD costs 41M+7S+22a, DBL costs 26M+8S+2m c +25a, mADD costs 32M+5S+22a and mDBLADD costs 57M + 8S + 42a. Here mADD and mDBLADD compute P ⊕ Q and [2]P ⊕ Q where P is projective and Q (which is typically a fixed lookup table element) is normalized. We assume a fixed, signed window approach to computing one-dimensional scalar multiplications in J C ; our experiments with 252-bit scalars found the window width w = 5 to be optimal. In this case the scalar multiplications need 8 × DBL + 7 × mADD operations to build the lookup table of 16 elements, 1I + 39M to normalize its entires, 200 × DBL+ 50 × mDBLADD operations in the main loop, and 1I + 4M to normalize the output. For the two-dimensional scalar multiplications, we assume the standard approach to such multiexponentiations (cf. [18, §3]). Our experiments showed that w = 2 is optimal here, meaning the lookup table contains 16 elements [u]P + [v]Q for 0 ≤ u, v ≤ 3. In this case two-dimensional scalar multiplications need 3 × DBL + 10 × mADD operations to build the lookup table, 1I+36M to normalize its entries, 125×DBL+125×mDBLADD operations in the main loop, and 1I + 4M to normalize the output.

We only present counts for the dominant operations, i.e., the (multi)scalar multiplications, below. For the exact costs of compression and decompression, see

Lemma 2 .

 2 With the notation above: Suppose (a, b) = (0, 0), and write o, e, m for the last three terms of C d (a, b), so C D (A, B) = (. . .

  Suppose we have a, b, c, d, e, and f in F q \ {0} such that A := a + b + c + d , B := a + b -c -d , C := a -b + c -d , D := a -b -c + d

- 1 :

 1 K fast C → K gen C is defined by the matrix (M (a, b, c, d) t ) -1 . Since evaluating these linear transformations only involves multiplications by constants, evaluating τ to map a point from K gen C to K fast C (or τ -1 to map a point from K fast C back to K gen C ) costs 16M + 12a. Since the coordinates on both K gen C and K fast C are projective, we can save 1M by scaling one of the entries in the transformation matrix to 1.

  Let m be a non-negative integer of bitlength β, and P an element of G. Algorithm 1 computes [m]P using one call to Project, one call to xDBL, β -1 calls to xDBLADD, and one call to Recover.

		and an element P
		of G Output: R = [m]P
	1 xP ← Project(P ) ; 2 (t1, t2) ← (xP , xDBL(xP )) ; 3 for i = β -2 down to 0 do 4 if mi = 0 then
	5 6	(t1, t2) ← xDBLADD(t1, t2, xP ) ; else
	7 8	(t2, t1) ← xDBLADD(t2, t1, xP ) ; end
	9 end
	10 R ← Recover(P, t1, t2) ; 11 return R ;
	Algorithm 1: A one-dimensional uniform scalar multiplication template,
	based on the Montgomery ladder
	Lemma 1.

Table 3 .

 3 Costs of operations for projective Montgomery models By 2 = x(x 2 + Ax + 1), where mA denotes multiplications by A and (A+2)/4, and mB denotes multiplications by B. A point is normalized if its projective Z-coordinate has been scaled to 1.

	Algorithm M S mA mB a I Conditions
	ADD 1 1	0	0	5 1 P and Q normalized
	Project 0 0	0	0	0 0 -
	xDBL 2 5	1	0	0 0 -
	xADD 4 2	0	0	6 0 -
	xADD 3 2 xDBLADD 6 4	0 1	0 0	6 0 P ⊖ Q normalized 8 0 -
	xDBLADD 5 4 Recover 13 1	1 1	0 1	8 0 P ⊖ Q normalized 8 1 P normalized

Table 4 .

 4 Costs of operations for projective short Weierstrass models y 2 = x 3 + ax + b, where m b denotes multiplications by b, 2b and 4b. A point is normalized if its projective Z-coordinate has been scaled to 1.

	Algorithm M S ma m b a I Conditions
	ADD 1 1 0	0 4 1 P and Q normalized
	Project 0 0 0	0 0 0 -
	xDBL 2 5 1	2 8 0 -
	xADD 7 2 1	1 6 0 -
	xADD 6 2 1 xDBLADD 9 7 2	1 4 0 P ⊖ Q normalized 3 12 0 -
	xDBLADD 8 7 2 Recover 11 2 1	3 12 0 P ⊖ Q normalized 1 7 1 P normalized

1. Let P be a point in E(F q ) \ E

[START_REF] Daniel | Differential addition chains[END_REF]

, and let m be a positive β-bit integer. Then Algorithm 1 computes [m]P in

(8β + 5)M + 7βS + 2βm a + 3βm b + (12β + 3)a + 1I .

2. Let P and Q be points in E(F q ) \ E

[START_REF] Daniel | Differential addition chains[END_REF]

, and let m and n be positive β-bit integers. Then Algorithm 3 computes [n]P ⊕ [n]Q in (14β + 12)M + (9β + 3)S + (3β + 1)m a + (4β + 1)m b + (16β + 11)a + 2I .

Table 5 .

 5 Costs of operations for fast Kummers, where mc is used to denote multiplications by the theta constants (x0, x ′ 0 , etc) or by the curve constants (fi). Normalization refers to the respective projective coordinates.Theorem 4. Let J C be the Jacobian of a genus 2 curve admitting a fast Kummer surface, as in §7.1.1. Let P be a point in J C (F q )\ J C[START_REF] Daniel | Differential addition chains[END_REF], and let m be a positive β-bit integer. Let P and Q be points in J C (F q ) \ J C[START_REF] Daniel | Differential addition chains[END_REF], and let m and n be positive β-bit integers. Then Algorithm 3 computes [m]P ⊕ [n]Q in

	Algorithm M S mc a I Conditions
	ADD 22 2 0 27 1 P and Q normalized
	Project 11 1 3 12 1 -
	xDBL 0	8 6 16 0 -
	xADD 14 4 3 24 0 -
	xADD 7 xDBLADD 17 9 6 32 0 -4 3 24 0 x(Q ⊖ P ) fixed
	xDBLADD 10 9 6 32 0 x(Q ⊖ P ) fixed Recover 133 12 7 97 1 P normalized
		Then
	Algorithm 1 computes [m]P in
	(10β + 134)M + (9β + 12)S + (6β + 10)m c + (32β + 93)a + 2I .

Table 6 .

 6 Costs of key generation, signing and verification in the above Schnorr signature scheme.

	Algorithm	Method	M	S	mc	a	I
	Key Generation	7554 2 2654 2312 1546 8221 2 fixed window mult. on JC 8629 2131 424 Algorithm 1
	Signing	7454 2 2654 2280 1522 8157 2 fixed window mult. on JC 8525 2099 416 Algorithm 1
	Verification	multiscalar mult. on JC Algorithm 3	10813 2074 256 4478 3325 2308 14272 2 8670 2

case (0, 1, 0) : (M, (E, O)) ← (xQ, xDBLADD(xQ, xP , x⊖)) ; case (0, 1, 1) : (O, (E, M)) ← (x⊕, xDBLADD(xQ, x⊕, xP )) ; case (1, 0, 0) : (O, (E, M)) ← (x⊕, xDBLADD(xP , x⊕, xQ)) ; case (1, 0, 1) : (M, (E, O)) ← (xP , xDBLADD(xP , xQ, x⊖)) ; case (1, 1, 0) : (O, (E, M)) ← (x⊕, xDBLADD(x⊕, xP , xQ)) ; case (1, 1, 1) : (O, (E, M)) ← (x⊕, xDBLADD(x⊕, xQ, xP )) ;endsw for i ← β -2 down to 0 do switch ci do case (0, 0, 0, 0) : (O, (E, M)) ← (xADD(O, E, x⊕), xDBLADD(E, M, xQ)) ; case (0, 0, 0, 1) : (O, (E, M)) ← (xADD(O, E, x⊕), xDBLADD(E, M, xP )) ; case (0, 0, 1, 0) : (O, (E, M)) ← (xADD(O, E, x⊕), xDBLADD(M, E, xQ)) ; case (0, 0, 1, 1) : (O, (E, M)) ← (xADD(O, E, x⊕), xDBLADD(M, O, xP )) ; case (0, 1, 0, 0) : (O, (E, M)) ← (xADD(O, E, x⊕), xDBLADD(M, O, xQ)) ; case (0, 1, 0, 1) : (O, (E, M)) ← (xADD(O, E, x⊕), xDBLADD(M, E, xP )) ; case (0, 1, 1, 0) : (O, (E, M)) ← (xADD(O, E, x⊕), xDBLADD(O, M, xQ)) ; case (0, 1, 1,1) : (O, (E, M)) ← (xADD(O, E, x⊕), xDBLADD(O, M, xP )) ; case (1, 0, 0, 0) : (O, (E, M)) ← (xADD(O, E, x⊖), xDBLADD(E, M, xQ)) ; case (1, 0, 0, 1) : (O, (E, M)) ← (xADD(O, E, x⊖), xDBLADD(E, M, xP )) ; case (1, 0, 1, 0) : (O, (E, M)) ← (xADD(O, E, x⊖), xDBLADD(M, E, xQ)) ; case (1, 0, 1, 1) : (O, (E, M)) ← (xADD(O, E, x⊖), xDBLADD(M, O, xP )) ; case (1, 1, 0, 0) : (O, (E, M)) ← (xADD(O, E, x⊖), xDBLADD(M, O, xQ)) ; case (1, 1, 0, 1) : (O, (E, M)) ← (xADD(O, E, x⊖), xDBLADD(M, E, xP )) ; case (1, 1, 1, 0) : (O, (E, M)) ← (xADD(O, E, x⊖), xDBLADD(O, M, xQ)) ; case (1, 1, 1, 1) : (O, (E, M)) ← (xADD(O, E, x⊖), xDBLADD(O, M, xP )) ; endsw end switch (m0, n0) do case (0, 0) : R ← Recover(S, E, O) ; case (0, 1) : R ← Recover(P, M, O) ; case (1, 0) : R ← Recover(P, M, E) ; case (1, 1) : R ← Recover(S, O, E) ; endsw return R ;

We suppose that the inputs x(Q), x(Q ⊕ P ), and x(Q ⊖ P ) are normalized here, to simplify the exposition. In practice, we use projectivized forms of these formulae (which corresponds to replacing ξ2, ξ3, and ξ4 with ξ2/ξ1, ξ3/ξ1, and ξ4/ξ1, etc.), in order to handle the non-normalized inputs that we encounter at the end of our addition chains. The input point P in JC can be presumed to be normalized, since it is the input to the scalar multiplication routine.

E.g., compute √ u and 1/v via w ← uv 2 , w ← w -1/2 , then ( √ u, 1/v) ← (uvw, uw).

This faster compression/decompression answers a question posed by Bernstein[START_REF] Daniel | Elliptic vs. hyperelliptic, part I. Talk at ECC 2006[END_REF].

Appendix A; we omit details of point validation in J C and K fast C , both of which are essentially for free. For simplicity, we compare key generation, signing, and verification operations in §8.2 assuming that no precomputation is performed offline, noting that, if space permits it, all of these operations can take advantage of precomputation in practice.

Schnorr signatures on the Gaudry-Schost Jacobian

Let the public generator, P , be any point of order N in J C (F q ). Including P , generic elements We summarize the costs of key generation, signing and verification in Table 6. In addition to the speed benefits of Algorithm 1 and Algorithm 3 in this scenario, we reiterate that their working on K fast C avoids the side-channel issues discussed in §6.1.

A Point compression in genus 2

Here we show how to compress points on K gen C , K fast C and J C . We use E to denote a fixed exponentiation in the underlying finite field, e.g., to compute a square root. It is commonly the case that E ≈ I.

Compressing points on the general Kummer K gen C . Generic points P = (ξ 1 : ξ 2 : ξ 3 : ξ 4 ) in K gen C (F q ) ⊂ P 3 (F q ) have ξ 1 = 0, so we can compute a normalized representative (1 : k 2 : k 3 : k 4 ), where each k i = ξ i /ξ 1 , at a cost of 3M + 1I. We can further compress P to the data (k 2 , k 3 , bit), where bit is a single bit, as follows: the defining equation of K gen C (Eq. ( 4)) is quadratic in ξ 4 , so

which means that k 4 can be recovered (during decompression) from k 2 , k 3 , and the "sign" of the square root. To compress, after computing (k 2 , k 3 , k 4 ), we compute

) and set bit to be the sign bit 5 of 2K 2 k 4 + K 1 . Computing the sign bit costs 7M + 2S + 11a.

To decompress (k 2 , k 3 , bit) to a point on K gen C , we first compute

We can then use a 5 When q is a large prime, the sign bit of an element is typically chosen as its parity when represented as an integer in [0, q -1).

simultaneous inversion-and-square-root 6 to compute both K 2 1 -4K 0 K 2 and 1/(2K 2 ) using one field exponentiation, before recovering k 4 via ( 16). The decompressed point is then (1 : k 2 : k 3 : k 4 ); the decompression costs 25M + 4S + 24a + 1E.

Compressing points on the fast Kummer K fast C . As with the general Kummer K gen C , generic points P = (X : Y : Z : T ) in K fast C (F q ) ⊂ P 3 (F q ) have a nonzero first coordinate X, so we might begin compressing P by normalizing the X-coordinate to 1. But the defining equation of K fast C (Eq. ( 15)) is quartic in all of its variables, which suggests that compressing another coordinate would require solving an (unwieldy) during decompression. A faster approach to achieving further compression 7 is to map (X : Y : Z : T ) to the general Kummer, and then to normalize before performing compression/decompression there. The cost of compressing P = (X : Y : Z : T ) into two F q elements and a single bit is therefore 25M + 2S + 23a + 1I, and the cost of decompression 40M + 4S + 36a + 1E.

Compressing points on the Jacobian. To compress the four Mumford coordinates, we follow Stahlke's technique [START_REF] Stahlke | Point compression on Jacobians of hyperelliptic curves over Fq[END_REF]. For general genus 2 curves, compression of x 2 + a 1 x + a 0 , b 1 x + b 0 into (a 1 , a 0 , bit 1 , bit 0 ) costs 3M + 1S + 4a, while decompression costs 36M + 5S + 45a + 2E. Compression costs the same for genus 2 curves in Rosenhain form, but decompression is significantly faster and requires 18M + 4S + 27a + 2E. In the latter case, we compress by setting bit 1 and bit 0 as the least significant bits of 4(a

) and b 1 respectively. For decompression, we recover b 1 and b 0 from (a 1 , a 0 , bit 1 , bit 0 ) by first computing

, and B = 2(f 4 a 0 (2a 0 -a 2 1 ) + a 0 (f 3 a 1 -2f 2 ) + a 1 (f 1 + a 0 (a 2 1 -3a 0 ))), and using bit 1 to choose z 0 as the correct root z 0 = √ B 2 -4AC. We then set z 1 = (z 0 -B)/(2A); here the square root and the inversion of 2A can again be combined into one exponentiation. We can then recover b 1 as b 1 = f 4 (a 2 1 -a 0 ) -a 1 (f 3 + a 2 1 -2a 0 ) + f 2 + z 1 , using bit 0 to choose the root. Finally, we recover b 0 as b 0 = (a 0 (f 4 a 1 -f 3 -q 2 + a 0 ) + qz 1 + f 1 )/(2b 1 ), noting again that the square root and the inversion can be computed simultaneously.