Explicit Formulas for Elligator 2
Conditional move & swap
We recommend constant time
implementations to avoid leaking secrets.
So instead of writing regular conditionals we use
CMOVE(x, y, condition) operation is a conditional move.
It means “if
condition, assign y to
if condition x = y
CSWAP(x, y, condition) operation is a conditional swap.
It means “if
condition, swap the values of x and
if condition x, y = y, x
CSWAP() must be implemented in a way that doesn’t use
actual branches or indices,
to make sure they run in constant time.
Inverse square root
inv_sqrt(x) function implements inverse square root.
Given a non-zero square s and a non-square n:
inv_sqrt(0) = (0, True)
inv_sqrt(s) = (√(1/x), True)
inv_sqrt(n) = (√(k/x), False), where k = √−1 if q ≡ 5 (mod 8) (Curve25519), and k = −1 if q ≡ 3 (mod 4) (Curve448).
Recommended way to implement
q ≡ 5 (mod 8) (Curve25519):
inv_sqrt(x) ir = x^((p - 5) / 8) quartic = x × ir^2 m_sqrt_m1 = quartic == -1 or quartic == -sqrt_m1 is_square = quartic == -1 or quartic == 1 or x == 0 CMOVE(ir, ir × sqrt_m1, m_sqrt_m1) return ir, is_square
sqrt_m1 = √−1 = abs((2(q−1) / 4)).
p = 2255−19 and
sqrt_m1 = 196811613
Recommended way to implement
q ≡ 3 (mod 4) (Curve448):
inv_sqrt(x) ir = x^((p - 3) / 4) legendre = x × ir^2 is_square = legendre != -1 return ir, is_square
In such cases, start by recovering a square root; guidance for this can be found in I-D.draft-irtf-cfrg-hash-to-curve-16, appendix I. Once one of the two square roots has been determined, compute the multiplicative inverse thereof to get the inverse square root. This is not very optimised, but works. Optimisation for any individual case is left as an exercise for the reader.
The direct map
Compute the point (u, v) from the representative r:
- w = −A / (1 + Z r2)
- e =
legendre(w3 + A w2 + B w)
- u = e w − (1 − e) (A / 2)
- v = −e √(u3 + A u2 + B u)
- P = (u, v)
When B = 1 (this is true of all Montgomery curves compatible with Elligator 2), the recommended way to implement the above is:
direct_map(r) u = r^2 t1 = u × Z v = t1 + 1 t2 = v^2 t3 = A^2 t3 = t3 × t1 t3 = t3 - t2 t3 = t3 × A t1 = t2 × v t1, s? = inv_sqrt(t3 × t1) u = u × Zu v = r × Zv CMOVE(u, 1, s?) CMOVE(v, 1, s?) v = v × t3 v = v × t1 t1 = t1^2 u = u × -A u = u × t3 u = u × t2 u = u × t1 t1 = -v CMOVE(v, t1, is_square XOR is_negative(v)) P = (u, v) return P
Where A, is a curve constant, Z an arbitrary non-square parameter, and Zu and Zv constants derived from Z. In practice, Z is chosen in a way that minimises computation: either a small number, or something that eliminates Zu and Zv.
In the context of key exchange with X25519 and X448, v can be ignored entirely. This saves a couple operations:
direct_map(r) u = r^2 t1 = u × Z v = t1 + 1 t2 = v^2 t3 = A^2 t3 = t3 × t1 t3 = t3 - t2 t3 = t3 × A t1 = t2 × v t1, s? = inv_sqrt(t3 × t1) u = u × Zu CMOVE(u, 1, s?) t1 = t1^2 u = u × -A u = u × t3 u = u × t2 u = u × t1 return u
|q||Prime field||2255 − 19|
These are the values suggested in the original Elligator paper, and used in most of the implementations we know of. A good alternative for Z, used by ristretto255, is Z = √−1, such that Zu = Zv = 1.
|q||Prime field||2448 − 2224 − 1|
The inverse map
If the point P = (u, v) satisfies the following:
- u ≠ −A;
- −Zu(u + A) is a square;
- if v = 0, then u = 0 as well.
Then the representative of P is computed thus:
- r = √(−u / (Z (u + A))) when v is non-negative,
- r = √(−(u + A) / (Z u)) otherwise.
The recommended way to test the eligibility of a point (u, v) and compute its representative r is:
inverse_map((u, v)) t = u + A r = -Z × u r = r × t r, s? = inv_sqrt(r) if not s?, return nothing # no representative CMOVE(u, t, is_negative(v)) r = u × r t = -r CMOVE(r, t, r.is_negative()) return r
The constants A and Z are the same as the direct map.
Converting from Edwards to Montgomery
The conversion formulas are known, but to avoid extraneous divisions we generally take advantage of projective coordinates (X, Y, Z) such that:
- x = X/Z
- y = Y/Z
The twisted Edwards (x, y) are converted to the Montgomery u coordinate thus:
ed25519_to_curve25519(X, Y, Z) T = Z + Y u = Z - Y u = 1 / u u = u * T return u
This one requires two conversions: from Ed448Goldilocks to an isogenous Edwards curve, then from that Edwards curve to Curve448. The isogeny conversion is as follows:
isogeny_to_ed(X, Y, Z): X2 = X^2 Y2 = Y^2 du = (Z^2 * 2) - X2 - Y2 V = Y2 - X2 D = du * V U = V * X * Y * sqrt_d2 V = (Y2 + X2) * du return (U, V, D)
sqrt_d2 = 2√−39081 = 197888467
Once we’ve done the above then added the low order point, we can convert to Montgomery (this is slightly different from Curve25519):
ed_to_curve448(X, Y, Z) T = Y + Z u = Y - Z u = 1 / u u = u * T return u
Adding a low order point
Curve25519 and Curve448 have few low order points (8 and 4 respectively). This enables a couple optimisations.
Curve25519 doesn’t need scalar multiplication to select its low order point. We could instead select them with a pre-computed table and constant time lookup, but even that is overkill. Ed25519 has 4 generators of the low order group (4 points of order 8). Let’s (somewhat arbitrarily) chose the generator H = (x, y) where both x and y are positive (below q/2). A naive look up table would look like this:
0.H = ( 0 , 1) 1.H = ( x , y) 2.H = ( sqrt(-1), -0) 3.H = ( x , -y) 4.H = (-0 , -1) 5.H = (-x , -y) 6.H = (-sqrt(-1), 0) 7.H = (-x , y)
x = 143993178
We can instead reduce the pre-computed table to only 3 field elements: x, y, and √−1. First, we need a selection function:
select(x, k, i) r = 0 CMOVE(r, k, (i / 2) % 2 == 1) # bit 1 CMOVE(r, x, (i / 1) % 2 == 1) # bit 0 CMOVE(r, -r, (i / 4) % 2 == 1) # bit 2 return r
Finally, we can compute the scalar multiplication i.H by selecting the right low order point:
select_lop(i): lx = select(x, sqrt_m1, i ) ly = select(y, 1 , i+2) return (lx, ly)
Curve 448 has 4 low order points: (0, 1), (1, 0), (0, -1), (-1, 0). Among them (1, 0) and (-1, 0) have order 4. Let’s (somewhat arbitrarily) decide that H = (1, 0). Here we don’t even need a clever way to compute i.H, because:
- (x, y) + 0.H = (x, y)
- (x, y) + 1.H = (y, −x)
- (x, y) + 2.H = (−x, −y)
- (x, y) + 3.H = (−y, x)
All we have to do is conditionally negate and swap x and y:
add_lop(x, y, i): l = (i//1) % 2 == 1 # bit 0 h = (i//2) % 2 == 1 # bit 1 CSWAP(x, y, l) CMOVE(x, -x, h) CMOVE(y, -y, l XOR h) return x, y