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() and CSWAP().

The CMOVE(x, y, condition) operation is a conditional move. It means “if condition, assign y to x”:

if condition
    x = y

The CSWAP(x, y, condition) operation is a conditional swap. It means “if condition, swap the values of x and y”:

if condition
    x, y = y, x

CMOVE() and 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

The 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 inv_sqrt(x) if 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

Where sqrt_m1 = √−1 = abs((2^{(q−1) / 4)}). For Curve25519, p = 2²⁵⁵−19 and sqrt_m1 = 19681161376707505956807079304988542015446066515923890162744021073123829784752

Recommended way to implement inv_sqrt(x) if 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

The direct map

Compute the point (u, v) from the representative r:

w = −A / (1 + Z r²)
e = legendre(w³ + A w² + B w)
u = e w − (1 − e) (A / 2)
v = −e √(u³ + A u² + 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 Z_u and Z_v constants derived from Z. In practice, Z is chosen in a way that minimises computation: either a small number, or something that eliminates Z_u and Z_v.

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

Curve25519 parameters

`q`	Prime field	2²⁵⁵ − 19
`A`	Curve constant	486662
`B`	Curve constant	1
`Z`	Non-square	2
`Z_u`	−`Z`√−1	18533721865243085798171333894366869895742859300972501694240749857708905250445
`Z_v`	√`Z_u`	19681161376707505956807079304988542015446066515923890162744021073123829784751

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 Z_u = Z_v = 1.

Curve448 parameters

`q`	Prime field	2⁴⁴⁸ − 2²²⁴ − 1
`A`	Curve constant	156326
`B`	Curve constant	1
`Z`	Non-square	−1
`Z_u`	−`Z`	1
`Z_v`	√`Z_u`	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

Curve25519

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

Curve448

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)

where sqrt_d2 = 2√−39081 = 197888467295464439538354009753858038256835152591059802148199779196087404232002515713604263127793030747855424464185691766453844835192428.

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

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)

Where x = 14399317868200118260347934320527232580618823971194345261214217575416788799818 and y = 2707385501144840649318225287225658788936804267575313519463743609750303402022.

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)

Curve448

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