Elliptic Curves 101

Building on Modular Math, we can now tackle Elliptic Curves. Elliptic curves are really nice because they allow for many of the same security properties as other systems, but require less storage.

Pre-Requisites

Fortunately, the first two articles of Elliptic Curve Cryptography: a gentle introduction - Andrea Corbellini are an exceptionally good primer; better than I would have created. So, rather than trying to re-write that material (and doing a worse job of it), just read the above, and then come back here for a worked example with a toy implementation.

Toy Implementation

The Weierstrass Normal Form of Elliptic curves is

\(y^2 = (x^3 + a*x + b) \bmod{p}\)

Bitcoin uses secp256k1, which looks like

\(y^2 = (x^3 + 7) \bmod{p}\)

using a huge prime number¹ as the modulus. For our toy implementation, we're going to use a much smaller prime number: 97. A point [x, y] is on the curve if it satisfies our equation. For instance, here is [1, 28].

\(\begin{align*} y^2 &= (x^3 + 7) \bmod{p} \\ 28^2 \bmod 97 &= 1 + 7 \\ 8 &= 8 \end{align*}\)

After that, we need to code our group operation, which we're going to call multiply.

Multiplication

It's relatively straightforward and copies the math from Andrea Corbellini.

export interface ECPoint {
  x: number
  y: number
}

export const curveMod = 97
export const curveA = 0
export const curveB = 7

// y^2 = x^3 + curveA * x + curveB mod curveMod

export const identity: ECPoint = { x: -1, y: -1 }

export function multiply(a: ECPoint, b: ECPoint): ECPoint {
  if (a.x === identity.x) {
    return b
  }

  if (b.x === identity.x) {
    return a
  }

  if (a.x === b.x && a.y !== b.y) {
    return identity
  }

  let slope = 0

  if (a.x === b.x) {
    slope = (3 * Math.pow(a.x, 2)) * modInverse(2 * a.y, curveMod)
  } else {
    slope = (b.y - a.y) * modInverse(b.x - a.x, curveMod)
  }

  const x = mod(Math.pow(slope, 2) - a.x - b.x, curveMod)
  const y = mod(slope*(a.x - x) - a.y, curveMod)

  return { x, y }
}

The preamble represents the identify laws:

• A • identity = A

• identity • A = A

• A • -A = identity

Then, we need to know if we are squaring a point (A • A) or multiplying two distinct points (A • B), since that affects the slope.

If we are squaring a point, the slope will be the derivative:

\(\begin{align*} y^2 &= x^3 + 7 \\ 2y * \dfrac{dy}{dx} &= 3x^2 \\ \dfrac{dy}{dx} &= \dfrac{3x^2}{2y} \end{align*}\)

It works the same way in modular arithmetic except instead of dividing, we

multiply by the modular inverse. Our slope is

slope = 3x^2 = modInverse(2y, 97)

For two distinct points, we can calculate the slope regularly.

slope = (b.y - a.y) / (b.x - a.x)

Division is off the table, so we use `modInverse` again:

slope = (b.y - a.y) * modInverse(b.x - a.x, 97)

Once we know the slope, we can write the equation for the line:

y = slope * x + b

Which must hold true for both of our points. We can substitute to solve for `b`

a.y = slope * a.x + b
b   = a.y - slope * a.x

So now we have

y = slope * x + (a.y - slope * a.x)            // line equation
0 = x^3 + 7 - y^2                              // curve equation
0 = x^3 + 7 - (slope*x + b)^2                  // substitute out y
0 = x^3 + 7 - ((slope*x)^2 + 2b*slope*x + b^2) // distribute
0 = x^3 + 7 - (slope*x)^2 - 2b*slope*x - b^2
0 = x^3 - slope^2*x^2 - 2b*slope*x + (7 - b^2)

This gives us our equation in standard polynomial form and lets us find the roots. We know that there are three roots: our known `a.x` and`b.x` and the solution we’re looking for.

We can write a cubic polynomial in roots form like this:

(x - root1) * (x - root2) * (x - root3) = 0

Since we already know two our our roots, we can substitute:

(x - a.x)*(x - b.x)*(x - solution) = 0
(x^2 - b.x*x - a.x*x + a.x*b.x)*(x - solution) = 0
(x^2 - (b.x + a.x)*x + a.x*b.x)*(x - solution) = 0
x^3 - solution*x^2 - (b.x + a.x)*x^2 + (b.x + a.x)*solution*x + a.x*b.x*x - a.x*b.x*solution = 0
x^3 - (solution + b.x + a.x)*x^2 + ((b.x + a.x)*solution + a.x*b.x)*x - a.x*b.x*solution = 0

Now we have two different representations of the same polynomial, both in standard form. Since they are the same polynomial, all of their coefficients must also be the same, which allows us to greatly simplify.

1                              = 1          // coefficient of x^3

solution + b.x + a.x           = slope^2    // coefficient of x^2

(b.x + a.x)*solution + a.x*b.x = -2b*slope  // coefficient of x

-a.x*b.x*solution              = 7 - b^2    // the free term

We can do some algebra:

solution.x = slope^2 - b.x - a.x                  // coefficient of x^2

solution.x = (-2b*slope - a.x*b.x) / (b.x + a.x)  // coefficient of x

solution.x = (7 - b^2) / (-a.x*b.x)               // the free term

Of these, the x^2 term is by far the easiest to calculate. That means for any multiplying any two points a and b, the x-value of the third point is

slope^2 - a.x - b.x

For modular arithmetic, we just need to take the mod at the end:

solution.x = mod(Math.pow(slope, 2) - a.x - b.x, curveMod)

Then we can substitute to get y with the line equation:

solution.y = slope * x + (a.y - slope * a.x) // line equation
solution.y = slope * solution.x + a.y - slope * a.x
solution.y = slope*(solution.x - a.x) + a.y

Finally, we need to apply the third step: translating it across the y-axis. With the real numbers, we would multiply by -1 and call it a day. In modular arithmetic we also need to take the mod at the end.

final.y = mod(-1 * slope*(solution.x - a.x) - a.y, 97)
final.y = mod(slope*(a.x - solution.x) - a.y, 97)

That was the last piece to understand our (in hindsight) compact algorithm:

if (a.x === b.x) {
  slope = (3 * Math.pow(a.x, 2)) * modInverse(2 * a.y, curveMod)
} else {
  slope = (b.y - a.y) * modInverse(b.x - a.x, curveMod)
}

const x = mod(Math.pow(slope, 2) - a.x - b.x, curveMod)
const y = mod(slope*(a.x - x) - a.y, curveMod)

return { x, y }

Exponentiation

We also want to be able to call

exponentiate(a: ECPoint, k: number): ECPoint

in order to represent

a*a*a*a...

k times. Here's the simple, brute-force approach:

function exponentiate(a: ECPoint, k: number): ECPoint {
  let result = a
  for (let i = 0; i < k - 1; i++) {
    result = multiply(result, a)
  }

  return result
}

It just works! It is also miserably slow. Here is what the execution log looks like:

result = a^1
result = a^1 * a = a^2
result = a^2 * a = a^3
...
result = a^(k-1) * a = a^k

Much better is to take advantage of doubling. Say that we wanted a^55. We can express that as

a^55 = a^(1 + 2 + 4 + 16 + 32)
     = a * a^2 * a^4 * a^16 * a^32

Here's an execution log that allows us to do this:

let result = 1
let double = a
result = result * double = a
double = double * double = a * a = a^2
result = result * double = a * a^2 = a^3
double = double * double = a^2 * a^2 = a^4
result = result * double = a^3 * a^4 = a^7
double = double * double = a^4 * a^4 = a^8
double = double * double = a^8 * a^8 = a^16
result = result * double = a^7 * a^16 = a^23
double = double * double = a^16 * a^16 = a^32
result = result * double = a^23 * a^32 = a^55

This would have previously taken 54 multiplication operations. Now, it is taking 10. Abstractly, the brute-force approach takes k-1 operations and the doubling strategy takes 2*log_2(k) in the worst case.

In order to make doubling work, we need a way to factor 55 into powers of 2. Here is one such procedure:

55 = 1 + 54
   = 1 + 2 * 27
   = 1 + 2 * (1 + 26)
   = 1 + 2 + 2*26
   = 1 + 2 + 2(2*13)
   = 1 + 2 + 4*13
   = 1 + 2 + 4*(1 + 12)
   = 1 + 2 + 4 + 4*12
   = 1 + 2 + 4 + 4*(2*6)
   = 1 + 2 + 4 + 8*6
   = 1 + 2 + 4 + 8*(2*3)
   = 1 + 2 + 4 + 16*3
   = 1 + 2 + 4 + 16*(1+2)
   = 1 + 2 + 4 + 16 + 16*2
   = 1 + 2 + 4 + 16 + 32

Whenever we encounter an *odd* number, we use that term, since odd numbers can always be represented as “1 + 2*n”, and the “1 + ...” part means that when we distribute our power of 2, it sticks around.

Here is how we can turn this into an algorithm:

function exponentiate(a: ECPoint, k: number): ECPoint {
  let result = identity
  let double = a

  while (k > 0) {
    if (mod(k, 2) === 1) {
      result = multiply(result, double)
    }

    double = multiply(double, double)
    k = Math.floor(k / 2)
  }

  return result
}

Generators and Trap Door Functions

Here's the idea: say that we all publicly agree to use {x: 1, y: 28} as our base point. I calculate

exponentiate({x: 1, y: 28}, k)

for some secret value k and get {x: 5, y: 61}. Are you able to tell what k I used? You could keep guessing every possible k until you found the answer. There are more sophisticated approaches, but none that I am aware of that get you below linear time with respect to k. If there are 2 digits of possible K values, you will probably have to guess and check check 2 digits of possible answers.

If there are, instead, 100 digits of possible `k` values, you will be checking for a long time.

Meanwhile, if I tell you that my my secret k value was 35, it is quick to verify that I am honest. This is the foundation of a Trapdoor Function - it is easy to go through in one direction, but difficult in the other.

The concept of group theory Generators is a little too abstract for our purposes, but in practical terms, we can pick a base point like {x: 1, y: 28} and keep applying the operation to itself and if we eventually wind up back at itself, we have a cyclic group - that's the case here:

exponentiate({x: 1, y: 28}, 80) === {x: 1, y: 28}

This means

exponentiate({x: 1, y: 28}, 85) === exponentiate({x: 1, y: 28}, 6)

since we keep looping back.

That is all we're looking for in a generator in this context, so we can code it in a few lines:

const generator: ECPoint = {x: 1, y: 28}

function g(n: number): ECPoint {
  return exponentiate(generator, n)
}

We could have picked whichever point we wanted, so long as it satisfies our curve equation. I just set x = 1 and then solved for y via guess-and-check. To do this in a reasonable amount of time with larger numbers, we can use the [Tonelli-Shanks algorithm.

1. 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F