A Prime-Order Group with Complete Formulas from Even-Order Elliptic Curves

. This paper describes a generic methodology for obtaining unified, and then complete formulas for a prime-order group abstraction homomorphic to a subgroup of an elliptic curve with even order. The method is applicable to any curve with even order, in finite fields of both even and odd characteristic; it is most efficient on curves with order equal to 2 modulo 4, dubbed “double-odd curves”. In large characteristic fields, we obtain doubling formulas with cost as low as 1M+5S, and the resulting group allows building schemes such as signatures that outperform existing fast solutions, e.g. Ed25519. In binary fields, the obtained formulas are not only complete but also faster than previously known incomplete formulas; we can sign and verify in as low as 18k and 27k cycles on x86 CPUs, respectively.


Introduction
Many cryptographic protocols, starting with the classic Diffie-Hellman key exchange [DH76] and Schnorr signatures [Sch90], can be expressed over a prime-order group abstraction.Elliptic curves defined over finite fields have a group structure which provides the required properties: • Elements can be serialized into a canonical representation, and deserialized efficiently.
• The group law (hereafter denoted additively) can be efficiently computed between any two elements.Similarly, the opposite of an element is easily obtained.
• Curves with prime order (thus being cyclic groups) can be found.
• Solving discrete logarithm in that group can be made computationally infeasible with curves of moderate size.
An elliptic curve over a finite field can always be described as a set of points (x, y) that are solution to a degree-3 polynomial equation known as the Weierstraß equation.The group law can be expressed as a pair of rational functions over the coordinates of the operands, and yielding the two coordinates of the result.Unfortunately, these expressions work only for most curve points; there are a few exceptional inputs for which they lead to divisions by zero: • The neutral element of the group is the "point-at-infinity" which does not have defined affine coordinates x and y.Any application of the group law that involves this point as input (adding the point-at-infinity to another point) or as output (adding a point to its opposite) is thus an exceptional case.
• Adding a point to itself is also an exceptional situation: the normal law expression involves computing the slope of the line joining the two points, which must be replaced with the tangent to the curve when the two points are equal to each other.
Such exceptional points are problematic for cryptographic implementations, because their proper handling requires testing for the exceptional cases and using alternate formulas, which would naturally lead to data-dependent variations in execution time and other similar side-channel leaks.Alternatively, an implementation would run the normal and exceptional formulas systematically, and perform a side-channel-free selection of the appropriate result, but that would imply a substantial increase in the computing cost.A better solution would be to find an alternate representation of points, and formulas that do not have exceptional inputs; following the terminology in [BL07], we will say that formulas with no exceptional case are complete, while formulas for which only the point-at-infinity is exceptional are unified (i.e.unified formulas at least properly handle the addition of a point with itself) 1 .Projective coordinates (X:Y :Z), such that x = X/Z and y = Y /Z, allow a representation of the point-at-infinity with a triplet of coordinates such that Z = 0.It has been shown by Arène, Kohel and Ritzenthaler [AKR12] that this is not, in all generality, sufficient to achieve formula completeness: any set of formulas will have exceptional inputs when considered over the algebraic closure of the field over which the curve is defined.However, for practical usages, we only use finite fields, and complete elliptic curve point formulas are then feasible, even if the same formulas would fail if the curve were lifted into an extension of the base field.Renes, Costello and Batina [RCB16] have previously published formulas for short Weierstraß curves over fields of caracteristic 5 or more, using projective coordinates; the formulas are complete as long as the difference between the two input operands is not a point of order two.Classic standard elliptic curves such as NIST's P-256 have odd order, and thus do not have any point of order two.
The Renes-Costello-Batina formulas have a non-negligible runtime cost (12M in general).Faster operations can be achieved with twisted Edwards curves [BBJ + 08, HWCD08], an alternate normal form of curves with formulas that lower the cost of point addition to 8M, and, more importantly, can be complete (depending on the exact curve definition).The well-known Curve25519, as used in Ed25519 signatures [BDL + 11], leverages this curve representation.The order of a twisted Edwards curve, however, is necessarily a multiple of four; such a curve cannot be a prime-order group.At best, we can have a curve of order f r for a (large) prime r and a (small) cofactor f , and try to work in the cyclic subgroup of points of r-torsion.Input points might however lie outside of that subgroup, leading to potential trouble in some protocols, because that allows the existence of distinct curve points which are still "equivalent" in the sense that they differ only by a small order point [CJ19].Even in the simple case of signatures, the non-trivial cofactor f can lead to different implementations disagreeing on whether a given signature (maliciously crafted by the signer) is valid or not, which can break consensus protocols [CGN20].
A prime-order group abstraction, avoiding cofactor issues, can be built on top of a (twisted) Edwards curve with cofactor f = 4 with the Decaf construction [Ham15].An extension of that method called Ristretto is applicable to cofactor f = 8, which covers the case of Curve25519 [ALdV, dVGH + 23].Decaf and Ristretto offer all the features expected from a prime-order group abstraction, albeit with relatively complex encoding and decoding procedures.
In this paper, we describe a general methodology to obtain a prime-order group abstraction out of elliptic curves with even order, with efficient formulas.The core idea consists in representing an r-torsion point P (for a prime r) by the point P + N , with N being a point of order two on the curve.The idea maps especially well to curves with order equal to 2 modulo 4, dubbed double-odd curves; however, it can be applied to other elliptic curves with an even order.With suitable choices of coordinate systems, coupled with a validation process for incoming points to verify that they have the expected P + N format, we obtain complete formulas that provide performance on par with or even faster than existing fast curves such as Curve25519, with canonical encoding into and decoding from a short, compressed format.We cover all finite field characteristics (2, 3, 5 and more).We thus obtain appropriate prime-order groups on which complex protocols can be built.As an illustration, we implemented both key exchange and signatures; with binary curve GLS254, we obtain short signatures (48 bytes at the 128-bit security level), with signature generation and verification in as little as 18k and 27k cycles on an Intel x86 CPU, respectively, which is 3 to 4 times faster than the best reported Ed25519 timings on the same platform.

Notations
In all of the following, we work in a finite field F q , with q = p m for a prime p (the field characteristic) and an integer m ≥ 1.The case p = 2, dubbed "binary field", is handled in section 5; the rest of this paper uses p ≥ 3. When m = 1, and thus q = p is prime, this is called a "prime field", which is canonically defined as integers taken modulo q.
We denote by z ∈ QR the fact that z is a square in F q , i.e. there exists some element s ∈ F q such that s 2 = z.s is a square root of z, denoted √ z.A non-zero square has two square roots, and the notation √ z does not distinguish between the two; square root computations appear only in contexts where an extra normalization step is applied to choose a specific root.Non-square z are denoted z / ∈ QR.Note that zero is a square.In binary fields, every element is a square and has a single square root, and the QR notation does not apply.
Elliptic curves are defined through equations that link point coordinates with some constant values, usually denoted a and b.These constants are assumed to be chosen such that multiplication by a or b is inexpensive.Practical implementations are expected to use some kind of fractional representation (e.g.projective coordinates) in which all divisions are mutualized into a single inversion that is typically performed only when serializing a point into bits; thus, most of the cost of individual operations (point additions, doublings...) can be approximated by the number of multiplications (M) and squarings (S) that they involve, while neglecting extra "inexpensive" operations such as additions, subtractions, and multiplications by constants derived from a and b.Squarings are typically faster than multiplications, though the ratio between the two depends on the used hardware platform, and the notion of the individual cost of an elementary operation is not well-defined in modern pipelined processors, especially superscalar CPUs that support out-of-order execution.As a rule of thumb, a squaring cost can be estimated to be about 80% the cost of a multiplication.
In the case of binary curves (section 5), the method is applicable to existing standard curves, and some standard NIST curves (the B-* non-Koblitz curves) use pseudorandom values for the b constant, not chosen to make multiplications fast.To cover that extra cost, we denote by m b the cost of multiplying by b (or a constant derived from b) in estimates for binary curves.

Geometrical Description
The core idea is to use a point of order two as pivot.Let F q be a finite field, and E an elliptic curve defined over that field with a degree-3 equation (i.e. a Weierstraß curve).Let n be the order of E; by Hasse's theorem, we know that n is close to q (specifically, We suppose that n is even, and r is a large prime divisor of n.The points of r-torsion in E are: where O is the point-at-infinity (the neutral element of the group).E[r] is a subgroup with order r or r 2 ; the latter case may happen only if r 2 divides n, and, in general, we prefer to avoid it, except in the context of pairing-friendly curves.We will hereafter suppose that r 2 does not divide n.The subgroup E[r] is the prime-order group structure over which we wish to build our abstraction.The group law over the curve is defined as follows: • The point-at-infinity O is the neutral element.
• For any given point P = (x, y) on the curve, there is a single other point P ′ with the same x coordinate; we define that other point to be the opposite of P , denoted −P .
• For two points P and Q on the curve, the line (P Q) intersects the curve in a third point R; we define P + Q = −R.The point R is obtained from −R by mirroring it against the x axis.
This definition highlights the exceptional points, where the formulation above does not work properly: the point-at-infinity has no defined coordinates, making the line (P O) ill-defined; and if P and Q are the same point, the line (P Q) is similarly underspecified (we use the tangent to the curve in that case).Since we supposed that n is even, we can write n = f r for some even cofactor f .There must be at least one point of order two, which we call N (since E[2] over the algebraic closure of F q is homomorphic to Z 2 × Z 2 , there are either one or three points of order two).By definition, N is not part of E[r].The main idea is to represent a point P ∈ E[r] by the point P + N (which is not in E[r]).This has the following advantages: • The neutral O is now represented by the point O + N = N , which has defined coordinates.
• The sum of P and Q, represented by P + N and Q + N , respectively, will be itself represented by (P + Q) + N .The latter point can be computed as: with the important remark that if P and Q are in and, in particular, Q + N cannot be the same point as P .
In this way, we can get at least unified formulas: when adding two points P and Q, known through their representants P + N and Q + N , we first apply the group law on P + N and N , to yield the point P , then we apply the law again, between P and Q + N ; in neither case, the two operands can be the same point, which fully avoids the exceptional case of adding a point to itself.
A remaining potential issue is when the first operand is O; if P = O and is represented by P + N = N , then the addition of N to N is an exceptional case by itself, both because it adds a point to itself, and because the output is the point-at-infinity.However, as will be detailed in the next section, that case can be eliminated in the formulas themselves.

Analytical Description
With appropriate changes of variables that do not modify the curve's algebraic shape (i.e.reversible affine transforms in the plane), the elliptic curve E can always be described by a Weierstraß equation; the elements of E are the point-at-infinity, and the set of points (x, y) ∈ F q × F q such that: for some constants a 1 , a 2 , a 3 , a 4 and a 6 .For a point P = (x, y) ∈ E, its opposite is −P = (x, −y − a 1 x − a 3 ).The geometrical definition translates into the following, when adding points P 1 = (x 1 , y 1 ) and P 2 = (x 2 , y 2 ) together: • The slope of the line (P 1 P 2 ) is computed as: If P 1 and P 2 have the same x coordinate, then the formula above fails.In that case, then either y 2 = −y 1 − a 1 x − a 3 , which implies that P 2 = −P 1 , and their sum is then the point-at-infinity O; or P 1 = P 2 ̸ = −P 1 , and we replace the line (P 1 P 2 ) with the tangent to the curve at that point, and its slope is: • The coordinates of P 3 = (x 3 , y 3 ) = P 1 + P 2 are then obtained as: These equations are the "classic affine addition formulas", which are detailed in many textbooks (e.g.[Was08]).
From that point, we consider that F q has characteristic 3 or more; the case of fields with characteristic 2 will be covered in section 5. We can therefore divide values by 2.
Since we assumed that the curve E has even order, and chose a point N = (x N , y N ) of order two, we can do a further change of variable (x, y) → (x + x N , y + y N ), which corresponds to moving the coordinate reference point to N itself.This does not change the general shape of the equation, even though the constants are modified; we can now assume, without loss of generality, that N = (0, 0).Substituting these coordinates into the Weierstraß equation yields a 6 = 0. Similarly, N being of order two, it must be that −N = N , and this implies that a 3 = 0. Since we consider here a field of characteristic distinct from 2, we can apply a further change of variable y → y + (a 1 x)/2, which modifies the equation into the following: 1 /4 and b = a 4 .The following notes apply: • This shape is not the classical "short Weierstraß" equation.It can be obtained under the hypothesis that there is a point of order two, i.e. that the curve order is even.
• In descriptions of elliptic curves, the case of a field of characteristic 3 is normally distinguished, due to the inability to divide by 3 in such a field.As a consequence of the assumption of an even curve order, we do not need to make that distinction here; in this paper, fields of characteristic 3 can use the same treatment and formulas as fields of characteristic 5 and more.
• The group law definition requires the existence of a well-defined tangent on each point, i.e. that the two partial derivatives 2y and 3x 2 + 2ax + b do not vanish simultaneously on any curve point (if such a case happens, the curve is said to be singular).With the curve shape above, the curve is non-singular if and only if b ̸ = 0 and a 2 − 4b ̸ = 0.
• The j-invariant of the curve is: • N is the only curve point such that x = 0.
For a curve point P = (x, y) ̸ = N, O, the point P + N is obtained as: P + N = (b/x, −by/x 2 ).Applying this formula into the group law definition, a straightforward calculation yields generic addition formulas.Given points P 1 , P 2 ∈ E[r], represented by points P 1 + N = (x 1 , y 1 ) and P 2 + N = (x 2 , y 2 ), we want to obtain the representant of their sum P 3 + N = (P 1 + P 2 ) + N = (x 3 , y 3 ).We first suppose that P 1 , P 2 and P 3 are all distinct from O, which implies that x 1 and x 2 are non-zero.As explained above, under these conditions, no exceptional case occurs in the addition formulas, and we get the following expressions: During the calculation, the expressions were simplified by removing a factor x 1 x 2 common to the numerator and denominator; this removal happens to also make the formulas compatible with the cases involving O. Indeed, it can be verified that if x 1 = 0 or x 2 = 0, then the correct result is obtained; similarly, if x 1 = x 2 and y 1 = −y 2 , i.e. adding a point to its opposite, the point N is correctly obtained as the representant of the sum O.These formulas are thus complete -provided that the two inputs are indeed the correct representants of two elements of E[r].The problem of verifying that an input point is a correct representant of an element of E[r] is covered in the section 3.3.Some sub-cases of the curve equation y 2 = x 3 + ax 2 + bx have already been studied in previous works (e.g.Doche-Icart-Kohel double-oriented curves use b = 16a [DIK06]; Castryck and Decru suggest b = −1 in the context of isogeny-based cryptography [CD20]).To our knowledge, none envisioned using point N as pivot to avoid exceptional cases in formulas.

Point Validation and Double-Odd Curves
In order to provide the expected prime-order group abstraction, it is necessary to make sure that externally received inputs (e.g.public keys) are correct representants of elements of E[r].We suppose here that we have received a point Q = (x, y) and that we have verified in some way that x and y fulfill the curve equation y 2 = x 3 + ax 2 + bx (the point is on the curve).
Torsion Check.A conceptually simple method to validate the point Q is to multiply it by the scalar r; for points of r-torsion, this must yield the neutral O (represented by N ).A potential difficulty is that the formulas shown in section 3.2 are guaranteed complete only as long as the operands are correct representants of r-torsion points, which is not yet known to be the case for the point Q.Analysis of the formulas shows that they fail in a single case: when inputs (x 1 , y 1 ) and (x 2 , y 2 ) are such that (x 1 , y 1 ) = (x 2 , y 2 ) + N .In such a case, all numerators and denominators vanish in the expressions of x 3 and y 3 .A consequence is that in actual implementations of the formulas, using any fractional representation (e.g.projective coordinates), an exceptional case will yield all-zero coordinates, and that situation will be "sticky": subsequent point additions using the formulas will remain all-zero values.The exact details vary depending on the concrete choice of coordinate system, but in general one can expect the following validation process to work: 1. Verify that the provided x and y coordinates fulfill the curve equation.
2. Assuming that (x, y) = P + N for some curve point P , compute rP + N with any double-and-add variant leveraging the formulas from section 3.2 with a fractional/projective concrete representation of coordinates.
3. If an invalid (all-zero) representation is obtained, then an exceptional case was hit and the point (x, y) is not valid.Otherwise, if no such case was obtained but rP + N is not the point N itself, then (x, y) is not valid.If rP + N yields a proper representation of N , then the input (x, y) is a valid representant of a point This process can work for any elliptic curve with an even order.It is unfortunately relatively expensive (the computation of rP can use the fact that r is not secret and is known in advance, thus allowing the precomputation of an efficient addition chain, but the cost will still be commensurate with that of a generic multiplication of a point by a scalar).However, for some curves, a substantially faster process can be used, as detailed below.
Point Halving and Double-Odd Curves Suppose that the cofactor is a power of two, i.e. that the whole curve order is equal to 2 t r for an odd prime r and an integer t ≥ 1.Any curve point can be uniquely written as the sum of a point of r-torsion and a point of 2 t -torsion.Let 2 h be the maximum order of all 2 t -torsion points (we have h ≤ t, but not necessarily h = t).Suppose that a given point Q can be halved h times, i.e. that Q = 2 h R for some other curve point R (not necessarily unique).We can write R as R = P + T , with P ∈ E[r] and T ∈ E[2 t ]; thus: Conversely, if Q is a point of r-torsion, then it can be halved indefinitely, by multiplying it with the inverse of 2 modulo r.We therefore obtain a test for membership to the subgroup of r-torsion points: these are exactly the points that can be halved at least h times.Point halving, in general, is a somewhat expensive operation, since it involves two square roots, and each square root in the field is (roughly) equivalent to about 1/10th of the cost of a point multiplication by a scalar.Moreover, each halving yields at least two solutions, and one must select the "right one", which is complicated if E[2 t ] is not cyclic.However, there are curves for which the validation process is efficient, namely if t = 1.These are curves of order 2r, i.e. twice an odd integer; we call such curves double-odd: Definition 1.A double-odd elliptic curve is a non-singular elliptic curve defined over a given finite field, such that its order is equal to 2 modulo 4.
An elliptic curve of equation y 2 = x 3 + ax 2 + bx in a field of characteristic 3 or more is double-odd if and only if both b and a 2 − 4b are non-squares in the field.Indeed: • If a 2 − 4b ∈ QR then the polynomial X 3 + aX 2 + bX has three roots, leading to three points of order two in the curve.These three points, together with O, make a subgroup homomorphic to Z 2 × Z 2 , whose order is 4 and must then divide the curve order.
• If a 2 − 4b / ∈ QR but b = c 2 for some constant c, then a 2 − 4b = (a + 2c)(a − 2c), which implies that one of a + 2c and a − 2c must be a square.We can assume without loss of generality that a + 2c ∈ QR (otherwise, we just replace c with the other root −c).Then it can be easily verified that T = (c, c √ a + 2c) is a valid curve point whose double is N ; this is therefore a point of order 4, and the curve order is a multiple of 4.

• Conversely, if b /
∈ QR and a 2 − 4b / ∈ QR, then N is the only point of order two (the polynomial X 3 + aX 2 + bX has a single root, which is zero); E not being double-odd would then imply the existence of a point M = (x m , y m ) such that 2M = N , which would imply that M + N = −M .As was pointed out previously, the x coordinate of −M is the same as the x coordinate of M (i.e.x m ), and the x coordinate of M + N is b/x m , leading to x m = b/x m , which is not possible when b / ∈ QR.
This is straightforwardly derived from the classic doubling formulas; one may also apply the formulas from section 3.2 and obtain 2P + N = (b/x ′ , −by ′ /x ′2 ).The x coordinate of the double of any point is therefore a quadratic residue.Conversely, on a double-odd curve (thus with b / ∈ QR), if P = (x, y) then the x coordinate of P + N is b/x, which is a square if and only if x is not a square.This leads to the following property: in a double-odd curve with equation y 2 = x 3 + ax 2 + bx, the curve points segregate into: • the point-at-infinity O; • the unique point N of order two, whose x coordinate is zero; • the non-infinity points of r-torsion, whose x coordinates are non-zero squares in the field; • the non-infinity points of order exactly 2r, whose x coordinates are non-squares in the field.
This finally yields an efficient point validation test for double-odd curves: a point (x, y) on the curve is the correct representant P + N of a point P ∈ E[r] if and only if either x = 0 (in which case P = O), or x / ∈ QR.This validation can be implemented with a single Legendre symbol computation, which can be done efficiently, e.g. with a binary GCD variant [Por20a] or a similar plus/minus algorithm [AHST23].
The fact that x ∈ QR if and only if (x, y) is the double of another rational point can also be deduced by noticing that x is equal to the unreduced 2-Tate pairing ⟨N, (x, y)⟩ 2 .The reduced 2-Tate pairing (after raising to exponent (q − 1)/2) is thus the Legendre symbol of x.The bilinearity of the pairing then implies the result.

Compression and Canonical Decoding
Using the pivot point N allows for a straightforward encoding of a subgroup element into a compressed format consisting of a single field element.Consider a point P ∈ E[r], represented by P + N = (x, y).The line that joins N and P + N intersects the curve in a third point, which is Therefore, for a given slope w ∈ F q , the line with slope w and passing through N can intersect the curve in at most one correct representant P + N of a point P ∈ E[r].We can thus encode a point P + N = (x, y) by mapping it to the line slope w = y/x (for the point N itself, we can use the value w N = 0, since the other points with y = 0, if any, have order two and thus are not correct representants P + N for any P ∈ E[r]).
Decoding w into the matching point representant P + N involves the following: 1.If w = 0 then the matching point is N .We can now assume that w ̸ = 0, and the solution (x, y), if any, must have x ̸ = 0 and y ̸ = 0.
2. Since w = y/x, we can rewrite the curve equation into x 2 +(a−w 2 )x+b = 0, a degree-2 equation in x that can be solved by computing its discriminant ∆ = (a − w 2 ) 2 − 4b and extracting its square root.If the discriminant is not a square then there is no solution and the input w is invalid.
3. Once x is obtained, the corresponding y is computed with y = wx.The quadratic equation above, in general, has two solutions, only one of which (at most) being a correct representant of an r-torsion point; thus, point validation should be performed to identify the correct solution.
In the case of a double-odd curve, the process in steps 2 and 3 simplifies into the following: the two solutions to the quadratic equation, if they exist, are such that their product is b, which is not a square for a double-odd curve.Thus, exactly one of the two solutions for x will be a non-square; this solution should be selected, and no further point validation is necessary.The cost of point decoding from the compressed format for a double-odd curve is mostly that of a square root and an additional Legendre symbol to select the correct solution.This process ensures a canonical encoding: a group element (represented by a point P + N ) can be encoded into a unique field element w, and the decoding process unambiguously reconstructs the corresponding source element, rejecting invalid encodings.Once decoded, the element is guaranteed to be a correct representant of an r-torsion point, and the formulas from section 3.2 can be applied without any exceptional case, thanks to their completeness.

Point Doubling and Isogenies
For this section, we denote the curve parameters explicitly: E(a, b) is the curve of equation y 2 = x 3 + ax 2 + bx.We use (x, w) coordinates, with w = y/x; the w value is the slope of the line from N to (x, y), and we already encountered it in section 3.4.Since {O, N } is a subgroup of E(a, b), we can use Vélu's formulas [V 71] to obtain an isogeny from E(a, b) into a dual curve, whose kernel is exactly {O, N }.Doing it twice, we obtain two isogenies: (We also applied an extra scaling by 1/4 in the second isogeny; Vélu's formulas would otherwise lead us to ) is double-odd, then w ̸ = 0 for all points other than N and O, and E(−2a, a 2 − 4b) is also double-odd.Moreover, for any point P ∈ E(a, b), then ψ ′ 1/2 (ψ(P )) = 2P .These isogenies can be translated into efficient point-doubling formulas with the following characteristics: • For the purposes of point doublings, we can represent a curve point P in Jacobian (x, w) coordinates with a triplet (X:W :J) such that x = X/J 2 and w = W/J, with J ̸ = 0.The special points N and O are represented with (0:W :0) and (W 2 :W :0) for any W ̸ = 0, respectively.
• Since we want to work with representants P + N , we need to include the overhead of converting the source to Jacobian (x, w) coordinates from whatever coordinate system we are otherwise using for the curve, and converting back the result with the extra addition of N after the doubling.In a sequence of doublings, "internal" doublings in Jacobian (x, w) coordinates avoid these conversion costs, so that the overhead is only a one-time cost for any sequence of successive doublings.Conversions and addition of N can be merged with the first and last doublings in a sequence.
• In Jacobian (x, w) coordinates, a point doubling and addition of N can be computed in cost 2M + 5S generically.If q = 3 mod 4, the cost can be 1M + 6S.If q = 3 mod 4, a = −1 and b = 1/2, the cost is down to 2M + 4S.The best case is when a = 0, for which a doubling (without the addition of N ) can be done in cost 1M + 5S.Moreover, all these formulas are complete.Details on the formulas are included in appendix B.

Fractional (x,u) Coordinates
The formulas in section 3.2 use the classic affine (x, y) coordinates.We can obtain other formulas amenable to faster implementations by replacing the y coordinate with another value, namely u = x/y.We choose this value u instead of its inverse w = 1/u, used above for compression, because we can naturally extend u to the value zero for the point N , where both x and y are zero: indeed, w is the slope of the line from N to a point P , and for N itself, that line should be the tangent to the curve on N , which is vertical; thus, its inverse should intuitively go to zero when reaching N .Replacing y with x/u in the previous formulas and simplifying leads to the following formulas.For two points P 1 and P 2 in E[r], represented by P 1 + N = (x 1 , u 1 ) and P 2 + N = (x 2 , u 2 ), respectively, their sum is represented by P 3 + N = P 1 + P 2 + N = (x 3 , u 3 ) with: These formulas are complete, as long as the two inputs are correct, i.e. have been validated to be proper representants of r-torsion points.In fact, in all generality, the formulas above have been derived under the assumption that P 1 + N ̸ = ±P 2 , and P 1 and P 2 are distinct from N and O.It can be verified that they still work when P 1 + N and/or P 2 + N is equal to N , and the point validation process from section 3.3 ensures that the other cases cannot happen.
In a practical implementation, we can use a fractional representation: a point (x, u) is represented by a quadruplet (X:Z:U :T ) such that Z ̸ = 0, T ̸ = 0, x = X/Z and u = U/T .This leads to general point addition formulas with cost 10M; sequences of n successive point doublings can be computed with cost as low as 3M + n(1M + 5S) for some specific curves.These processes are detailed in appendix C.

Double-Odd Jacobi Quartic
In section 3, we described a generic process by which a prime-order group abstraction could be extracted from an elliptic curve with an even order.In case the curve is double-odd, the point decompression and validation is reasonably efficient, though somewhat slower than the equivalent process in, for instance, Ristretto.The derived formulas in fractional (x, u) coordinates are also slightly less efficient than for twisted Edwards curves, except for long sequences of successive doublings.We can, however, do better, with another change of coordinates.

Map to a Jacobi Quartic
As previously, we consider a finite field F q of characteristic 3 or more, and an elliptic curve E of even order, expressed with the non-short Weierstraß equation y 2 = x 3 + ax 2 + bx for two constants a and b.We now focus on a double-odd curve, i.e. its order is 2r for some odd integer r; as explained in section 3.3, this means that both b and a 2 − 4b are non-squares.Point-at-infinity is denoted O; the point N = (0, 0) is the only point of order two on the curve.For all other curve points, both x and y coordinates are well-defined and non-zero.
Given a point (x, y) on the curve (other than N and O), we define the (e, u) coordinates as follows: Applying the curve equation, the coordinate e can also be written as: This last formula can also be computed for point N , yielding e = −1.We formally define the (e, u) coordinates of N and O to be (−1, 0) and (1, 0), respectively.
The mapping (x, y) → (e, u) is reversible, by noticing that: We can thus use (e, u) coordinates to represent all curve points without any loss of information.
In (e, u) coordinates, the curve equation becomes: which is a form known as the (extended) Jacobi quartic, first studied by Jacobi in the 19th century [Jac29].In the usual formulation, the equation is denoted Y 2 = DX 4 + 2AX 2 + 1, with subvariants depending on whether D is a square or not.In our case, e and u coordinates play the role of Y and X, respectively, and the D constant is a 2 − 4b, which is a non-square for all double-odd curves.The mapping goes in both directions: a Jacobi quartic Y 2 = DX 4 + 2AX 2 + 1 is turned into a curve of equation y 2 = x(x 2 + ax + b) with a = −A and b = (A 2 − D)/4 by setting x = (Y + 1 + AX 2 ) and y = x/X.Thus, the Jacobi quartic Y 2 = DX 4 + 2AX 2 + 1 is a double-odd curve if and only if D / ∈ QR and A 2 − D / ∈ QR.There exists some extensive literature on Jacobi quartics and their formulas [BJ03, CC86, Jac29, WW27].The following classic point addition formulas can also be derived straightforwardly from the (x, u) formulas shown in section 3.6; for points P 1 = (e 1 , u 1 ) and P 2 = (e 2 , u 2 ), their sum P 3 = (e 3 , u 3 ) = P 1 + P 2 can be computed as: These formulas are complete for all curve points: it can be easily verified that they also work for the cases on which the (x, u) formulas would have failed (note that since a 2 − 4b / ∈ QR, the denominator 1 − (a 2 − 4b)u 2 1 u 2 2 cannot be zero).In (e, u) coordinates, if P = (e, u) then P + N = (−e, −u) and −P = (e, −u) (this explains why the coordinates of O were chosen to be (1, 0), since N has coordinates (−1, 0)).
It should be noted that the above formulas were considered in the case of double-odd curves, but they are still complete for curves whose order is a multiple of 4, as long as a 2 − 4b / ∈ QR, i.e. that N is the only point of order 2. If b ∈ QR then there are points of order 4 (see section 3.3) whose x coordinate is a square root of b; for such points, the e coordinate is equal to zero.On double-odd curves, e ̸ = 0 for all curve points.

Prime-Order Group and Point Compression
Since the formulas over the Jacobi quartic are complete on a double-odd curve, we can redefine the prime-order group abstraction with a different description that leads to a more efficient point compression and decompression process.Namely, we consider the quotient group E/{O, N }, which is homomorphic to E[r]; elements of the quotient group are pairs {P, P + N }, and any two such pairs are either identical or disjoint.Equivalently, we are still working with E[r] but for each point P ∈ E[r] we allow two possible representants, which are P and P + N (whereas in section 3 we enforced the use of only P + N as a valid representant).We do not really care which representant we are using, since both work; we only need a normalizing rule for canonical encoding, which we will now define.
Let sign be an arbitrary sign function over F q , i.e. a function such that: • For any z ∈ F q , sign(z) ∈ {0, 1}.
There are many such functions; a convenient one for prime fields (fields of integers modulo a given prime q) is to represent the value z as an integer i ∈ [0, q − 1] range, and use sign(z) = i mod 2 (the "least significant bit" of the value).A value z is said to be "negative" if sign(z) = 1, or "non-negative" if sign(z) = 0. Encoding of a point P = (e, u) into an element of F q uses the following process: 1.If sign(e) = 1, then return −u, otherwise return u.
Since P + N = (−e, −u), P and P + N encode into the same output value.An equivalent description of encoding is that for a given element {P, P + N } of the quotient group, the two points P and P + N are the two possible representants, and exactly one of them has a non-negative e coordinate; we use the u coordinate of that specific representant.When encoding the group neutral {O, N }, we get zero.Decoding of a given value u into a point P is performed by recomputing a matching e coordinate, as follows:

If ∆ /
∈ QR, then there is no solution and u is invalid.Otherwise, let e = √ ∆.

Return (e, u).
This process illustrates that decoding must work for any validly encoded point, but also that it is unambiguous: if there is a matching e coordinate, then there can be at most two solutions for e, but they are opposite to each other and thus cannot be both non-negative, and the output of the decoding process is necessarily a representant of the same quotient group element as was used as input to the encoding process.
Starting from (e, u) affine coordinates, the encoding process cost is trivial; realistically, an implementation will use a fractional representation of the coordinates, and the cost of the encoding is the cost of normalizing that representation to affine coordinates, i.e. mostly one inversion in the field.This is faster than the encoding process in Decaf or Ristretto, where a combined square root and inversion must be performed: when working in a prime field, that operation requires an exponentiation with an exponent of about the same size as the modulus, whereas a single inversion, as needed for double-odd Jacobi quartics, can be performed more efficiently with a binary GCD derivative [BY19,Por20a].For decoding, a square root computation is needed, and should offer performance similar to that of Decaf or Ristretto, albeit with a somewhat simpler process since there is no combined square root/inversion operation; the decoding process inherently outputs affine (e, u) coordinates.
It may be noted that all of the above was described in the context of double-odd curves, but can be potentially applied to curves whose order is a multiple of 4, in which case either or both of the following apply: • There are two other points of order two on the curve, with y = 0 but x ̸ = 0 (in (x, y) coordinates).For these two points, the u coordinate is not well-defined; moreover, the addition formulas are incomplete, since they fail for such points as input or as output.This case corresponds to a 2 − 4b ∈ QR.
• There are two points of order four; if one of them is T then the other is T + N .These two points have the same x coordinate, and for both e = 0.They thus have two distinct coordinates u, but both have a non-negative e, and the process described above is ambiguous.This case corresponds to b ∈ QR.
The second situation can be solved with an extra normalization rule.The problem, in that case, is that the encoding process cannot choose which of T or T + N to use, and the encoding is not canonical; the decoding process still works properly, in that it will return either T or T + N but both represent the same quotient group element.In order to obtain a canonical encoding, we only need to apply an extra normalization rule, e.g.: if e = 0 and sign(u) = −1, then replace u with −u.With this extra rule, we obtain a canonical encoding of the quotient group for all even-ordered curves, as long as they contain a single point of order two (this would apply to both Curve448 and Curve25519, for instance).Of course, if the curve is not double-odd, but its order is a multiple of 4, then the quotient group has even order and we will not obtain a prime-order group abstraction that way.

Isogenies and Point Equality
In section 3.5, we defined the isogenies ψ 1 and ψ ′ 1/2 , over the curve points in (x, w) coordinates.The formulas can be straightforwardly transformed into (e, u) coordinates: For a double-odd curve, e ̸ = 0 for all curve points, and these expressions are well-defined without any exceptional case.If a canonical encoding is defined, as in section 4.2, then two points can be compared with each other by encoding both and checking whether the outputs are identical.This can be inefficient, though, as encoding typically involves normalization from a fractional representation to affine coordinates, hence at least one field inversion.A faster method can be obtained by using the ψ 1 isogeny, whose kernel is {O, N }; for any points P 1 and P 2 on the curve, the two points represent the same quotient group element if and only if ψ 1 (P 1 ) = ψ 1 (P 2 ).Moreover, ψ 1 (P 1 ) is necessarily an r-torsion point in the destination curve (E(−2a, a 2 − 4b)), and no two distinct elements of the subgroup of r-torsion points may share the same u coordinate.
We thus only need to compute the u coordinates of ψ 1 (P 1 ) and ψ 1 (P 2 ) and compare the two values together, which can be done while staying in the fractional domain.As expressed above, the two values are −u 1 /e 1 and −u 2 /e 2 , respectively; therefore, the points P 1 and P 2 represent the same quotient group element if and only if u 1 e 2 = u 2 e 1 .

Extended (e,u) Coordinates
The formulas described above can be transformed into an implementation by using several possible representations.Since the addition formulas shown in section 4.1 were already known, several previous papers explored them; the most efficient were found by Duquesne [Duq07], and also rediscovered by Hisil, Wong, Carter and Dawson [HWCD09] under different notations that ultimately resolve to the same elementary operations.We use notations close to the ones used by Hisil et al: a point (e, u) is represented by a quadruplet of values (E:Z:U :T ) such that Z ̸ = 0, e = E/Z, u = U/Z and U 2 = T Z (which implies that u 2 = T /Z).The cost is 8M + 3S; this is similar to the 10M cost for (x, u) coordinates (section 3.6).On the other hand, we can apply to double-odd curves in Jacobi quartic form the same improvements on sequences of doublings as in (x, u) coordinates, with a lower per-sequence overhead: • A single doubling can be done in cost 2M + 5S for all curves.If q = 3 mod 4, or if a = 0, then the cost can be lowered to 1M + 6S.
• The single doubling formulas internally use Jacobian (x, w) coordinates, and additional doublings can be performed with the costs expressed in section 3.5, down to 1M + 5S per doubling on curves with a = 0.
All the relevant formulas are described in appendix D. Compared with twisted Edwards curves, the double-odd Jacobi quartics offer slower generic addition but faster doublings: on twisted Edwards curves, generic point addition has cost 8M with extended coordinates [HWCD08], but single doublings in the same coordinate system have cost 4M + 4S (additional doublings can be performed at cost 3M + 4S).The overall gain or loss of Jacobi quartics depends on how many operations of each type are used in a given protocol.The multiplication of a group element by a scalar is a very common primitive, and it is typically implemented with a window-optimized double-and-add algorithm; for such an implementation, there will be four or five point doublings for every generic addition, and the faster doubling formulas make Jacobi quartics faster than twisted Edwards curves for this task.

j255e and jq255s
As a concrete example of double-odd Jacobi quartics amenable to fast software implementations, we define two curves at the usual "128-bit" security level.In both cases we use fields of integers modulo the prime q = 2 255 − m for a small integer m, for which fast implementations are feasible, especially on low-power microcontrollers, as long as m is small enough (roughly speaking, best performance is achieved with m < 2 15 ).
The jq255e curve uses a = 0, so as to leverage the fastest doublings.Such a curve has j-invariant equal to 1728.We need q = 1 mod 4 to avoid the curve being supersingular; we prefer to set q = 5 mod 8 so that square roots can still be computed with relative ease (using Atkin's formulas [Atk92]).For such curves, trying b = ±2 exhausts all possibilities; we thus look for the lowest m = 3 mod 8 such that q = 2 255 − m is prime and y 2 = x 3 + bx defines over F q an elliptic curve of order 2r for a prime r (which will be close to 2 254 ).The first hit is for m = 18651 and b = −2.Using the methodology described here, this yields a group of prime order r (which is very slightly lower than 2 254 ).As a nice bonus, such a curve is a GLV curve [GLV01], with a very efficient endomorphism that can be used to achieve substantial speed-ups for some operations, in particular multiplication of a point by a scalar (the endomorphism is simply (e, u) → (e, u √ −1)).The jq255s curve is an alternate choice such that q = 3 mod 4, a = −1 and b = 1/2: for such curves, doublings in Jacobian (x, w) coordinates can be done in cost 2M+4S.Again, we select the lowest m such that q = 2 255 − m is prime and the equation y 2 = x 3 − x 2 + 1/2 defines a curve over F q of order 2r for a prime r; the first match is m = 3957.The point of jq255s is to demonstrate that double-odd Jacobi quartics can achieve performance at least on par with twisted Edwards curves without requiring the extra structure of a GLV curve; in practical situations, though, jq255e should be preferred.
For both curves, we define and use the prime-order quotient groups, as explained in the previous sections.Over these groups, we can define efficient key exchange (Diffie-Hellman [DH76]) and signature (Schnorr [Sch90]) protocols; for the latter, we leverage the curiously neglected remark from Schnorr that the "challenge" part only needs to be half the size of the group, thus yielding signatures of size 48 bytes at the 128-bit security level, which is shorter than the commonly encountered 64 bytes from Ed25519 [BDL + 11] or standard ECDSA on the P-256 curve [ECD23] 2 .Half-size challenges still ensure the expected security [NSW09].A full specification of the jq255e and jq255s groups, along with the key exchange, signature, and hash-to-curve protocols, is available at: https://c2sp.org/jq255 Open-source implementations and a summary of the formulas are available from the dedicated Web site at: https://doubleodd.group/As an example of the achieved performance, on an Intel Core i5-8259U x86 processor ("Coffee Lake" core, in 64-bit mode), we achieve, with an optimized (and fully constanttime) implementation in Rust, a general jq255e point multiplication by a scalar in about 69k clock cycles, and a signature verification in 82k cycles.For jq255s, which does not benefit from a GLV endomorphism, the values are 103k and 83k, respectively.For Curve25519, the corresponding values are 103k and 108k 3 .The faster signature verifications for jq255e and jq255s come mostly from the use of shorter challenge values (which also make signatures shorter); the compared Ed25519 implementation leverages the Antipa et al optimization [ABG + 06, Por20b] but is still more expensive because of its longer challenge values.

Binary Curves
So far we have considered the case of finite fields with characteristic 3 or more.We now explore how our generic methodology can be applied to elliptic curves defined over fields of characteristic two; such fields and curves are often said to be "binary".
2 In the case of Ed25519, the longer signature is explicitly meant to support batch verification; the lack of use of compressed challenges is still a common feature of most Schnorr-derived schemes, and rarely explained.
3 An x-only point multiplication routine is often used with Curve25519, under the name "X25519"; it is sufficient for some protocols, e.g.most usages of Diffie-Hellman key exchange.X25519 can be faster than generic point multiplication; 88k cycles have been reported on a CPU type equivalent to our test system [x2523].However, X25519 does not validate that incoming points are part of the expected subgroup, admits multiple equivalent inputs, and does not yield the complete output point to support further arbitrary group operations; it thus departs from the prime-order group abstraction.We consider here a full point multiplication, which is somewhat more expensive.

Binary Fields
Let q = 2 m .We can define the field F q as the ring of polynomials with binary coefficients (F 2 [z]) with all computations performed modulo a given polynomial M of degree m; any M can be used as long as it is irreducible.All finite fields with the same cardinal are isomorphic to each other, and the isomorphisms can be computed in practice; therefore, no choice of M is better or worse than any other from a security point of view.We can thus choose a modulus M that favours performance (e.g. with a very low Hamming weight).We recall here some relevant properties of binary fields: • Squaring is a field automorphism: for any field elements x and y, (xy) 2 = x 2 y 2 and (x + y) 2 = x 2 + y 2 .The same applies to raising to power 2 n for any integer n, including negative values (e.g. Every field element is a square and has a unique square root.Raising to the power 2 n is a linear operation (when considering F q as a vector space of dimension m over F 2 ) and can be implemented efficiently.
• Inversion can be performed with a variant of Fermat's little theorem described by Itoh and Tsujii [IT88].This does not make inversions fast enough to contemplate using them for each curve point addition, but it still makes them fast enough to justify some practices, e.g.normalizing to affine coordinates a "window" of values for an optimized double-and-add algorithm.
• The trace of a field element x is: The trace is linear; it is always equal to 0 or 1; the trace of 0 is 0. For any x, Tr(x 2 ) = Tr(x) 2 = Tr(x).If m is odd then Tr(1) = 1.The trace of an element x can always be efficiently computed because it is equal to the sum of a few specific coefficients of x as a binary polynomial in z (which coefficients are part of the trace expression depends on the choice of the field modulus M ).Correspondingly, there is always at least one field element of trace 1 with minimal Hamming weight (i.e. a single bit).
• Quadratic equations can be reduced to either a square root computation, or to the equation x 2 + x = d for some field element d.If Tr(d) = 1 then the equation has no solution; if Tr(d) = 0 then there are two solutions; if x is a solution then x + 1 is the other solution.When m is odd, a solution can be obtained with the halftrace: In fact, for an odd m, we always have Since the halftrace is linear, it can be computed with a simple bit-by-bit process (or even faster with look-up tables to process bits by chunks).

Structure of Binary Curves
Over a binary field, elliptic curves with an odd order are supersingular, which is usually a problem for security, because discrete logarithm over supersingular curves can normally be reduced to the much easier problem of discrete logarithm in a finite field through pairings.Practical binary elliptic curves thus have even order.
A non-supersingular binary curve can always be represented, through changes of variables, with a short Weierstraß equation: for two constants A and B.Moreover, the change of variable y → y + Cx, for any constant C, modifies the equation by replacing A with A + C 2 + C, but keeping B unchanged.Since C 2 + C can be any field element of trace zero, we can always replace A with any other constant with the same trace.In particular, when Tr(A) = 0, we can replace A with zero, and otherwise we can arrange for A to have minimal Hamming weight.The ten standard NIST binary curves (B-* and K-* curves with field degrees m ranging from 163 to 571 [NIS23]) systematically use an odd m, and constants A equal to zero or one.
The curve order is n = 2 t r for some integer t ≥ 1, and odd integer r.Moreover, the subgroup of 2 t -torsion points is always cyclic, and there are 2 t−1 generators; we call T one of these generators (arbitrarily chosen).There is a single point of order two, which is N = 2 t−1 T = (0, √ B).When Tr(A) = 1, t = 1, i.e. the curve has order 2r (this is then a double-odd curve).
We can apply the methodology exposed in section 3. The following notes apply: • A change of variable is applied, to move N to (0, 0).Namely, we add √ B to the y coordinate, and we get the equation: In this new equation, adding N to point (x, y) yields: • Point validation involves, as in section 3.3, checking whether a given point can be halved t times.A point (x, y) can be halved if and only if Tr(x) = Tr(A); computing the half-point involves solving a quadratic equation, which is feasible at moderate cost (e.g. with the halftrace when m is odd); see [Knu99] for details.In general, given a point P = (x, y), we can not only check whether it is an r-torsion point, but we can even compute the unique point Q and integer k such that Q ∈ E[r], 0 ≤ k < 2 t , and P = Q + kT .When Tr(A) = 1, this is especially easy since T = N and k = 1 − Tr(x).
Point validation can therefore be extended into full-curve computations by splitting points into their r-torsion and 2 t -torsion parts, and performing computations over the latter with simple integer additions modulo 2 t .For a prime-order group abstraction this is not needed; we only need to ensure that a given point is equal to P + N for P an r-torsion point.
• As in section 3.4, we can compress a point P + N = (x, y) into the value w = y/x; decompression can be coupled with point validation, and requires one inversion and solving t quadratic equations.
• The use of P + N instead of P immediately yields unified formulas for point additions, but they are not especially fast, and some exceptional cases remain when the neutral point appears as operand or output.To solve both issues, we change the coordinate system; this will be described in the next section.

(x,s) Coordinates
We consider the field F q with q = 2 m , and the elliptic curve of equation y 2 + xy = x 3 + ax 2 + bx for two constants a and b (with b ̸ = 0).The curve has order 2 t r for an odd integer r (the curve is usually chosen so that r is prime, but this is not strictly necessary here).An r-torsion point P is represented by the point P + N , with N = (0, 0) being the unique point of order two on the curve.For such a point P + N = (x, y), we define the extra coordinate s: s = y + x 2 + ax + b Given x, s can always be computed from y and vice versa; we can thus use (x, s) coordinates without any loss of information.When x ̸ = 0, we have s = y 2 /x, and y/x = s/x (which can be used for point compression).
Applying the s coordinates to the unified formulas obtained by using our methodology yields the following formulas.Given P 1 + N = (x 1 , s 1 ) and P 2 + N = (x 2 , s 2 ), which represent r-torsion points P 1 and P 2 , respectively, the representant P 3 +N = (P 1 +P 2 )+N = (x 3 , s 3 ) is such that: These formulas were derived from unified formulas, but it can be verified that they also work for the neutral N , both as operand and as output.These formulas are thus complete (for correct representants of r-torsion points).
In (x, s) coordinates, the neutral is N = (0, b); the opposite of P + N = (x, s) is −P + N = (x, s + x).A practical implementation may represent a point P + N = (x, s) with a quadruplet (X:S:Z:T ), such that Z ̸ = 0, x = √ bX/Z, s = √ bS/Z 2 and T = XZ (the factors " √ b" are used here to reduce the number of multiplications by √ b in the formulas).
Since our methodology applies to pre-existing standard curves, and the NIST B-* curves use pseudorandom B constants that were not chosen so that multiplication by B is inexpensive, we also include such multiplications in cost estimates; specifically, we denote by m b the cost of multiplying by With our chosen representation of points as quadruplets (X:S:Z:T ), we obtain generic point addition complete formulas with cost 8M + 2S + 2m b (when a = 0, this is reduced to 7M + 2S + 2m b ); this is better than the previously best known formulas, which were furthermore incomplete (using (x, λ) coordinates, with cost 11M + 2S [OLAR13]).We can also compute a sequence of n point doublings in cost n(2M + 4S + 2m b ) + 2S + 6m b , or in cost n(3M + 4S + m b ) + 3m b (which one is faster depends on the implementation platform, the constant b, and the number of successive doublings n).All formulas are detailed in appendix E.

Application to GLS254
We applied the (x, s) formulas to the GLS254 curve.That curve was defined by Aardal and Aranha [AA22], following up on previous work by Oliveira et al [OLAR14], and leveraging an extensive security analysis by Hankerson, Karabina and Menezes [HKM09], showing that the chosen parameters are not susceptible to the GHS attack [GHS02] even though the curve uses a composite extension degree m = 254.The GLS254 curve uses the GLS structure [GLS09], which offers an efficient endomorphism that can be used to optimize some operations such as multiplication of a point by a scalar.
The base field F q , with q = 2 254 , is defined as a tower of field extensions: • A field of order 2 127 is obtained by considering binary polynomials in z, modulo M = z 127 + z 63 + 1.
• Elements of F 2 254 are x = x 0 + ux 1 , where x 0 and x 1 are elements of F 2 127 , and u is a formal value such that u 2 + u = 1.
The trace (in F 2 254 ) of x 0 + ux 1 is equal to the trace of x 1 (in F 2 127 ), which is itself equal to the least significant bit of x 1 (its degree-zero coefficient as a binary polynomial of degree at most 126).Inversions in F 2 254 can be reduced to inversions in F 2 127 by noticing that (x 0 + ux 1 )(x 0 + x 1 + ux 1 ) = (x 2 0 + x 0 x 1 + x 2 1 ), which is an element of F 2 127 .A quadratic equation in F 2 254 can similarly be solved by computing two halftraces in F 2 127 .
The curve equation parameters are A = u and B = 1 + z 27 .The value B was chosen so that multiplications by B are especially efficient.The curve order is 2r, with r ≈ 2 253 .It should be noted that Tr(A) = 1; point validation is thus especially simple since it suffices to check that x has trace 0 (to ensure that a point is a P + N representant of an r-torsion point P ), and no actual point halving is required.
For better performance, when converting to (x, s) coordinates, we first raise everything to the power 4. Raising to power 4 is a field automorphism; the curve algebraic structure is thus preserved.With this step, we define b = B 2 , and √ b = B = 1 + z 27 ; this operation thus ensures that multiplications by √ b are fast.The curve endomorphism ζ is the following.For point P represented by P + N = (x, s), the endomorphism output is point P ′ represented by P ′ + N = (x ′ , s ′ ) with: where ϕ(v) = v 2 127 (Frobenius endomorphism on the field as a degree-2 extension of F 2 127 ).In the quadruplet representation (X:S:Z:T ), this yields: which can be computed with only a few field additions, which are just bitwise XORs.Using these formulas, our implementation (fully constant-time, in Rust) achieves signature generation and verification in 18k and 27k cycles on our test Intel Core i5-8259U x86 processor, respectively (again using 48-byte Schnorr signatures).This is 3 to 4 times faster than the fastest Ed25519 implementations on that same platform.This code leverages the carryless multiplication opcode (pclmulqdq) that originally appeared on x86 CPUs along with the AES hardware implementation (because it helps with the GCM authenticated encryption mode).On platforms that do not have such an opcode, we can use integer multiplications and masking, as originally suggested by Knuth [Knu69] (exercise 4 of section 4.6, page 363, solution on page 536).On a RISC-V SiFive-U74 core, with fast 64-bit multiplications but no carryless multiplication, we sign and verify in 378k and 636k cycles with GLS254, compared to 158k and 304k for Ed25519 on the same platform; the GLS254 performance is still sufficient for most uses of signatures, but not a record breaker.

Conclusion
We presented a generic method to avoid most or all of exceptional cases in point addition formulas on large classes of elliptic curves.The method can be applied to any curve with even order, on any kind of finite field (including binary fields); it offers secure prime-order group abstractions free of cofactor issues and amenable to constant-time implementations.The formulas are also as fast as, or faster than, existing fast curves, so that the security of formula completeness and removal of cofactors does not come at the price of reduced performance.
Existing standard binary curves (e.g. the NIST B-* and K-* curves) can immediately benefit from our formulas.For curves in odd characteristic fields, many standard curves have a prime order (e.g.P-256 or secp256k1) and our methods cannot be used with them; curves meant to be used in Montgomery or Edwards form (e.g.Curve25519) are compatible with our methodology, though they are not optimal in the sense that point validation can be expensive.In odd characteristic fields, the optimal case is double-odd curves, which is why we defined two such curves (jq255e and jq255s).
All our treatment considered the case of a curve subgroup of prime order r; however, our formulas also work with any odd composite r.A potential application is CSIDH, which, as defined in [CLM + 18], works in a finite field F q with q = 3 mod 8, over which (supersingular) curves of order q + 1 are defined.For such curves, r = (q + 1)/4 is a product of many small primes.Moreover, the suggested implementation of CSIDH already uses a curve equation with form y 2 = x 3 + ax 2 + x, which matches our equation from section 3.2 (using b = 1).Whether our methodology can help with speed and/or convenience of CSIDH implementations remains to be explored.

A Conventions for Formulas
In subsequent sections, detailed formulas are presented for operations on elliptic curves.The formulas use pseudocode in Python/Sage syntax; the code can actually run, provided that an appropriate finite field implementation is included under the name Fq: • The function Fq, when applied to an integer, must return an object representing that integer reduced modulo the field characteristic.Field elements must implement the usual arithmetic operators (addition +, subtraction -, multiplication *, division /, exponentiation **, and comparison operators == and !=).Field elements must also provide the methods is_zero() (which returns True if and only if the element is zero), is_square() (which returns True if and only if the element is a square in the field), and sqrt() (which returns a square root for the element, assuming one exists).
• Division by 2 (halving) is assumed to be fast, and denoted with the division operator.Practical implementations will use a dedicated function (for prime fields, this is typically done with a right-shift, followed by adding (q + 1)/2 if the dropped bit had value 1).
• The curve parameters are defined by the field elements a and b, corresponding to constants a and b, respectively.Some functions use extra constants derived from a and b, e.g.b ′ = a 2 − 4b, under the name bp.Each function header comment recalls the needed constants (if any).
• The sign() function must return the conventional "sign" (as defined in section 4.2) for a field element, as an integer of value 0 or 1.
For a prime field of modulus q, the Sage Zmod primitive is appropriate.For instance, for the curve jq255s, the following definitions provide the field and constants: (Constants b ′ , α, β and γ are all used by some of the functions listed in the following sections; for this curve, they all exist and the relevant functions can be used.)Curve points are represented by tuples of field elements, for parameters to function calls, and values returned by functions.The pseudocode aims at clarity of exposition and does not try to achieve constant-time operations; however, conditional jumps are used only in compression-related functions, and can be implemented with constant-time conditional setters at negligible computational overhead.
Appendices B, C and D use these conventions.They do not necessarily require the curve to be double-odd, unless explicitly specified.Appendix E covers the case of binary curves and uses its own conventions.

B Formulas for Jacobian (x,w) Coordinates
We use the conventions expressed in appendix A; this section is for a field of characteristic 3 or more.
In (x, w) coordinates, the y coordinate of a point (x, y) is replaced with w = y/x.We use a three-value representation (X:W :J) with J ̸ = 0, x = X/J 2 and w = W/J (it is called "Jacobian" because the J factor is equivalent to curve isomorphisms, as in classic Jacobian coordinates).The point-at-infinity O and the point of order two N do not have a well-defined w coordinate, but they can be represented as triplets (W 2 :W :0) and (0:W :0) for any W ̸ = 0, respectively.We present here four doubling formulas, presented as pseudocode (Python/Sage syntax); one computes 2P for an input point P , while the three others compute 2P + N for the input P .
The expressions given in table 1 work for all points P and P + N for any P ∈ E[r]; they can fail in the presence of a point of order two which is not N .If there is only one point of order two in the curve (i.e. a 2 − 4b / ∈ QR), then the formulas are complete; this includes the case of double-odd curves.The derived formulas will be used only on points on which exceptional cases are not encountered.For any such point P : Combining the expressions and seeking some computational shortcuts yields the functions presented thereafter.The function double_addN_xwj_generic() computes 2P + N for an input point P , on all curves (figure 1).If 4b − a 2 is a square, then we can define the constant γ = ( √ 4b − a 2 )/2 and use it in a slightly faster formula, with cost 1M + 6S (instead of 2M + 5S).This is exposed in function double_addN_xwj_gamma() (figure 2).In the expected case of a double-odd curve, for which a 2 − 4b / ∈ QR, γ exists if and only if q = 3 mod 4. When a = −1 and b = 1/2, some computational shortcuts can apply, resulting in a doubling process with cost 2M + 4S.This is illustrated in function double_addN_xwj_s() (figure 3).This case applies to the double-odd curve jq255s (defined in section 4.5).When a = 0, the curve has equation y 2 = x 3 + bx and doubling formulas with cost 1M + 5S apply; this is function double_xwj_e() (figure 4).Take care that for this function, the point 2P is returned, not 2P + N .

C Formulas for Fractional (x,u) Coordinates
In fractional (x, u) coordinates, an r-torsion point P is represented by the point P + N , whose coordinates (x, u) are themselves represented by the quadruplet (X:Z:U :T ) for Z ̸ = 0, T ̸ = 0, x = X/Z and u = U/T .The neutral element uses X = U = 0. Points are serialized by compressing them into w = 1/u = T /U (if U = 0 then w is conventionally set to 0).The decompression process rebuilds a matching (X:Z:U :T ) representation, taking care to verify that the obtained point is indeed a correct P + N representant of an r-torsion point P .Functions compress_xu() and decompress_doubleodd_xu() (figure 5) implement this process in the specific case of a double-odd curve, where point validation is easy.return False E = D.sqrt() x = (C + E)/2 if x.is_square(): x = x -E return (x, Fq(1), Fq(1), w) For general point addition, two extra constants α = (4b − a 2 )/(2b − a) and β = (a − 2)/(2b − a) are used.If the curve is such that a = 2b, then α and β are not defined; a workaround is to move to the isomorphic curve y 2 = x 3 + (aϵ 2 )x 2 + (bϵ 4 )x by multiplying x and u by ϵ 2 and 1/ϵ, respectively, for any ϵ ̸ = 0 (e.g.ϵ = 2).The function add_xu() (figure 7) computes (P 1 + P 2 ) + N , given P 1 + N and P 2 + N , assuming that both inputs are correct representants of r-torsion points.Negation is obtained by negating the u coordinate (figure 8).Subtraction is obtained by combining addition with negation of the second operand.A sequence of n point doublings can always be computed by calling function double_xu repeatedly, but a faster method, especially for long sequences, is to arrange for the first doubling to produce an output in Jacobian (x, w) coordinates, then performing (n − 2) doublings in that coordinate system (using the functions described in appendix B, which are faster), and making the last doubling with formulas that output the final result in fractional (x, u) coordinates.Function xdouble_xu (figure 10) implements this process.It internally uses a doubling function in Jacobian (x, w) coordinates, which can be any of the four functions from appendix B that applies to the used curve (in particular, both functions that return 2P and functions that return 2P + N can be used, since this is for the "internal" doublings).For some curves, the per-sequence overhead can be reduced with some curve-specific optimizations.If a = −1 and b = 1/2 (as in the jq255s curve), function xdouble_xu_s() (figure 11) saves one squaring (the saving is in the final doubling).If a = 0 (e.g. for the jq255e curve), then two squarings can be avoided with a slightly different process: only the ψ 1 isogeny is applied (i.e."half" of the first doubling), then n − 1 doublings in the dual curve E(−2a, a 2 − 4b) (for which −2a = 0), and finally θ ′ 1/2 is applied with conversion back to fractional (x, u) coordinates.This is shown in function xdouble_xu_e() (figure 12).

D Formulas for Extended (e,u) Coordinates
Extended (e, u) coordinates use the Jacobi quartic form.As explained in section 4, this really applies only to double-odd curves, for which the efficient compression and decompression process yields a proper prime-order group abstraction.We therefore only consider double-odd curves in this section.
In (e, u) coordinates, an r-torsion point P can be represented by either P or P + N , and both representants work equally well.The coordinates are represented by a quadruplet (E:Z:U :T ) such that Z ̸ = 0, e = E/Z, u = U/Z and U 2 = T Z.For all points, E ̸ = 0.The neutral uses U = T = 0 and E = ±Z.Compression and decompression use functions compress_doubleodd_eu() and decompress_doubleodd_eu(), respectively (figure 13), using the method exposed in section 4.2.Negation is done by function neg_eu() (figure 16); a subtraction can then be computed as the combination of a negation and an addition.Function xdouble_eu_generic() (figure 17) performs a sequence of n successive point doublings, and works on all double-odd curves; by setting n = 1, a single doubling is obtained.This process is simpler than in (x, u) coordinates because each point has two representants and we can use either, including for the output, which avoids the need for a special process for the last doubling.This is also why a dedicated function for single doublings is not needed.Function xdouble_eu_generic() performs the first doubling with a special sequence that also converts its output to Jacobian (x, w) coordinates, and the remaining n − 1 doublings use the functions specified in appendix B. For many curves, doublings in Jacobian (x, w) coordinates can be performed more efficiently than the generic 2M+5S procedure, and similar optimizations are applicable to the first doubling.Function xdouble_eu_gamma() (figure 18) handles curves such that −b ′ = 4b − a 2 ∈ QR (for a double-odd curve, b ′ / ∈ QR, hence this condition is equivalent to working in field F q with q = 3 mod 4); function xdouble_eu_s() (figure 19) handles the specific case of (a, b) = (−1, 1/2), which applies to curve jq255s; function xdouble_eu_e() (figure 20) is for curves with a = 0, e.g.jq255e.Function equals_xs() (figure 22) compares two points together; it uses the fact that s/x is the slope of the line from N to P + N and thus uniquely identifies a correct representant of an r-torsion point.The formula completeness can be verified by noticing that there is a single valid point P + N such that x = 0 (this is N itself), and also at most one valid point P + N such that s = 0 (if P + N is such a point, then the other one is −P , which is not a valid representant of an r-torsion point).Point addition is performed by function add_xs() (figure 23).The cost is 8M+2S+2m b .If the curve constant a is zero, then one multiplication is saved (the value T1T2 does not need to be computed) and the cost is lowered to 7M + 2S + 2m b .As was pointed out in section 5.2, if Tr(a) = 0, then we can always apply a curve isomorphism that sets a to zero while leaving b unchanged. .In this coordinate system, efficient formulas with cost 2M+4S+2m b are applied.Since these formulas are applicable to any non-supersingular binary curve in short Weierstraß format, they could be used outside of the (x, s) formalism; they are complete on the whole curve, provided that the point-at-infinity (O) is represented as (X:0:0:0) for some X ̸ = 0. Conversions to the short Weierstraß curve and back add some per-sequence overhead, for a total cost of n(2M + 4S + 2m b ) + 2S + 6m b .

Figure 5 :Figure 6 :
Figure 5: Point compression and decompression in fractional (x, u) coordinates for a double-odd curve.

#Figure 13 :Figure 14 :Figure 15 :
Figure 13: Point compression and decompression in extended (e, u) coordinates for a double-odd curve.

#Figure 16 :
Figure 16: Point negation in extended (e, u) coordinates for a double-odd curve.

#Figure 17 :
Figure 17: Sequence of point doublings in extended (e, u) coordinates for a double-odd curve.