Elliptic curve point multiplication

Point at infinity ${\displaystyle {\mathcal {O}}}$ is the identity element of elliptic curve arithmetic. Adding it to any point results in that other point, including adding point at infinity to itself. That is:

${\displaystyle {\begin{aligned}{\mathcal {O}}+{\mathcal {O}}={\mathcal {O}}\\{\mathcal {O}}+P=P\end{aligned}}}$

Point at infinity is also written as 0.

Point negation is finding such a point, that adding it to itself will result in point at infinity (${\displaystyle {\mathcal {O}}}$).

${\displaystyle {\begin{aligned}P+(-P)={\mathcal {O}}\end{aligned}}}$

For elliptic curves of the form E: y² = x³ + ax + b, negation is a point with the same x coordinate but negated y coordinate:

${\displaystyle {\begin{aligned}(x,y)+(-(x,y))&={\mathcal {O}}\\(x,y)+(x,-y)&={\mathcal {O}}\\(x,-y)&=-(x,y)\end{aligned}}}$

With 2 distinct points, P and Q, addition is defined as the negation of the point resulting from the intersection of the curve, E, and the straight line defined by the points P and Q, giving the point, R.[1]

${\displaystyle {\begin{aligned}P+Q&=R\\(x_{p},y_{p})+(x_{q},y_{q})&=(x_{r},y_{r})\end{aligned}}}$

Assuming the elliptic curve, E, is given by y² = x³ + ax + b, this can be calculated as:

${\displaystyle {\begin{aligned}\lambda &={\frac {y_{q}-y_{p}}{x_{q}-x_{p}}}\\x_{r}&=\lambda ^{2}-x_{p}-x_{q}\\y_{r}&=\lambda (x_{p}-x_{r})-y_{p}\\\end{aligned}}}$

These equations are correct when neither point is the point at infinity, ${\displaystyle {\mathcal {O}}}$, and if the points have different x coordinates (they're not mutual inverses). This is important for the ECDSA verification algorithm where the hash value could be zero.

Where the points P and Q are coincident (at the same coordinates), addition is similar, except that there is no well-defined straight line through P, so the operation is closed using a limiting case, the tangent to the curve, E, at P.

This is calculated as above, taking derivatives (dE/dx)/(dE/dy):[2]

${\displaystyle \lambda ={\frac {3x_{p}^{2}+a}{2y_{p}}}}$

where a is from the defining equation of the curve, E, above.

The simplest method is the double-and-add method,[3] similar to square-and-multiply in modular exponentiation. The algorithm works as follows:

To compute sP, start with the binary representation for s: ${\displaystyle s=s_{0}+2s_{1}+2^{2}s_{2}+\cdots +2^{n-1}s_{n-1}}$, where ${\displaystyle s_{0}~..~s_{n-1}\in \{0,1\},n=\lceil \log _{2}{s}\rceil }$.

• Iterative algorithm, index increasing:

   let bits = bit_representation(s) # the vector of bits (from LSB to MSB) representing s
   let res = ${\displaystyle {\begin{aligned}{\mathcal {O}}\end{aligned}}}$ # point at infinity
   let temp = P # track doubled P val
   for bit in bits: 
       if bit == 1:            
           res = res + temp # point add
       temp = temp + temp # double
   return res

• Iterative algorithm, index decreasing:

   let bits = bit_representation(s) # the vector of bits (from LSB to MSB) representing s
   let i = length(bits) - 2
   let res = P
   while (i >= 0): # traversing from second MSB to LSB
       res = res + res # double
       if bits[i] == 1:            
           res = res + P # add
       i = i - 1
   return res

Note that both of the iterative methods above are vulnerable to timing analysis. See Montgomery Ladder below for an alternative approach.

• Recursive algorithm:

  f(P, d) is
     if d = 0 then
         return 0                         # computation complete
     else if d = 1 then
         return P
     else if d mod 2 = 1 then
         return point_add(P, f(P, d - 1)) # addition when d is odd
     else
         return f(point_double(P), d / 2)   # doubling when d is even

where f is the function for multiplying, P is the coordinate to multiply, d is the number of times to add the coordinate to itself. Example: 100P can be written as 2(2[P + 2(2[2(P + 2P)])]) and thus requires six point double operations and two point addition operations. 100P would be equal to f(P, 100).

This algorithm requires log₂(d) iterations of point doubling and addition to compute the full point multiplication. There are many variations of this algorithm such as using a window, sliding window, NAF, NAF-w, vector chains, and Montgomery ladder.

In the windowed version of this algorithm,[3] one selects a window size w and computes all ${\displaystyle 2^{w}}$ values of ${\displaystyle dP}$ for ${\displaystyle d=0,1,2,\dots ,2^{w}-1}$. The algorithm now uses the representation ${\displaystyle d=d_{0}+2^{w}d_{1}+2^{2w}d_{2}+\cdots +2^{mw}d_{m}}$ and becomes

  Q ← 0
  for i from m to 0 do
      Q ← point_double_repeat(Q, w)
      if di > 0 then
          Q ← point_add(Q, diP) # using pre-computed value of diP
  return Q

This algorithm has the same complexity as the double-and-add approach with the benefit of using fewer point additions (which in practice are slower than doubling). Typically, the value of w is chosen to be fairly small making the pre-computation stage a trivial component of the algorithm. For the NIST recommended curves, ${\displaystyle w=4}$ is usually the best selection. The entire complexity for a n-bit number is measured as ${\displaystyle n+1}$ point doubles and ${\displaystyle 2^{w}-2+{\tfrac {n}{w}}}$ point additions.

In the sliding-window version, we look to trade off point additions for point doubles. We compute a similar table as in the windowed version except we only compute the points ${\displaystyle dP}$ for ${\displaystyle d=2^{w-1},2^{w-1}+1,\dots ,2^{w}-1}$. Effectively, we are only computing the values for which the most significant bit of the window is set. The algorithm then uses the original double-and-add representation of ${\displaystyle d=d_{0}+2d_{1}+2^{2}d_{2}+\cdots +2^{m}d_{m}}$.

  Q ← 0
  for i from m downto 0 do
      if di = 0 then
          Q ← point_double(Q)
      else 
          t ← extract j (up to w − 1) additional bits from d (including di)
          i ← i − j
          if j < w then
              Perform double-and-add using t 
              return Q
          else 
              Q ← point_double_repeat(Q, w)
              Q ← point_add(Q, tP)
  return Q

This algorithm has the benefit that the pre-computation stage is roughly half as complex as the normal windowed method while also trading slower point additions for point doublings. In effect, there is little reason to use the windowed method over this approach, except that the former can be implemented in constant time. The algorithm requires ${\displaystyle w-1+n}$ point doubles and at most ${\displaystyle 2^{w-1}-1+{\tfrac {n}{w}}}$ point additions.

In the non-adjacent form we aim to make use of the fact that point subtraction is just as easy as point addition to perform fewer (of either) as compared to a sliding-window method. The NAF of the multiplicand ${\displaystyle d}$ must be computed first with the following algorithm

   i ← 0
   while (d > 0) do
       if (d mod 2) = 1 then 
           di ← d mods 2w
           d ← d − di
       else
           di = 0
       d ← d/2
       i ← i + 1
   return (di−1, di-2, …, d0)

Where the signed modulo function mods is defined as

   if (d mod 2w) >= 2w−1
       return (d mod 2w) − 2w
   else
       return d mod 2w

This produces the NAF needed to now perform the multiplication. This algorithm requires the pre-computation of the points ${\displaystyle \lbrace 1,3,5,\dots ,2^{w-1}-1\rbrace P}$ and their negatives, where ${\displaystyle P}$ is the point to be multiplied. On typical Weierstrass curves, if ${\displaystyle P=\lbrace x,y\rbrace }$ then ${\displaystyle -P=\lbrace x,-y\rbrace }$. So in essence the negatives are cheap to compute. Next, the following algorithm computes the multiplication ${\displaystyle dP}$:

   Q ← 0
   for j ← i − 1 downto 0 do
       Q ← point_double(Q)
       if (dj != 0)
           Q ← point_add(Q, djP)
   return Q

The wNAF guarantees that on average there will be a density of ${\displaystyle {\tfrac {1}{w+1}}}$ point additions (slightly better than the unsigned window). It requires 1 point doubling and ${\displaystyle 2^{w-2}-1}$ point additions for precomputation. The algorithm then requires ${\displaystyle n}$ point doublings and ${\displaystyle {\tfrac {n}{w+1}}}$ point additions for the rest of the multiplication.

One property of the NAF is that we are guaranteed that every non-zero element ${\displaystyle d_{i}}$ is followed by at least ${\displaystyle w-1}$ additional zeroes. This is because the algorithm clears out the lower ${\displaystyle w}$ bits of ${\displaystyle d}$ with every subtraction of the output of the mods function. This observation can be used for several purposes. After every non-zero element the additional zeroes can be implied and do not need to be stored. Secondly, the multiple serial divisions by 2 can be replaced by a division by ${\displaystyle 2^{w}}$ after every non-zero ${\displaystyle d_{i}}$ element and divide by 2 after every zero.

It has been shown that through application of a FLUSH+RELOAD side-channel attack on OpenSSL, the full private key can be revealed after performing cache-timing against as few as 200 signatures performed.[4]

The Montgomery ladder[5] approach computes the point multiplication in a fixed amount of operations. This can be beneficial when timing, power consumption, or branch measurements are exposed to an attacker performing a side-channel attack. The algorithm uses the same representation as from double-and-add.

  R0 ← 0
  R1 ← P
  for i from m downto 0 do
      if di = 0 then
          R1 ← point_add(R0, R1)
          R0 ← point_double(R0)
      else
          R0 ← point_add(R0, R1)
          R1 ← point_double(R1)

      // invariant property to maintain correctness
      assert R1 == point_add(R0, P)
  return R0

This algorithm has in effect the same speed as the double-and-add approach except that it computes the same number of point additions and doubles regardless of the value of the multiplicand d. This means that at this level the algorithm does not leak any information through branches or power consumption.

However, it has been shown that through application of a FLUSH+RELOAD side-channel attack on OpenSSL, the full private key can be revealed after performing cache-timing against only one signature at a very low cost.[6]

Rust code for Montgomery Ladder:[7]

/// Constant operation point multiplication. 
/// NOTE: not memory safe.
/// * `s`: scalar value to multiply by
/// * multiplication is defined to be P₀ + P₁ + ... Pₛ
fn sec_mul(&mut self, s: big) -> E521{
    let mut r0 = get_e521_id_point();
    let mut r1 = self.clone();
    for i in (0..=s.significant_bits()).rev()  {
        if s.get_bit(i) {
            r0 = r0.add(&r1);
            r1 = r1.add(&r1.clone());
        } else {
            r1 = r0.add(&r1);
            r0 = r0.add(&r0.clone());
        }
    }
    r0 // r0 = P * s
}