The series consists of four parts; each part uses the concepts discovered in the previous parts:

The Magic of Elliptic Curves
The Magic of Elliptic Curves Combined with Finite Fields
Using The Magic of Elliptic Curves to Sign and Verify Messages
ECDSA Implementation in Python Using Zero Dependencies [you're here]

The bottlenecks:

We need to be able to perform basic arithmetical operations on very large numbers. In programming we can easily operate on large numbers, using bignum arithmetic. Just google it, and you will find out how easy it is. So our programming language must support it, or we should use some external package to work with it. In the examples of this part, I will use Python, which supports bignum arithmetic out of the box. For the Live Demo, I will use JavaScript, and there we will need the BigNumber.js package.
The other bottleneck that we will encounter is finding the multiplicative inverse of a very large number. Obviously, brute force is not going to work. The multiplicative inverse can be found by the Extended Euclidean algorithm, which has the complexity of O(log(n)). Python (3.8+) can do it out of the box with its built-in pow function:

def find_inverse(number, modulus):
    return pow(number, -1, modulus)

If you need the actual implementation of the algorithm, check my Live Demo!

Let’s start writing our code!

We need one simple thing, related to the elliptic curve: Point. Let’s define a class Point. In its constructor, we should make check, whether the point lies on the curve:

class Point:
    def __init__(self, x, y, curve_config):
        a = curve_config['a']
        b = curve_config['b']
        p = curve_config['p']

        if (y ** 2) % p != (x ** 3 + a * x + b) % p:
            raise Exception("The point is not on the curve")

        self.x = x
        self.y = y
        self.curve_config = curve_config

We need to be able to compare two points, add them together, and multiply them by an integer.

Let’s add a method to check if two points are equal:

def is_equal_to(self, point):
        return self.x == point.x & self.y == point.y

Now let’s implement add method, which returns a new Point as the result of addition:

def add(self, point):
        p = self.curve_config['p']

        if self.is_equal_to(point):
            slope = (3 * point.x ** 2) * find_inverse(2 * point.y, p) % p
        else:
            slope = (point.y - self.y) * find_inverse(point.x - self.x, p) % p

        x = (slope ** 2 - point.x - self.x) % p
        y = (slope * (self.x - x) - self.y) % p
        return Point(x, y, self.curve_config)

All the formulas are listed in Part 2.

Now let’s implement the multiply method:

The most straightforward implementation would be this:

def multiply(self, times):
    point = self
    for i in range(times - 1):
        point = point.add(self)
    return point

But let’s say we need to multiply our point by a big number: 115792089237316195. Even if we had the speed of 1 billion additions per second, this would take 3.6 years to calculate this point!

And this is not even a big number for us! Here is a big number:

115792089237316195423570985008687907852837564279074904382605163141518161494337

Calculating the point in this way would take billions of billions of billions of billions… of years!

We can define that the efficiency of this algorithm above is O(n), which is of no use for our purposes. If you remember, there is an easy way to achieve O(log2(n)) complexity by continuously doubling our point:

2P = P+P

4P = 2P + 2P

8P = 4P + 4P

16P = 8P + 8P

32P= 16P + 16P

64P = 32P + 32P

And so log2(115792089237316195) = 56

log2(115792089237316195423570985008687907852837564279074904382605163141518161494337) = 256

So we don’t need billions of billions of billions… of years. We just need 256 operations to get to this large point!

Just one moment: to efficiently multiply by values that are not a degree of 2, it’s reasonable to store all the previous values, and then combine the results together.

For example, if we need to get 100P, we can no longer double 64P. Neither we can add points one by one: potentially this would take billions of billions of years on larger numbers. What’s reasonable to do instead, is:

96P = 64P + 32P

100P = 96P + 4P

So for that purpose, we need to store all the previous P’s and afterwards efficiently use them.

So here is an efficient implementation:

def multiply(self, times):
        current_point = self
        current_coefficient = 1

        pervious_points = []
        while current_coefficient < times:
            # store current point as a previous point
            pervious_points.append((current_coefficient, current_point))
            # if we can multiply our current point by 2, do it
            if 2 * current_coefficient <= times:
                current_point = current_point.add(current_point)
                current_coefficient = 2 * current_coefficient
            # if we can't multiply our current point by 2, let's find the biggest previous point to add to our point
            else:
                next_point = self
                next_coefficient = 1
                for (previous_coefficient, previous_point) in pervious_points:
                    if previous_coefficient + current_coefficient <= times:
                        if previous_point.x != current_point.x:
                            next_coefficient = previous_coefficient
                            next_point = previous_point
                current_point = current_point.add(next_point)
                current_coefficient = current_coefficient + next_coefficient

        return current_point

Thus we’ve got a super efficient implementation! And now we can perform all the needed operations on an elliptic curve.

Let’s define secp256k1:

secp256k1_curve_config = {
    'a': 0,
    'b': 7,
    'p': 115792089237316195423570985008687907853269984665640564039457584007908834671663
}
x = 55066263022277343669578718895168534326250603453777594175500187360389116729240
y = 32670510020758816978083085130507043184471273380659243275938904335757337482424
n = 115792089237316195423570985008687907852837564279074904382605163141518161494337
g_point = Point(x, y, secp256k1_curve_config)

I’m using only decimal numbers in our examples because they’re intuitive for a human.

So far we’ve implemented everything that we discussed in Part 2. Now let’s implement the actual digital signature algorithm, described in the Part 3.

Sign method of ECDSA:

def sign_message(message, private_key):
    k = random.randint(1, n)
    r_point = g_point.multiply(k)
    r = r_point.x % n
    if r == 0:
        return sign_message(message, private_key)
    k_inverse = find_inverse(k, n)
    s = k_inverse * (message + r * private_key) % n
    return r, s

Verify method of ECDSA:

def verify_signature(signature, message, public_key):
    (r, s) = signature
    s_inverse = find_inverse(s, n)
    u = message * s_inverse % n
    v = r * s_inverse % n
    c_point = g_point.multiply(u).add(public_key.multiply(v))
    return c_point.x == r

Let’s pick some random number as our private key, for example, 123456789012345.

Let our message be 12345.

Do you remember how to get PublicKey from PrivateKey?

private_key = 123456789012345  # any random integer
public_key = g_point.multiply(private_key)
message = 12345  # any integer

Now let’s sign and try to verify:

signature = sign_message(message, private_key)
print('Signature: ', signature)
print('Is valid: ', verify_signature(signature, message, public_key))

It works! You can try to corrupt the signature or the original message and make sure that our algorithm works properly.

Here is the complete code:

import random

def find_inverse(number, modulus):
    return pow(number, -1, modulus)

class Point:
    def __init__(self, x, y, curve_config):
        a = curve_config['a']
        b = curve_config['b']
        p = curve_config['p']

        if (y ** 2) % p != (x ** 3 + a * x + b) % p:
            raise Exception("The point is not on the curve")

        self.x = x
        self.y = y
        self.curve_config = curve_config

    def is_equal_to(self, point):
        return self.x == point.x and self.y == point.y

    def add(self, point):
        p = self.curve_config['p']

        if self.is_equal_to(point):
            slope = (3 * point.x ** 2) * find_inverse(2 * point.y, p) % p
        else:
            slope = (point.y - self.y) * find_inverse(point.x - self.x, p) % p

        x = (slope ** 2 - point.x - self.x) % p
        y = (slope * (self.x - x) - self.y) % p
        return Point(x, y, self.curve_config)

    def multiply(self, times):
        current_point = self
        current_coefficient = 1

        pervious_points = []
        while current_coefficient < times:
            # store current point as a previous point
            pervious_points.append((current_coefficient, current_point))
            # if we can multiply our current point by 2, do it
            if 2 * current_coefficient <= times:
                current_point = current_point.add(current_point)
                current_coefficient = 2 * current_coefficient
            # if we can't multiply our current point by 2, let's find the biggest previous point to add to our point
            else:
                next_point = self
                next_coefficient = 1
                for (previous_coefficient, previous_point) in pervious_points:
                    if previous_coefficient + current_coefficient <= times:
                        if previous_point.x != current_point.x:
                            next_coefficient = previous_coefficient
                            next_point = previous_point
                current_point = current_point.add(next_point)
                current_coefficient = current_coefficient + next_coefficient

        return current_point

secp256k1_curve_config = {
    'a': 0,
    'b': 7,
    'p': 115792089237316195423570985008687907853269984665640564039457584007908834671663
}
x = 55066263022277343669578718895168534326250603453777594175500187360389116729240
y = 32670510020758816978083085130507043184471273380659243275938904335757337482424
n = 115792089237316195423570985008687907852837564279074904382605163141518161494337
g_point = Point(x, y, secp256k1_curve_config)

def sign_message(message, private_key):
    k = random.randint(1, n)
    r_point = g_point.multiply(k)
    r = r_point.x % n
    if r == 0:
        return sign_message(message, private_key)
    k_inverse = find_inverse(k, n)
    s = k_inverse * (message + r * private_key) % n
    return r, s

def verify_signature(signature, message, public_key):
    (r, s) = signature
    s_inverse = find_inverse(s, n)
    u = message * s_inverse % n
    v = r * s_inverse % n
    c_point = g_point.multiply(u).add(public_key.multiply(v))
    return c_point.x == r

# test starts here
private_key = 123456789012345  # any random integer
public_key = g_point.multiply(private_key)
message = 12345  # any integer

signature = sign_message(message, private_key)
print('Signature: ', signature)
print('Is valid: ', verify_signature(signature, message, public_key))

So the implementation of the entire ECDSA algorithm took just 100 lines of code! And it’s perfectly working. This is absolutely the same algorithm as the one used in Bitcoin!

Live Demo

As I promised at the beginning of this article, here is the live demo using only the concepts and formulas described in the article. Just a couple of notes:

Initially, we could only sign integer messages. But in the demo, you can choose to apply a hash function (sha256) to your message. Thanks to it, a message can be a string.
Bitcoin uses slightly different formats of public keys and signatures.
Never use it in a production environment! It is not safe. For production, you must only use well-tested solutions.

The end of the series

I hope this series of articles was very useful for you. At least, I did my best to make it useful. Feel free to share it with friends or use any piece of it anywhere. Just please leave a link to the original article.

Feel free to contact me and ask questions:

exemak@gmail.com
t.me/exemak
Mikhail Karavaev

Blog

Bitcoin Signatures From Scratch (4/4): ECDSA Implementation in Python Using Zero Dependencies

Mikhail Karavaev