Arithmetic on Quantum Computers: Addition, Faster

22 min readOct 11, 2018

In this story, I will be exploring a method to increase the efficiency of addition on a quantum computer using the quantum Fourier Transform, and implementing it on a quantum computer using QISkit, IBM’s quantum computing API for Python.

In case you missed part one of this series, an article on a simple adder circuit that also serves as a gentle introduction to quantum computing, check it out here.

Let’s get started.

Further down the rabbit hole

In my last article we walked through a simple adder circuit — while it did the job, it was not any different from addition implemented on a classical computer. Quantum computers are not much use if they only do what can be done on a classical computer, so, this time, we will use the ‘quantumness’ of qubits to make an addition program that runs faster and uses less qubits for the same inputs.

Before that, however, we will build an understanding of the Fourier transform. This is not necessary reading; you can skip to the implementation of the addition program without any loss, although a grasp of the underlying principles while writing the addition program will prove helpful.

As we move further into the world of quantum computing, it is only natural that algorithm implementation become more ‘quantum’ through the use of quantum mechanical effects like superposition and entanglement. Along with those large and complicated sounding words comes a lot of math as well, so be prepared: quantum computing is not for the faint-hearted.

Don’t bother asking for recipes anymore

If you have some fruit, milk, and a recipe, you can make a milkshake. But what if you want to know what was put into a milkshake made by someone else? In other words, can you find a recipe for a milkshake just by drinking it? That’s a little harder than making a milkshake.

If you had a machine that could separate the milkshake into its constituents using ‘filters’, you could construct a recipe from a sample of the milkshake. The machine would not work all the time, though. If you try feeding it with a brownie, it won’t work, because an intermediary process (baking) changes the ingredients, making it impossible to retrieve them. Only foods which do not undergo baking, frying, or any other non-reversible transformation can be filtered by the machine.

Additionally, the filters must be set to capture all ingredients in the food the machine is fed with. If you try to analyse an Oreo milkshake, but forget to set a filter to capture Oreos, the machine will not give you a complete recipe.

The Fourier transform works similarly to our hypothetical machine. The are are only a few differences:

Instead of taking food as an input, it takes a wave.
These waves could be any kind of wave, like sound, radio, or anything else.
In the place of a recipe, the Fourier transform separates the complicated input wave into many simpler waves.

What is the difference between complicated and simple? Look at this:

This wave looks periodic, but it is hard to describe. Once we pass it through our Fourier transform machine, we get the three waves that add up to make the one above:

The three output waves look a lot simpler.

But how does the machine actually do this?

Wrapping waves around circles

Imagine a simple cosine wave, shifted by one unit upwards, with a frequency of 3 cycles per second:

Now take the wave, and wrap it around a circle so that the height of the wave at any point of time is the distance from the origin for the new graph. This wrapping can take place at any speed we want it to; let’s start with 1 rotation per second, meaning our new graph covers 0.25 cycles of the sine graph:

If you know what polar coordinates are, the new graph can be represented by the equation below:

Suppose that the graph is now actually some sort of wire, and that we want to find how the center of mass of the wire changes as we change the number of rotations we make per second. To give you an idea of what it looks like, here are a few pictures of the wave wrapped at different rotation speeds. The purple dot represents the center of mass:

1 rotation per second

2 rotations per second

2.25 rotations per second

2.5 rotations per second

2.75 rotations per second

3 rotations per second

As you can see, the graph lines up nicely when the rotation frequency is 3 times per second, and the maximum x-value of the center of mass is reached at this point.

Remember the equation for our original wave?

Our wrapped graph also has a rotation frequency of 3 cycles per second; is something special happening?

It turns out that when the rotation frequency is equal to the frequency of the function, the graph wrapped around a circle lines up to form a single shape (in this case a cardioid) and the x-value of the point giving the center of mass is at its maximum. This is the working principle of the Fourier transform; by wrapping a wave around a circle and looking at the center of mass of the wrapped wave, we can find out the periodicity of its different components.

Pretty cool, right?

Now, we can involve a little more math; this will give us the reason why the Fourier transform can be implemented so easily using qubits, and an understanding of how the Fourier transform is implemented on a quantum computer.

Rotations, wrapping and Euler’s formula

In the previous section, we learnt that wrapping a wave around a circle with different rotation speeds helps us compute the Fourier transform of that wave. But how do we describe this action, mathematically?

We need some expression that is equivalent to a rotation about the origin of the complex plane (with real numbers on the x-axis, and imaginary ones on the y-axis) that changes as we change some parameter, so that we can decide by how much and how fast we want to rotate about the circle. Euler’s formula is a beautiful and very convenient way of expressing rotations:

This is equivalent to a rotation by angle θ on the complex unit circle. If we replace θ with t, then every 2π seconds one complete rotation would be made. To simplify our job, we can move the 2π to the exponent, so that every second we move around the circle once:

If we introduce another variable s, we can control how fast we move around the circle:

This allows us to move around the circle s times every second. But this only allows us to rotate anticlockwise, so to rotate clockwise we add a minus sign to the exponent.

Now to wrap a wave around a circle at a frequency s, all we have to do is multiply the function f(t) representing the wave with our rotation function:

Now we have to figure out a way to express the ‘center of mass’ of the function. The easiest way would be to sample from the wrapped up wave n times, and take the average of the values you got:

If we try to make our center of mass calculation more accurate by sampling more points, we can take the entire wrapped wave into account by taking the integral of its function, and dividing by the time interval that we are considering:

So now we have an expression for the center of mass for a wave wrapped around a circle. If we choose not to divide by the time interval, we get the Fourier transform of the wave for a frequency s:

Awesome!

What if we have the Fourier transform of a wave, and want to know what the original wave looked like? That would involve reversing what we did to get the Fourier transform, and so is called the inverse Fourier transform:

So from a time-intensity function f(t) used to represent a wave, we get a frequency-intensity function ĝ(f) called the Fourier transform of f(t), and we can reverse the process using the inverse Fourier transform.

How does this help us add two numbers?

Think of numbers as constant functions, like the ones below:

No matter what the input value is, the function would always give the same value. If we plot these two functions on a graph, we get two straight lines:

The lines can also be represented as part of infinitely large circles:

If we compute the Fourier transforms of these two ‘waves’, we get something pretty cool:

It’s pretty easy to see that for any nonzero s, the value of the Fourier transforms of both functions are zero. If we add the two Fourier transforms together:

One of the properties of the Fourier transform is that the sum of the Fourier transforms of two signals is equal to the Fourier transform of the sum of the signals:

This property is called linearity, and is what allows us to add using the Fourier transform.

Fourier transforms for computers

While thinking of a number as a constant function is fine for the purposes of our understanding, storing this function on a quantum computer would be a waste of qubits as we would store multiple instances of the same value corresponding to different points in time. Instead, we use a slightly different version of the Fourier transform that works with discrete data sets, called the discrete Fourier transform, and instead of using an integral we use the summation operator Σ.

Mathematically, this version of the Fourier transform is written as:

This looks hard, and is probably the seemingly hopeless kind of equation that makes people take up paintball instead of quantum computing, but it isn’t too hard to calculate:

y__k and x__j are complex numbers.
j and k are indices, ranging from 0 to N-1.
N is the length of the data set x.

A simple example will help our understanding. Given x = {1, 2}, can we calculate its discrete Fourier transform?

Here N is 2, as there are two elements in the data set x.
j and k range from 0 to 1.
x__0 is 1, and x__1 is 2.

Our goal is to find y__0 and y__1 :

Now we solve for y__1:

So the discrete Fourier transform of {1, 2} is {3/√2, -1/√2}. Not that hard, right?

Qubit superposition and the quantum Fourier transform

The quantum Fourier transform is almost the same as the discrete Fourier transform, with only one difference:

Instead of containing numbers, the data set is a vector, with each element in the vector representing the probability of the state being measured as that combination of 0 and 1, and is often called ψ instead of x .

It can be computed almost the same way as well for any qubit state ψ:

What do I mean by a vector representing the state of a qubit?

When you measure a qubit, it chooses either the |0> and the |1> state and shows you that value. The reason for this lies more in quantum mechanics than computing, so for now we’ll take it as a given. But can we predict which one the qubit will choose, based on its state?

Intuitively, the larger the value of α, the more likely it would be that the qubit picks the |0> state, and the larger the value of β, the more likely it would be to pick the |1> state. Close — the values of α*α and β*β represent the probability that the qubit will choose the |0> and the |1> state respectively. So now we can express the relationship between α and β:

Why the complex conjugate *? α and β need not be positive, or even real. Think about what happens when we simply use α and β, or their squares, as probabilities. Would these values fit the definition of a probability — a number between 0 and 1?

To make the concept of the quantum Fourier transform clear, we can use an example. Consider the 1-qubit state |ψ> = α|0> + β|1>. Can we compute its quantum Fourier transform?

We can create a matrix that represents the 1 qubit quantum Fourier transform operation:

If we add a 1/√2 factor to the matrix as well:

Solving for the matrix yields the below:

This matrix is extremely important in quantum computing, and is called the Hadamard matrix after its discoverer. It is implemented in QISkit like this, for the first qubit of quantum register a in a quantum circuit qc:

qc.h(a[0])

If we consider a 2 qubit state ψ:

Calculating the quantum Fourier transform of this state would give us F(ψ):

We can use a matrix to represent the 2 qubit quantum Fourier transform as well:

I will leave the derivation of this matrix to you.

Converting the quantum Fourier transform into gates

If you look into the matrices we derived in the previous section, you might realize that we never change the value of α*α + β*β ; it always remains 1. If you compute the determinants of the matrices, you will arrive at the same result: the absolute values of the determinants are equal to 1.

Since the magnitude of the state vector ψ remains unchanged, we can term the transformations as rotations on the complex plane; this is equivalent to performing a series of operations of the form e^(iθ) , where θ takes different values. If we go back to the matrices again, but this time represent each matrix coefficient as a rotation, we see some interesting patterns emerge:

We can see that for a general qubit state represented by a vector with N possible outcomes, the quantum Fourier transform can be given by this matrix:

We can verify this by calculating the matrix for a 3 qubit state:

The N-dimensional matrix is equivalent to the function representing the quantum Fourier transform:

If you remain unconvinced, try using both the matrix and the function to calculate the quantum Fourier transform of a simple 1 qubit state, say |ψ> = (1/1.41421)|0> + (1/1.41421)β|1>. The results from both methods will be the same.

We will use this definition of the quantum Fourier transform of a general N-dimensional vector representing the quantum state ψ to implement a quantum adder, by decomposing this matrix into a series of phase rotation gates, which can be implemented using the quantum R gate:

By combining many R gates with different values of k, we can implement any rotation around the complex unit circle.

Another important way we can simplify our implementation is using controlled-R gates:

These gates allow us to skip converting the second number into its Fourier transform by setting the bits of the second number as the control bits for the applied controlled-R gates. We also use these gates while converting a number into its quantum Fourier transform.

Think about what happens when we apply the gate above using a control bit that has a zero probability of being measured as 1; the target qubit would remain the same. This is the same behavior that we see when the exponent factor in the quantum Fourier transform formula is zero; the coefficient of the state vector we are dealing with remains unchanged.

Unfortunately, there is no controlled-R gate in QISkit, but since the gate is identical to the controlled-U1 gate, we can use that instead. This gate takes one parameter as input, but this time we can choose the angle by which the vector coefficient is rotated instead of a factor by which the angle 2 π divides.

The controlled-U1 gate can be applied onto the first qubit of a quantum register a in a quantum circuit qc to rotate it by θ radians, conditioned on the first qubit of a quantum register b like this:

qc.cu1(θ, a[0], b[0])

Now that we have the quantum Fourier transform down, and know what gates we can use to implement a quantum adder, we can start writing the code for one.

A faster quantum adder

The first step, as usual, is to take user input:

#Take two numbers as user input in binary form   
first = input("Enter a number with less than 10 digits.")
second = input("Enter another number with less than 10 digits.")

The reason the maximum length of the numbers is 9 is because the IBM simulator only supports up to 32 (qu)bits. If we assume that n is the length of both the numbers we want to add:

We need n + 1 qubits to hold the first number.
We need n + 1 qubits to hold the second number.
We need n + 1 classical bits to hold the final output.

In total, that makes it 3n + 3 (qu)bits, which we have to fit into 32:

3n + 3 ≤ 32

Solving for n gives n ≤ 9.67, or n ≤ 9, since we cannot have two third of a qubit.

We now have to process the input a little before we transfer the information to quantum registers:

l1 = len(first)
l2 = len(second)
#Making sure that 'first' and 'second' are of the same length 
#by padding the smaller string with zeros
if l2>l1:
    first,second = second, first
    l2, l1 = l1, l2
second = ("0")*(l1-l2) + second

We pad the smaller string to make sure it is the same length as the longer string, so that applying the quantum Fourier transform becomes easier. Now we have to define our addition procedure, which will be comprised of multiple functions.

The first function will be called createInputState(), and, as you might have guessed, converts the first number into a state appropriate for our addition process, by converting it into its quantum Fourier transform. We can perform this using a few simple rules:

Then take the index passed. Let’s call this number n.
Apply the Hadamard gate on the qubit in register a with index n.
Starting from the index of the qubit passed, count backwards until you reach zero. Let’s call this number m.
Each time you count down once, apply a controlled rotation gate, controlled by the qubit in register a with index n, on the qubit in register a with index m.
This rotation gate would be called with parameter equal to π/2^(n — m)

To understand what the function does more clearly, here is what happens when you pass the qubits of a 4 qubit register a to the function:

Quantum Fourier transform of a 4 qubit register.

The little numbers below the dots represent the angle, in radians, by which the qubit is conditionally rotated. But the controlled gates use dots on both ends, so how do we tell which qubit is the control, and which one is the target bit?

It’s simple. The qubit which has the dot with the number next to it is the control qubit, and the one without the number is the target qubit. If you compute the quantum Fourier transform matrix for a 4 qubit state, you will find that this gate sequence is equivalent. We can encode this behavior by first applying a Hadamard gate, then applying repeated controlled-U gates with argument decreasing by a factor of half with each iteration:


def createInputState(qc, reg, n, pie):
    """
    Apply one Hadamard gate to the nth qubit of the quantum register               
    reg, and then apply repeated phase rotations with parameters  
    being pi divided by increasing powers of two.
    """
    qc.h(reg[n])    
    for i in range(0, n):
        qc.cu1(pie/float(2**(i+1)), reg[n-(i+1)], reg[n])

The term phase rotation is simply another way to refer to a controlled-U gate.

Next we create a function evolveQFTstate() that performs controlled rotations on the qubits of the register holding the first number, controlled by the qubits of the second number. Writing a few rules for this function as well will help us:

Take the index passed, given by n.
Starting from the index of the qubit passed, count backwards until you hit -1, given by the number m. Before you count down once, apply a controlled rotation gate controlled by the qubit in register b with index m on the qubit in register a with index n.

We can visualize this for a four qubit register a, which holds the first number, and a register b which holds the second:

Visualization of the evolution of |F(ψ(a))> to |F(ψ(a+b))>

If we define our function according to these rules, we get this:

def evolveQFTState(qc, reg_a, reg_b, n, pie):
    """
    Evolves the state |F(ψ(reg_a))> to |F(ψ(reg_a+reg_b))> using the     
    quantum Fourier transform conditioned on the qubits of reg_b.
    Apply repeated phase rotations with parameters being pi divided 
    by increasing powers of two.
    """
    for i in range(0, n+1):
        qc.cu1(pie/float(2**(i)), reg_b[n-i], reg_a[n])

We also need a way to convert the sum a+b, stored in the register a, from its quantum Fourier transform ψ(a+b) to its output form, a+ b. To do this, we construct an inverse quantum Fourier transform function called inverseQFT().

This function would perform the same sequence performed in the createInputState() function, but in reverse order. This would involve a new set of rules:

Take the index passed, given by n.
Starting from -1, count up until you reach n. Let’s call the number you are at m.
Each time you count up once, apply a controlled rotation gate, controlled by the qubit in register a with index n, on the qubit in register a with index n-m-1.
Apply a Hadamard gate to the qubit in register a with index n.

We can visualize this for a 4 qubit register a:

Inverse quantum Fourier transform of a register.

If we write the corresponding code for our function, it would look like this:

def inverseQFT(qc, reg, n, pie):
    """
    Performs the inverse quantum Fourier transform on a register 
    reg. Apply repeated phase rotations with parameters being pi    
    divided by decreasing powers of two, and then apply a Hadamard 
    gate to the nth qubit of the register reg.
    """
    for i in range(0, n):
        qc.cu1(-1*pie/float(2**(n-i)), reg[i], reg[n])
    qc.h(reg[n])

Awesome! Now that we’ve defined our functions, we can go on to create the registers we need and store our numbers in them:

def add(first, second, n):
    pie = math.pi    a = QuantumRegister(n+1, "a") #Holds the first number
    b = QuantumRegister(n+1, "b") #Holds the second number

We then create a classical register cl to hold the final output and put all of the registers together in a quantum circuit qc:

    cl = ClassicalRegister(n+1, "cl") #Holds the final output
    qc = QuantumCircuit(a, b, cl, name="qc")

Now we have to store the two numbers in the registers a and b:

    #Flip the corresponding qubit in register a if a bit in the 
    #string first is a 1
    for i in range(0, n):
        if first[i] == "1":
            qc.x(a[n-(i+1)])    #Flip the corresponding qubit in register b if a bit in the 
    #string second is a 1
    for i in range(0, n):
        if second[i] == "1":
            qc.x(b[n-(i+1)])

Now we can call our functions:

    #Compute the Fourier transform of register a
    for i in range(0, n+1):
        createInputState(qc, a, n-i, pie)    #Add the two numbers by evolving the Fourier transform     
    #F(ψ(reg_a)) >to |F(ψ(reg_a+reg_b))>
    for i in range(0, n+1):
        evolveQFTState(qc, a, b, n-i, pie)    #Compute the inverse Fourier transform of register a
    for i in range(0, n+1):
        inverseQFT(qc, a, i, pie)

Then we measure the qubits in register a and store the results in the classical register cl:

    for i in range(0, n+1):
        qc.measure(a[i], cl[i])

Finally, we register and set our API token and url, and run our program:

    import Qconfig
    register(Qconfig.APIToken, Qconfig.config['url'])
    #Select backend and execute job
    result = execute(qc, backend='ibmq_qasm_simulator', 
                      shots=256).result()
    counts = result.get_counts("qc")
    print(counts)
    #Select result with maximum count
    output = max(counts.items(), key=operator.itemgetter(1))[0]

You might have noticed in that last article the shots parameter was set to 2, while this time it is as 512; this is to ensure that the result returned the maximum number of times is the correct one. This can be done if we select a sufficiently large number for the shots parameter, so that the most probable answer (the correct one) is also the one with the maximum count.

Why do we need to do this? More importantly, why are we getting multiple answers for the same operation?

The answer lies in a quantum mechanical property called decoherence. While the actual definition is a little more complicated, for us we can put it simply: over time, a quantum system lose its ‘quantumness’, eventually behaving akin to a classical system.

This messes up the superposition states of the system (which we use extensively) hence distorting the results we obtain slightly. Over longer runs, however, the most probable result (which is the correct one) will be the one that is returned the most.

The output I got when I ran the program with inputs 10101 (21) and 11010 (26) was 101111 (47), and can be seen below:

Enter a number with less than 10 digits. 10101
Enter another number with less than 10 digits. 11010
{'101111': 296, '001111': 216}
101111

We can also look at the gate sequence executed for this function call to understand better what operations took place:

The full program used in this tutorial can be found here.

Looking back

So, we’re done!

Admittedly, this tutorial was a little longer than it could have been — but that was only because we built a strong, intuitive understanding of the Fourier transform.

This is particularly useful as many algorithms covered later in this series utilize the Fourier transform in a similar manner as this addition implementation.

For the curious reader, quantum mechanics deals with the behavior of extremely small particles and systems like electrons and atoms, and is the basis on quantum computers are built and operate.

For a fun introduction to the topic, I recommend the book Quantum Mechanics: The Theoretical Minimum, which will leave you with a sound mathematical understanding of the principles of this exciting field.