## Introduction

Well, quantum computing for programmers anyway. I was curious as to how quantum computers worked and why they are so much better than classical computers for some problems. I turned to the fountain of all knowledge, the internet, and found only high level overviews, misinformation or very technical stuff. In the end, the most accessible thing I found was Quantum Computing for the Determined. It’s very good, but heavy going for a programmer like me. It didn’t seem like the basic concepts were all that difficult, and I still don’t think they are, so I figured the best way to try and understand this was to build a quantum computer simulator. If I could get it to run an algorithm that was exponentially faster on a quantum computer than a classical one, I might be some way to understanding quantum computing. The simulator ended up at less than a hundred (non blank) lines of python even though it’s written in a fairly verbose style. As a bonus, those hundred lines contain the implementation of a quantum algorithm that runs exponentially faster than a classical computer. Unfortunately it’s not possible to run quantum algorithms at full speed on a classical computer, so don’t be expecting to use it to break any encryption algorithms. I’m not really qualified to write this as my background is programming, so constructive criticism is most welcome. I’m not all sure how correct all this information is, but I have a reasonable degree of confidence in the simulator, as it runs a quantum algorithm and gives the expected result.

## Quantum Computing

You’ve probably heard of a Qubit. It’s the quantum equivalent of a classical bit. A classical bit can be either 1 or 0. A quantum Qubit can only ever be observed to be 1 or 0, but internally it looks quite different. Think of a clock face with only 8 hours on it, with 8 o’clock at the top. Imagine an hour hand on that clock. If the hour hand is at 4 or 8 o’clock (straight down or straight up), you’ll always observe a one on the Qubit. If the hour hand is at 2 or 6 o’clock (horizontal), you’ll always observe a zero. If the hour hand is at an odd numbered hour, it’ll be 50/50 whether you observe a 0 or a 1. This is known as a superposition of 1 and 0. Once you’ve observed a 1 or a 0 on the Qubit, you’ll always observe a 1 or 0 unless you do something to change it. In fact, it’s as if the act of looking at it moves the hour hand to 8 o’clock or 2 o’clock. We can’t ever see the internal state of a Qubit. That seems to be a fundamental part of quantum mechanics. We’ll only ever see 8 o’clock or 2 o’clock.

We can represent a single Qubit with a line pointing in a certain direction, like the hour hand on our clock.

The length of the line must always be 1 (corresponding to a probability of 1), so, from Pythagoras, the sum of the squares of x and y must be 1. The probability of a Qubit being in the 1 state is y squared, and the probability of it being in the 0 state is x squared.

## More Qubits

We’re going to represent the internal state of our quantum computer simulator as an array of numbers. We’re going to have 1 element for each possible state of the system, and it’s value will be the square root of the probability that the system is in this state. A single Qubit system has 2 possible states, 0 and 1, so it’ll be represented by a 2 element array. One element will correspond to x in the diagram, and the other will correspond to y.

Adding another Qubit doubles the number of possible states the system can be in. Think of it like this: The computer can be in all the states it was before, with the new Qubit set to 0. It can also be in all the states it was before with the new Qubit set to 1. So for a 2 Qubit system, we’d have a 4 element array, and for a 3 Qubit system, we’d have an 8 element array. For a N Qubit system we’d have a pow(2, N) element array.

A quantum computer can act on all of these states in parallel. This is where the power of quantum computers comes from, and makes simulating a quantum computer an NP complete problem. Note that it doesn’t mean quantum computers can solve NP complete problems, because we can only ever observe one of these states. Quantum algorithms can still employ some clever trickery to exploit some of this parallelism. We’ll see one of them in action in our simulator later.

The index to our state array will be the state of our computer. The square of the value of each element will be the probability that the computer is in that state. The initial state of our quantum computer must be known with certainty. This means one element of the state array will have a probability of 1, and all the other elements will have a probability of 0. The element with probability 1 represents the initial known state of the machine.

Our state array is private. This is the same in the real world. We have no way to determine the quantum state until we observe something, at which point a Qubit becomes 0 or 1 with complete certainty. For simplicity, we’re only going to allow observation of Qubits that are in a known state in our simulator. In a real quantum computer you can look at any Qubit, but if it is not in a known state, other Qubits may be affected. Consider for example:

State sqrt(probability) 00 0 01 0.5 10 0.5 11 sqrt(2)/2

State 00 and 10 are the only states with a 0 in the least significant bit. Since state 00 has 0 probability, if you happen to observe a 0 in the least significant Qubit, your computer will definitely be in state 10. So we now know that the most significant Qubit is 1. This is quantum decoherence, and the 2 qubits are said to be entangled with one another.

## Quantum Logic Gates

We’re only going to simulate one quantum logic gate in our system, as that’s all that’s required to implement a basic quantum algorithm that is exponentially faster than a classic computer. The gate is called a Hadamard gate. On a single Qubit, all it does is advance the clock by one hour and flip the hour hand. On a 2 Qubit system, our state array is already a 4 dimensional vector, so it’s not possible to visualise. A Hadamard gate applied to a Qubit in the 1 state will result in a Qubit in a state that is equally likely to be 1 or 0. This is called a superposition of 1 and 0 if you remember. I’m not sure why it flips the line. Maybe it’s so that applying it twice gets you back to your starting point, or maybe it’s just easier to implement them with this behaviour.

To apply a Hadamard gate to a single Qubit system with it’s state represented by x and y for 0 and 1 respectively, we transform x and y like this:

x = (x + y) / squareRoot(2) y = (x - y) / squareRoot(2)

To apply a Hadamard gate to a Qubit in a multi Qubit system, we need to find all the pairs of states that differ by only that Qubit, and apply the transformation above to each pair. The pairs correspond to the 0, or x, axis of the Qubit and the 1, or y, axis of the Qubit.

State Binary 0 000 1 001 2 010 3 011 4 100 5 101 6 110 7 111

For example, the pairs for the least significant bit are:

0, 1 (000, 001) 2, 3 (010, 011) 4, 5 (100, 101) 6, 7 (110, 111)

For the middle Qubit:

0, 2 (000, 010) 1, 3 (001, 011) 4, 6 (100, 110) 5, 7 (101, 111)

And for the most significant qubit:

0, 4 (000, 100) 1, 5 (001, 101) 2, 6 (010, 110) 3, 7 (011, 111)

## A Quantum Algorithm

Now all we need is a quantum algorithm to run on our simulator. We’re going to implement the Deutsch-Jozsa algorithm as it’s pretty simple (to implement, if not to understand). It’s not a very useful algorithm, but it is the basis of other quantum algorithms, such as Shor’s for prime factoring. RSA encryption relies on prime factoring being a hard problem, so Shor’s algorithm could be used to break RSA.

The Deutsch-Jozsa algorithm will take a function as it’s input. The function must take a number and map it to either 0 or 1. For every input, we’ll promise that the function will either:

- Always return 0
- Always return 1
- Return 1 or 0, but the same number of inputs will map to 0 as map to 1. For example, a function that returns 1 if the number is even, or 0 if it’s odd.

Deutsch-Jozsa will tell us if the function is constant or not. If you were to solve this problem on a classical computer, worst case you’d have to call the function for half the inputs + 1 to be sure which type of function you had. Deutsch-Jozsa will only call the function once (or rather, it will ask a quantum oracle to transform a carefully crafted superposition of possible states).

Our algorithm requires enough Qubits to represent the number passed to the function, and one extra. So if the input to our function was an 8 bit number, we’d need 9 bits for our algorithm. The extra Qubit will be the first Qubit, just for ease of implementation. This means that the even numbered indexes of the state array correspond to states where the extra bit is 0.

## The Quantum Oracle

For each state, the Quantum Oracle will call the function with all the bits of that state, except the extra bit, as input. If the function returns one, the Quantum Oracle will negate the probability root of that state.

Even though the Quantum Oracle and the Hadamard gate are implemented as a loop in our simulator, a quantum computer does all this in parallel. We have to be a bit clever to exploit this parallelism due to the restrictions on observing Qubits.

## Implementing The Quantum Algorithm

We’ll start off with our extra Qubit in the 1 state, and every other Qubit in the 0 state. Then we’ll apply the Hadamard gate to each Qubit to create a quantum superposition of every possible state. In other words, every state in which our quantum computer could be is equally likely. You can see this in the log line just before the Quantum Oracle is applied. Next we’ll apply the Quantum Oracle. The 1st Qubit is now irrelevant. Next we apply the Hadamard gate to each of the other Qubits, then we’ll observe them. If they’re all 0, our function was constant, otherwise it isn’t. I’m not going to offer a deeper understanding of how this algorithm works, as I don’t have one at the moment. Maybe I’ll spend more time figuring it out in the future, but for now, here’s the code. Use it at your own risk, but I’m pretty sure it won’t create a black hole or anything like that. It’s written in python and works on version 2.7.3 on fedora 17.

#!/usr/bin/python import math class QuantumComputerSimulator: SQRT_2 = math.sqrt(2) ONE_LOWER_TOLERANCE = 0.999 ONE_UPPER_TOLERANCE = 1.001 def __init__(self, numberOfQubits, initialState, verbose): self.__state = [0.0 for _ in range(0, pow(2, numberOfQubits))] self.__state[initialState] = 1.0 # Our quantum computer definitely starts in "initialState". self.__verbose = verbose self.log("Initial state") def applyHadamard(self, qubitIndex): "This is a constant time operation in a quantum computer." for index in range(0, len(self.__state), 2): zeroIndex = self.swapBits(0, qubitIndex, index) oneIndex = self.swapBits(0, qubitIndex, index + 1) zero = self.__state[zeroIndex] one = self.__state[oneIndex] self.__state[zeroIndex] = (zero + one) / self.SQRT_2 self.__state[oneIndex] = (zero - one) / self.SQRT_2 self.log("After Hadamard transformation on Qubit %s" % qubitIndex) def quantumOracle(self, function): "This is constant time on a quantum computer if function is constant time" for index in range(0, len(self.__state), 2): # We go in steps of 2 as the first qubit is not an input to our function result = function(index / 2) if result == 1: self.__state[index] = -self.__state[index] self.__state[index + 1] = -self.__state[index + 1] self.log("After quantum oracle") def readQubit(self, qubit): mask = 1 << qubit zeroProb = 0.0 oneProb = 0.0 for index in range(0, len(self.__state)): stateProb = self.__state[index] stateProb *= stateProb if index & mask: oneProb += stateProb else: zeroProb += stateProb if not self.isOne(oneProb + zeroProb): raise Exception("Something went wrong. Total probability for a Qubit should be 1.0. Could be a rounding error.") if self.isOne(oneProb): return True if self.isOne(zeroProb): return False # In reality, observing a qubit in a superposition (not definitely 1 or 0) would result in a 1 or 0, # and all entangled Qubits would need to be adjusted accordingly. raise Exception("Qubit is not in a known state.") def isOne(self, number): return number > self.ONE_LOWER_TOLERANCE and number < self.ONE_UPPER_TOLERANCE @staticmethod def swapBits(index1, index2, number): "Swap 2 bits in number, specified by index1 and index2" # mask for the required bits mask1 = number & (1 << index1) mask2 = number & (1 << index2) # Move them to the 1st bit bit1 = (mask1) >> index1 bit2 = (mask2) >> index2 # Zero out the bits in the number number ^= (mask1 | mask2) # Set the swapped bits number |= (bit1 << index2) number |= (bit2 << index1) return number def log(self, name): if self.__verbose: print "%s, %s" % (name, ", ".join(map(str, self.__state))) def always0(value): return 0 def always1(value): return 1 def isOdd(value): return (value & 1) def isEven(value): return (value ^ 1) & 1 def isConstantFunction(function, name, verbose): "This is the Deutsch-Jozsa algorithm" qubits = 3 computer = QuantumComputerSimulator(qubits, 1, verbose) for qubit in range(0, qubits): computer.applyHadamard(qubit) computer.quantumOracle(function) for qubit in range(1, qubits): computer.applyHadamard(qubit) functionChanges = False for qubit in range(1, qubits): functionChanges |= computer.readQubit(qubit) if verbose: print if functionChanges: print "%s is balanced" % name else: print "%s is constant" % name if verbose: print for function, name in [(always0, "Always0"), (always1, "Always1"), (isOdd, "isOdd"), (isEven, "isEven")]: isConstantFunction(function, name, True)