Boltzmann machines, spin, Markov chains and all that

The image above displays a set of handwritten digits on the left. They look a bit like being sketched on paper by someone in a hurry and then scanned and digitalized, not very accurate but still mostly readable – but they are artificial, produced by a neuronal network, more precisely a so called restricted Boltzmann machine.

On the right hand side, you see the (core part of) the code that has been used to produce this image. These are about forty lines of code, and there is some code around it which is not shown, but stripping off all the comments and boilerplate code, we could probably fit the algorithm into less than fifty lines of code.

I found this contrast always fascinating. Creating something that resembles handwritten digits sounds incredibly complex, but can be done with an algorithm that can  be quenched into a comparatively short program – how does that work? So I started to dig deeper, trying to understand neuronal networks in general and in particular Boltzmann machines – the mathematical foundations, the algorithm and the implementation.

Boltzmann machines are a rather special class of neuronal networks that have a reputation of being hard to train, but are important from a theoretical point of view, being closely related to seemingly remote fields like thermodynamics, statistical mechanics and stochastical processes. That makes them interesting, but also hard to understand. I embarked on that journey a couple of months ago, and I thought it would be nice to create a series of blog posts on this. My current thinking is to have one or two posts on each of the following topics over the next couple of weeks.

This is a lot of content, and please do not expect to see one post every other day. But maybe I should just start and we will see how it goes….

So let me try to roughly sketch what a Boltzmann machine is. First, it is a neuronal network. As such it is composed of units that can be compared to the neurons in the nervous system. Similar to a neuron, each unit has an input and an output.  The output of a neuron can serve as input for another neuron. In general, every neuron receives inputs from many other neurons and delivers outputs to many other neurons.

Perceptron

The diagram above shows a very simple neuronal network. It consists of four neurons. The three of them on the left serve as input to the network. Think of them as the equivalent of cells in – say – your visual cortex that are activated by some external stimulus. The cell on the right is the output of the network. Its activation is computed based on the outputs of the neurons on the left and some parameters of the network called the weights which model the strength of the connection between the neurons and which we denote by w_i, i = 0, 1, 2, as follows.

o = \sigma(w_0 x_0 + w_1 x_1 + w_2 x_2) = \sigma(w^T X)

Here \sigma is a function called the activation function and X is the vector that is formed by x_0, x_1 and x_2. There are a few standard choices for the activation function, a common one being the sigmoid function.

How is such a network applied? Let us look at this for a problem which is the “Hello World” of machine learning – recognition of handwritten digits. In that problem, you start with a collection of digitized images of handwritten digits, like those that are known as the MNIST database. These images have 28×28 = 784 pixels. We want to design a neuronal network that can classify these images. Such a network should consequently have 784 input units and 10 output units (bear with me that I did not produce a picture of that network, even though it would probably be fun to do this with Neo4J). We present an image to the network by setting the value of the input unit i to the intensity of pixel i. Our aim is to adjust the weights in such a way that if the image represents digit n, only output n is significantly different from zero. This will allow us to classify unknown images – we simply present the image to the network and then see which output is activated.

But how do we find the correct values for the weights? We need weights that connect each of the ten output units to each of the 784 input units, i.e. we have 7840 weights. Thus we are looking for a point in a 7840 dimensional vector space – not easy. Here the process of training comes into play. We take our set of sample images and present them to the network, with initially randomly chosen weights. The output will then not be what we want, but differ from the target output by an error. That error can be expressed as a function of the weights.  The task is then to find a minimum of the error function, and there are ways to do this, most notably the procedure which is known as gradient descent.

As the name suggests, we need the gradient of the error function for that purpose. Fortunately, for the type of neuronal network that we have sketched so far, the gradient can be computed fairly easily – in fact the activation function is chosen on purpose to make dealing with the derivatives easy. We end up with a comparatively simple training algorithm for this type of network and maybe I will show how a simple implementation in Python in a later post – but for now let us move on to Boltzmann machines.

Boltzmann machines or more precisely restricted Boltzmann machines (RBMs) are also composed of units and weights, but work a bit differently.  There are inputs, which are usually called visible units. But there are no classical outputs. Instead, there is a layer of units that is connected to the visible units and is called hidden units, as shown in the following picture.

RestrictedBoltzmannMachine

In the example of handwritten digits again, you would again have 784 visible units. However, the hidden units would not obviously correlate with the digit represented by the input. Instead, you would have a more or less arbitrary number of hidden units, say 300.

During training, you present the network one image, again by setting the values of the visible units (the input) to the pixel intensities. You then compute the value of the hidden units – but you do not do this deterministically, but bring in some random element. Roughly speaking, if the combined input to a hidden unit is p, you set the value of the unit to one with probability p. Then this process is repeated, this time starting from the hidden units. This gives you certain values for the visible units. You then compare this value to the original value and try to adjust the weights such that you get as closely as possible to the original input (this is not exactly true and a not very precise description of one of the possible learning algorithms called contrastive divergence, but we will get more precise later on).

If you succeed, you will be able to reconstruct the value of the visible units from the values of the hidden units. But there are less hidden units, so that the network has apparently learned a condensed representation of the input that still captures the essence of the input. In the case of digits, you can visualize the state of the network after training and obtain something like

ToyRBM_Weights

These are visualizations of the weights connected to some of the hidden units of a Bolzmann machine after the training phase. We can see that some of the units have appearantly “learned” some characteristic features of the digits, like the vertical stroke that appears in the digits one and seven. These features can be used for several purposes.

First, we can use the features as input for other neuronal networks. We now have only 300 inputs instead of originally 784, and that might simplify the problem a bit and make the process of training the next network easier.

Or, we could start with some random values for the hidden units and calculate the resulting values of the visible units to create sample images – and in fact this is basically what I have done to produce the samples at the top of this post (using an algorithm called PCD, but we will get to this).

Note that Boltzmann machines differ significantly from the type of neural network that we have considered earlier. One major difference is that to train a Boltzmann machine, you do not need to know the digit that the image represents. At no point have we used the information that the images in our database represent ten different digits – that information was not used in the design of the network nor during the training. This approach to machine learning is called unsupervised learning and is obviously very versatile – you do not have to tell the machine what the structure of the data is, it will detect the features independently.

Second, a Boltzmann machine can not only classify images, but can also create images that resemble a given set of input data. That can be used to reconstruct partially available images or to create new images from scratch – networks with this ability are called generative networks.

The price we have to pay for this is that Boltzmann machines are hard to train. The point is that we can still define some sort of error function, but we cannot easily calculate its gradient – there is no straightforward analytic expression for it that could be evaluated within a reasonable amount of time (if you know statistical physics a bit, that might remind you of the partition function, and that is more than pure coincidence, as I will show you in a later post). So we need to approximate the gradient. Technically, the gradient is an integral

\int f(x) dx

for some known, but complicated function f. Even if we cannot find an analytic expression for this, we can try to approximate it. Your first idea might be “Riemann sums”, but it turns out that this is not a good idea, as our function lives in the space of all weights which has a very high dimension. Instead, we will use an approach called Monte Carlo integration where we represent the integral as an expectation value, draw a sample and approximate the expectation value by the sample average. This is where stochastical methods like Markov chains will come into play.  And finally we will see that the behaviour of our network during training has some striking analogies with the behaviour of certain physical systems like solids exposed to a magnetic field at low temperatures, which are described by a model called the Ising model, and learn how techniques that physicists have  developed for this type of problems apply to neuronal networks.

That is it for today – I hope I could give you a rough idea of what is ahead of us. At least I hope that you start to be curious how all this works out – so looking forward to the next post where I will start with some background from physics.

Keys in the bitcoin network: the public key

In my last post, we have looked in some detail at the private key – how it is generated and how it can be decoded and stored. Let us now do the same with the public key.

Recall that a public key is simply a point on the elliptic curve SECP256K1 that is used by the underlying ECDSA algorithm – in fact it is obtained by multiplying the generator point on the curve by our private key. As any point on the curve, it therefore has an x-coordinate and a y-coordinate, both being 32 bytes unsigned integers.  So one way to encode the public key would be as follows.

  • take the x-coordinate as a point, represented by an integer smaller than p
  • convert this into a 32 byte hexadecimal string, using for instance big endian encoding
  • do the same for the y-coordinate
  • and concatenate these two strings to obtain a single 64 byte hexadecimal string

This encoding is simple, but it has a drawback. Remember that we encode not just a random pair of integers, but a point on the curve, so the x-coordinate and y-coordinate are related by the curve equation

y^2 = x^3 + ax + b

Thus given x, we almost know y – we know the square of y modulo p, and there can be at most two different roots of this equation. So we could reconstruct y if we have x and an additional bit that tells us which of the two solutions we need.

Let us now assume that p is odd. If y is a solution of the equation for a given value of x, then p – y (which is -y modulo p) is the second solution. As p is odd, exactly one of the two numbers y and p – y is even. We can therefore use an additional bit that is equal to y modulo 2 to distinguish the two solutions. It is convention to store this bit in a full additional byte, using the value 2 if y is even and the value 3 if y is odd, so that we obtain a representation of the public key (and in fact any other point on the curve) in at most 33 bytes: at most 32 bytes for the value of the x-coordinate and the additional byte containing the value of y modulo 2. This representation is called the compressed representation (see for instance the publication of the SECG, section 2.3).

If there is a compressed representation, you might expect that there is also an uncompressed representation. This is simply the representation that we have described above, i.e. storing both x and y, with an additional twist: to be able to distinguish this from a compressed representation that always starts with 0x02 or 0x03, a leading byte with value 0x04 is added so that the total length of an uncompressed representation is at most 65 bytes. Since version 0.6.0, the bitcoin reference implementation defaults to using compressed keys (see the function CWallet::GenerateNewKey).

Let us summarize what we have learned so far in a short Python code snippet that will take a private key (stored as integer in the variable d), calculate the corresponding point on the elliptic curve SECP256K1 using the ECDSA library and create a compressed representation of the result.

#
# Determine the public key from the
# secret d
#
import ecdsa
curve = ecdsa.curves.SECP256k1
Q = d * curve.generator
#
# and assemble the compressed representation
#
x = Q.x()
y = Q.y()
pubKey = x.to_bytes(length=32, byteorder="big")
pubKey = binascii.hexlify(pubKey).decode('ascii')
if 1 == (y % 2):
    pubKey = "03" + pubKey
else:
    pubKey = "02" + pubKey
print("Compressed key:  ", pubKey)

This way of encoding a public key is in fact not specific to the bitcoin network, but a standard that is used whenever a point on an elliptic curve needs to be encoded – see for instance RFC5480 by the IETF which is part of the X.509 standard for certificates.

However, this is still a bit confusing. If you known the nuts and bolts of the bitcoin protocol a bit, you will have seen that participants publish something that is called an address which is a string similar to

mx5zVKcjohqsu4G8KJ83esVxN52XiMvGTY

That does not look at all like a compressed or uncompressed public key. We are missing something.

The answer is that an address is in fact not a public key, but it is derived from a public key. More precisely, it is an encoded version of a hash value of the public key. So given the address, it is easy to verify that this address belongs to a certain public key, but it is very hard to reconstruct the public key given the address.

To understand the relation between a public key and an address better, it is again time to take a look at the source code of the reference client. A good starting point is the RPC method getnewaddress. This method is defined in the file wallet/rpcwallet.cpp and creates an instance of the class CBitcoinAddress which in turn is – surprise – derived from our old friend CBase58Data. The comments are quite helpful, and it is not difficult to figure out that a bitcoin address is obtained as follows from a public key.

  • create a hexadecimal compressed representation of the public key
  • apply a double hash to turn this into a sequence of 20 bytes – first apply the hash algorithm SHA256, then RIPEMD160 (this is called a Hash160 in the bitcoin terminology as the output will have 160 bits)
  • add a prefix to mark this as a public key address – the prefix is again defined in chainparams.cpp and and is zero for the main network and 111 for the test networks
  • take the hash256 checksum and append the first four bytes
  • apply Base58 decoding to the result

This is already very similar to what we have seen before and can be done in a few lines of Python code.

def hash160(s):
    _sha256 = hashlib.sha256(s).digest()
    return hashlib.new("ripemd160", _sha256).digest()
#
# Apply hash160
#
keyId = hash160(bytes.fromhex(pubKey))
#
# Append prefix for regtest network
#
address = bytes([111]) + keyId
#
# Add checksum
#
chk = hash256(address)[:4]
#
# and encode
#
address = btc.utils.base58Encode(address + chk)
print("Address:         ", address)

Heureka! If we run this, we get exactly the address mx5zVKcjohqsu4G8KJ83esVxN52XiMvGTY that the bitcoin client returned when we started our little journey at the beginning of the post on private keys.

As always, the full source code is also available on GitHub repository. If you want to run the code, simply enter

$ git clone https://github.com/christianb93/bitcoin.git
$ cd bitcoin
$ python Keys.py

That was it for today. We have now covered the basics of what constitutes participants in the bitcoin network. In the next few posts in this series. we will look at the second main object in the bitcoin world – transactions. We will learn how to interpret transactions, and will eventually be able to manually create a transaction to instruct a payment, sign it, hand it over to our test network and see how it is processed.

Keys in the bitcoin network: the private key

In my last post, I have shown you how arithmetic on elliptic curves can be used to create and verify digital signatures. We have seen that every party that creates a signature is represented by a private key – kept securely – and a public key, which is made available to everyone who wants to verify the signature. In a blockchain, digital signatures are used to verify ownership of bitcoins, and therefore private and public keys play a pivotal role in the bitcoin network. Bitcoin transaction outputs refer to public keys, and only the person that is in control of the matching private key can spend the bitcoin. Thus it is worth to take a closer look at how keys are represented in the bitcoin protocol.

Bitcoin uses the ECDSA signature algorithm to sign messages and verify signatures. Therefore a private key in the bitcoin network is simply an integer. More precisely, it is an integer between one and the order of the generator of an elliptic curve SECP256K1 (if that sounds like gibberish to you, you should read my previous post on this). Generating private keys is therefore very easy – you simply randomly select a 32 byte integer until you find one which is below that order (I am cheating a bit – obviously you need a good source of random numbers for this which makes it hard to predict the private key). I have done that for you:

d = 103028256105408389446438916672504271192164767440296751065327418112299269382535

Of course that is not a private key that I really use – it would not be smart to publish it if it were. But it is a perfectly valid private key.  Unfortunately, this is not the way how a private key is typically stored and presented by a bitcoin client. In fact, what I did to create this key is to run the commands


$ bitcoin-cli -regtest getnewaddress "myAccount"

mx5zVKcjohqsu4G8KJ83esVxN52XiMvGTY

$ bitcoin-cli -regtest dumpprivkey "mx5zVKcjohqsu4G8KJ83esVxN52XiMvGTY"

cVDUgUEahS1swavidSk1zdSHQpCy1Ac9XSQHkaxmZKcTTfEA5vTY

after installing the reference bitcoin core software on my computer. The last line is the private key – and that looks very much different from the number d above (you do not necessarily need a bitcoin installation to follow this post, but you will definitely need it to run all the examples in future posts, so this might be a good point in time to stop and install it. So go to the download page, get the version for your machine and install it. Do not forget to start the bitcoin daemon in regtest mode, if you want to avoid downloading the full blockchain with more than 200 GBytes at the time of writing. In my bitcoin configuration file bitcoin.conf located in the directory ~/.bitcoin, I have set the options

regtest=1
server=1

to tell bitcoin to accept RPC commands and to run in regression test mode. But back to keys now …).

The funny string that the bitcoin client will present you as the private key is in fact an encoded private key using the WIF format (wallet interchange format).  Let us try to understand how we can convert this into the number d displayed above.

The first thing you need to know about a WIF encoded private key is that it is encoded using the Base58 standard. Similar to Base64 (which an official IETF standard described in RFC4648), this is a standard to encode a number in a way that can easily be transmitted over channels like e-mail or even printed on paper without having to deal with binary values. Essentially, the idea is that we use an alphabet of 58 ASCII characters and to convert the number to the base 58, representing each digit by the corresponding character from this alphabet. In addition, there is some logic to handle leading zeros, more precisely to avoid that they are dropped during the conversion. If you want to see all this in detail, the authorative answer is (as always) the source code of the module base58c.cpp in the C++ reference implementation which is hosted on GitHub.

To do this in Python, we have – as always – several choices.  We can search for a library that performs Base58 encoding and decoding  – for instance https://github.com/keis/base58. For the sake of demonstration, I have created my own routines. To decode, i.e. to turn a Base58 string into a sequence of bytes, the following code will do.


BASE58_ALPHABET = '123456789ABCDEFGHJKLMNPQRSTUVWXYZabcdefghijkmnopqrstuvwxyz'
def base58_decode(s):
    #
    # Strip off leading 1's as these represent leading
    # zeros in the original
    #
    zeros = 0
    while (zeros < len(s)) and (s[zeros] == '1'):
        zeros = zeros + 1
    s = s[zeros:]
    #
    # We first turn the string into an integer
    #
    value, power = 0, 1
    for _ in reversed(s):
        value += power * BASE58_ALPHABET.index(_)
        power = power * 58
    #
    # Now convert this integer into a sequence of bytes
    #
    result = value.to_bytes((value.bit_length() + 7) // 8, byteorder='big')
    #
    # and append the leading zeros again
    #
    for _ in range(zeros):
        result = (0).to_bytes(1, 'big') + result
    return result

I have stored this in a module btc.utils for later use, which you can find, along with the other examples from this post, in my GitHub repository.

Let us now apply this to our example WIF file.

import binascii
import btc.utils

#
# The WIF encoded private key
#
wif = "cVDUgUEahS1swavidSk1zdSHQpCy1Ac9XSQHkaxmZKcTTfEA5vTY"
print("WIF:    ", wif)
#
# Convert into a sequence of bytes
#
b = btc.utils.base58_decode(wif)
#
# and into hex
#
h = binascii.hexlify(b).decode('ascii')
print("Hex:    ", h)

If you run this code, you will find that the output has 76 characters, i.e. 76 / 2 = 38 bytes. That cannot be quite right, because we expect that our private key is an integer with 32 bytes only. So there are six extra bytes. Where do they come from?

To get used to it, let us again try to find the answer in the source code of the reference implementation (you can browse the code in the GitHub repository online or (recommended) clone to obtain a local copy so that you can use tools like grep and their friends. As a starting point, remember that we have received our private key via the bitcoin-cli tool that communicates with the bitcoin server via RPC calls, and that we have used the RPC method dumpprivatekey. So let us search for that.


(cd bitcoin/src ; grep -R -I "dumpprivkey" *)

That will give you a few matches in two files,wallet/rcpwallet.cpp and wallet/rpcdump.cpp. If you open these two files and look at the code, you will find that the former refers to the latter. This function first retrieves the key from the wallet and then creates an instance of the class CBitcoinSecret (which is derived from CBase58Data and invokes its ToString() method to obtain the textual representation of the key that is then returned.

One more usage of grep will tell you that the code for these classes is located in the file base58.cpp which we have already met. The constructor calls CBitcoinSecret::SetKey() and the method ToString() is implemented in the base class, so we also need to look at CBase58Data::ToString(). Going carefully through this code, we find that the data is actually composed of four parts.

WIFFormat

The first byte is a version number which is used to distinguish a private key used for the testnet (239) from a key for the productive network (128) (look up the values in chainparams.cpp). The next few bytes are the actual secret d, encoded as a hexadecimal string using big endian encoding (i.e. the most significant octet first). The next byte is a flag that describes whether the public key that belongs to this private key should be stored in compressed format or not. We will get back to this point in a later post on public keys and addresses, for the time being we can safely ignore this byte.

The last four bytes are again interesting. They form a checksum for the remainder. These four bytes are obtained by applying what is called a Hash256 in the bitcoin language and which is just a double SHA256 hash, and then taking the first four bytes of the result. Thus in order to turn the WIF string into a number, we have to decode it using Base58, strip off the last four bytes, verify the checksum, strip off one additional byte, remove the trailing version number and convert the remaining hexadecimal string into an integer.

import binascii
import btc.utils
import hashlib

def hash256(s):
    return hashlib.sha256(hashlib.sha256(s).digest()).digest()

#
# The WIF encoded private key
#
wif = "cVDUgUEahS1swavidSk1zdSHQpCy1Ac9XSQHkaxmZKcTTfEA5vTY"
print("WIF:       ", wif)
#
# Convert into a sequence of bytes
#
b = btc.utils.base58_decode(wif)
#
# and into hex
#
h = binascii.hexlify(b).decode('ascii')
print("Hex:       ", h)
#
# Strip off checksum
#
chk = h[-8:]
print("Checksum:  ", chk)
h = h[:-8]
#
# and verify it
#
_chk = hash256(bytes.fromhex(h))[:4]
assert(_chk == bytes.fromhex(chk))
#
# Strip off version byte
#
print("Version:   ", int(h[:2], 16))
h = h[2:]
#
# and compression flag
#
h = h[:-2]

d = int.from_bytes(bytes.fromhex(h), "big")
print("Secret:    ", d)

If you run this, you should get the value for d that we started with at the beginning of this post.

We have now seen how private keys can be generated and encoded. Typically private keys are kept in a wallet, but they could also be printed out in their WIF encoded form and stored offline – you could even create a private key on a machine not connected to any network and store it in this way. This is called a paper wallet in the bitcoin world.

But what about the public key? If you want to receive bitcoins, the payee needs access to your public key or at leat to a condensed version of it called the bitcoin address.  We have seen that for the ECDSA algorithm, the public key can be calculated given the private key, and we will see that the address can in turn be calculated given the public key. In my next post, I will guide you through this process.

A primer on elliptic curve cryptography: practice

In the last post, we have looked a bit at the theory behind elliptic curves. In this post, we will now see how all this works down to earth and use Python to actually run some calculations.

The first thing that we need is an explicit formula for the addition of two points on an elliptic curve. We will not derive this here, but simply give you the result – see for instance [1] for more details. Given two points (x_1, y_1) and (x_2, y_2), the coordinates of their sum (x_3, y_3) can be determined as follows.


inv = inv_mod_p(x2 - x1, p)
x3 = ((y2 - y1)*inv)**2 - x1 - x2
y3 = (y2 - y1)*inv*(x1 - x3) - y1

Here we assume that we have a function inv_mod that will give us the inverse modulo some prime number p which is the number of elements of our base field. Typically, this can be done using the so-called extended euclidian algorithm.

However, there are a few special cases we need to consider. This is apparent if we look at this formula in more detail – what happens if the inverse does not exist because the points x_1 and x_2 are equal?

This happens if we want to add two points that have the same x-coordinate. There are two cases we need to consider. First, the y-coordinates could be equal as well. Then we are trying to add a point to itself, and we need to apply the formula provided in  [1]  for that special case. Or the y-coordinate of the second point is minus the y-coordinate of the first point. Then we try to add a point to its own inverse. The result will be the neutral element of the group which is usually called the point at infinity (there is a reason for this: if we embedd the curve into a projective plane, this point will in fact be the intersection of its completion with the line at infinity in the projective plane…).

To describe points on an elliptic curve, we need both coordinates, the x and the y-coordinate. Thus it makes sense to implement a class for this purpose. Instances of this class need to store the x- and y-coordinates and – for convenience – a boolean that tells us whether a point is the point at infinity. So our full code for this class looks as follows.


class CurvePoint:

    def __init__(self, x, y, infinity = False):
        self.x = x
        self.y = y
        self.infinity = infinity

    def __add__(self, other):
        #
        # Capture trivial cases - one of the points is infinity
        #
        if self.infinity:
            return other
        if other.infinity:
            return self
        #
        # First check whether we are adding or doubling
        #
        x1 = self.x
        x2 = other.x
        y1 = self.y
        y2 = other.y
        infinity = False
        if (x1 - x2) % p == 0:
            #
            # Are we talking about doubling or addition
            # of the inverse?
            #
            if (y1 + y2) % p == 0:
                infinity = True
                x3 = 0
                y3 = 0
            else:
                inv = inv_mod_p(2*y1, p)
                x3 = (inv*(3*x1**2 + a))**2 - 2*x1
                y3 = (inv*(3*x1**2 + a))*(x1 - x3) - y1
        else:
            #
            # Standard case
            #
            inv = inv_mod_p(x2 - x1, p)
            x3 = ((y2 - y1)*inv)**2 - x1 - x2
            y3 = (y2 - y1)*inv*(x1 - x3) - y1

        return CurvePoint(x3 % p, y3 % p, infinity)

As already mentioned, that assumes that you have a function inv_mod_p in your namespace to compute the inverse modulo p. It also assumes that the variables p, a and b that describe the curve parameters are somewhere in your global namespace (of course you could introduce a class to represent a curve that stores all this, but let us keep it quick and dirty at this point).

Now having these routines, we can actually do a few example and verify that the outcome is as expected. We use a few examples with low values of p from [1] .

#
# Define curve parameters
#

p = 29
a = 4
b = 20
#
# and add a few points
#
A = CurvePoint(5,22)
B = CurvePoint(16, 27)
O = CurvePoint(0,0,infinity=True)
C = A + B
assert(C.x == 13)
assert(C.y == 6)
assert(C.infinity == False)
C = A + A
assert(C.x == 14)
assert(C.y == 6)
assert(C.infinity == False)
A = CurvePoint(17,19)
B = CurvePoint(17,10)
C = A + B
assert(C.infinity == True)
A = B + O
assert(A.x == B.x)
assert(A.y == B.y)
assert(A.infinity == B.infinity)
A = O + B
assert(A.x == B.x)
assert(A.y == B.y)
assert(A.infinity == B.infinity)

For the sake of demonstration, we have shown how to build elliptic curve arithmetic from scratch. We could now proceed to implement multiplication and the ECDS algorithm ourselves, but as so often, there is a Python library that will do this for us. Well, there is probably more than one, but I like the Python ECDSA library maintained by Brian Warner.

The most basic classes in this library are – you might have guessed that- curves and points. Curves are initialized providing the basic parameters p, a and b. Then points are created by specifying a curve and the x and y coordinates of the points. Thanks to operator overloading, points can then be added and multiplied with integers using standard syntax. Here is a code snippet that reproduces the first of our examples from above using the ECDSA library.


import ecdsa

#
# Create a curve with parameters p,a and b
# with the ECDSA library
#
curve = ecdsa.ellipticcurve.CurveFp(p,a,b)
#
# Define two points and add them
#
A = ecdsa.ellipticcurve.Point(curve, 5, 22)
B = ecdsa.ellipticcurve.Point(curve, 16, 27)
C = A + B
assert(C.x() == 13)
assert(C.y() == 6)
assert(C != ecdsa.ellipticcurve.INFINITY)

Very easy – and very useful. And of course reassuring to see that we get the same result than our hand-crafted code for the arithmetic above gave us.

But this is not yet all, the library can of course do much more. Let us use it to create a signature. For that purpose, we obviously need a reasonable large value for the prime p, otherwise we could easily use a brute-force attack to determine our private key from the public key. In my previous post on the theoretical foundations, I have already mentioned the papers published by the SECG, the Standards for efficient cryptography group. This group has published some standard curves that we can use. One of them is the curve SECP256K1 which is a curve over a prime field F_p with

p = 2^{256} - 2^{32} - 2^{9} - 2^{8} - 2^{7}-2^{6} - 2^{4} - 1 = 115792089237316195423570985008687907853269984665640564039457584007908834671663

This curve is the curve which is used by the bitcoin protocol. As many other standardized curves, it is hard-coded in the ECDSA library.  To get it, use the following code.


#
# Get the standard curve SECP256K1
# and its parameters
#
curve = ecdsa.curves.SECP256k1
G = curve.generator
p = curve.curve.p()
a = curve.curve.a()
b = curve.curve.b()
n = G.order()

Let us now apply what we have learned about signatures. First, we need a private key and a public key determined from it. Thus we pick a random number that will be our secret and multiply the generator of the curve by it to get our public key.

#
# Determine a private key and a public key
#
d = ecdsa.util.randrange(n-1)
Q = d*G
pKey = ecdsa.ecdsa.Public_key(G, Q)
sKey = ecdsa.ecdsa.Private_key(pKey, d)

Next we need a hash value that we will sign. Usually, we would derive this value from a message using some cryptographic hash function like SHA256, but we will simply simulate this by drawing a random number h. We can then use the method sign of the ECDSA private key object to create a signature. This will return a signature object from which we can retrieve the values r and s.

h = ecdsa.util.randrange(n-1)
k = ecdsa.util.randrange(n-1)
signature = sKey.sign(h, k)
r = signature.r
s = signature.s

Let us verify that the algorithm really works as described  the previous post. So we first need to multiply our randomly chosen integer k with the generator of the curve, the number r should then be the x-coordinate of this point. We then invert k modulo n and multiply the result by h + dr. This should give us s.

_r = (k*G).x() % n
assert(_r == r)

w = inv_mod_p(k, n)
_s = ((h+d*r)*w) % n
assert(_s == s)

If you run this code, you will hopefully see that the assertions pass. Our code works, and produces the same result as the ECDSA library. That is already reassuring. Having gone so far, we can now of course also verify the signature – again we do this once using the ECDSA library and once using our own code.

#
# Now we manually verify the signature
#
w = inv_mod_p(s, n)
assert(1 == (w*s % n))
u1 = w * h % n
u2 = w * r % n
X = u1*G + u2*Q
assert(X.x() == r)

#
# Finally we verify the signature using the
# lib
#
assert(pKey.verifies(h, signature) == True)

That is it for today. We have seen how elliptic curves can be used in practice to create and verify digital signatures and have looked at the ECDSA library that offers ready made functions for that purpose. The full source code can by downloaded from GitHub.

In my next post, I will look at the way how private and public ECDSA keys appear in the bitcoin protocol. If you want to learn more and play around a bit with elliptic curves in the meantime, I recommend the online tool [2]. You might also want to take a look at Andrea Corbellinis excellent post on elliptic curve cryptography.

1. D. Hankerson, A. Menezes, S. Vanstone,
Guide to elliptic curve cryptography, Springer, New York 2004
2. Elliptic curve point addition online tool at https://cdn.rawgit.com/andreacorbellini/ecc/920b29a/interactive/modk-add.html