\documentclass[]{article} \begin{document} \title{Quantum Computing and Shor's Algorithm} \date{Original Spring 1999, Revised 2015} \author{Matthew Hayward} \maketitle \pagebreak \tableofcontents \pagebreak \section{Preface} This paper is intended to be a beginners introduction to the field of quantum computing. As I began to study this topic at the University of Illinois under the supervision of Dr\. Roy Campbell in 1999, there was a lack of introductory material to the topic. This paper will hopefully serve as an introduction to the rudiments of quantum computing and the specifics of Shor's algorithm for factoring large numbers. In order to get the most out of this paper, readers should be familiar with the following topics: \begin{itemize} \item Binary representation of numbers \item Complex numbers \item Vector mathematics \end{itemize} That being said there are appendices covering the bare essentials of those topics that are needed to fully comprehend this paper. The field of quantum computing has its own vocabulary, most of the novel terms used in this paper are described in the glossary at the end. \section{Introduction} Astounding advances have been made in the manufacture of computers between 1950 and 1999. The number of atoms needed to represent a bit in memory has been decreasing exponentially since 1950. Likewise the number of transistors per chip, clock speed, and energy dissipated per logical operation have all followed their own improving exponential trends. Despite these fantastic advances, the manner in which all computers function has been essentially identical: combinations of a few simple logical operations that process some number of operations per unit time. This rate of improvement in transistor size cannot be sustained much longer, at the current rate as of 199 in the year 2020 one bit of information will requite only one atom to represent it. The problem is that at that level of miniaturization the behavior of the components of a computer will become dominated by the principles of quantum physics. (Williams, Clearwater) With the size of components in classical computers shrinking to where the behavior of the components may soon be dominated more by quantum physics than classical physics researchers have begun investigating the potential of these quantum behaviors for computation. Surprisingly it seems that a computer whose components are able to function in a quantum are more efficient than any classical computer can be. It is the physical limitations of the classical computer, and the possibilities for the quantum computer to perform certain useful tasks more rapidly than any classical computer which drive the study of quantum computing. \section{The Classical Computer} \subsection{Turing Machines} One of the people most directly responsible for the current concept of computing machines is Alan Turing. Other early giants in the field of Computer Science include Kurt Godel, Emil Post, and Alonzo Church. Turing and others proposed mathematical models for computing which allowed for the study of algorithms and in absence of any particular computer hardware. These abstraction have proved invaluable in the field of computer science, allowing for general results about algorithms and their run time complexity to be obtained without having to consider a variety of potential computation engines on which those algorithms would run. Turing's model is called a Turing machine. The design of the Turing machine is the following: The machine consists of an infinite tape divided into cells, each cell can contain a 1, a 0, or a blank. There is also a head which can move about the tape, read a cell, write to a cell, change its internal state, or end computation. It is the internal state of this head that is the program which decides what action the head will take next. (Steane) As simple as this model appears, it encompasses all the functionality of today's most modern supercomputer. That is to say, given the time and memory, there is not a single computation that can be performed on a modern supercomputer that can not be performed on a Turing machine. The immense value of the Turing machine is that when examining what sort of functions can be computed by a computer, one need only to examine a Turing machine, and not one of the millions of potential computing devices to determine if a computation is possible. If a Turing machine can perform the computation, then it is computable. If a Turing machine can not, then the function can not be computed by any classical computer. (Shor) \subsection{Church-Turing Thesis} This bold claim, that any computer is essentially equivalent to a Turing machine grew out of contemporaneous work by Alonzo Church and Alan Turing, and is variously referred to as Church's Thesis, the Church-Turing Thesis, the Turing-Church thesis, the Church-Turing conjecture, and Turing's thesis. There are several equivalent formulations of the thesis stated at different times and to different degrees of rigor, but the central idea is that: \begin{quote} Every effectively calculable function can be produced by a Turing machine. \end{quote} The claim here is that anything which one would like to compute, which one reasonably could compute (for instance given pen and paper, and enough time), are precisely those things which a Turing machine can compute. The Church Turing thesis is perhaps best understood as a \textbf{definition} of the types of functions that are calculable in the real world - not as a theorem to be proven. As evidence for the suitability of this as a definition, multiple (indeed every one considered to far) distinct models of computation have been shown to be equivalent to the Turing model with regard to what functions they can compute. Beyond the mere ability of a computing device to calculate a function is the consideration of how long it would take to do so. This notion is generalized by the time complexity of a given algorithm, which considers the number of steps or operations an algorithm would have to perform relative to the length of the input in order to compute the function. There are various extensions of the Church-Turing thesis that address the relative efficiency of different models of computation, such as the \textbf{Complexity-Theoretic Church-Turing Thesis} states that: \begin{quote} A Turing machine can efficiently simulate any realistic model of computation. \end{quote} Here \emph{efficiently} means a Turing machine can simulate any realistic model of computation with at most polynomial time expansion of the computation time or other resources (such as memory). This hypothesis may be viewed either as a definition of what a realistic model of computation is, or as a statement of fact (which may turn out to be false) about realistic models of computation. There are no known models of computation which are thought to be realizable that violate this thesis. The discovery of a model which did violate this thesis would likely be an astonishing breakthrough in the fields of computational complexity, computer science, and/or physics. These theses gives us insight into the power of computing machines. If a computer theoretically can compute all the functions a Turing machine can compute (given enough memory and time) then it is as powerful as a Turing machine. \subsection{Running Time and Complexity Classes} Turing produced the Turing machine in an effort to answer the question ``Is there a mechanical procedure by which the truth or falsity of any mathematical conjecture can be decided?'' Turing reasoned that if there was such a method, and if it were truly mechanical in nature and requiring no creative effort, then it could be performed by a Turing machine. The proof offered by Turing that there is no such procedure is the same as the proof of Turing's Halting problem, which states that no Turing machine $H$ can in general be given any arbitrary Turing machine $X$ and its input ($I_x$) as its own ($H$'s) input, and determine whether $X$ will run to completion on $I_x$. The demonstration by Turing of a problem which can not be solved by a Turing machine suggests the division of functions into two groups: those that can be computed by a Turing machine and those that can not. By the Church-Turing thesis we can generally say that functions that are not Turing computable are simply ``not computable.'' It should be stressed at this point that this classification of a function that is not Turing computable as not computable in a general sense is due to the \textbf{definitional} aspect of the Church-Turing thesis. It is a logical possibility that some future computing machine will indeed be able to compute some function that a Turing machine can not - however no such machine is known, and the development of such a machine would be a large breakthrough. Arguments can be made that no such machine would be physically possible - however this is not a settled question. Determining if a function is computable at all is interesting, but this distinction is not fine enough to determine if a function can realistically be solved by a physical computer in a reasonable amount of time. If a computation will take billions of years to compute, it may as well be uncomputable from the perspective of the pragmatist. A algorithm can be characterized by the number of operations and amount of memory it requires to compute an answer given an input of size $n$. These characterizations of the algorithm determine what is called the algorithm's \emph{running time} or \emph{time complexity}. Specifically, the time complexity of an algorithm to compute a function is determined by looking at how the number of operations of the algorithm scale with the size of the input of the function it computes. We'll take the size of the input to be the number of bits needed to represent that number in binary. An algorithm which computes its output for an input of size $n$ using resources (computational steps, memory, etc.) that is bounded above by a polynomial function of $n$ are said to be of \emph{polynomial time}. For example if an algorithm with an input size of 10 bits took $10^{4} + 7*10^{2} + 1001$ operations to compute its answer it would be considered a polynomial time algorithm. Problems which can be solved in polynomial time or less are generally deemed to be \emph{tractable}. An algorithm which solves a problem in polynomial time may be referred to as \emph{efficient}. Problems which require more than polynomial time are usually considered to be \emph{intractable}, for example an algorithm which would require $2^{n}$ operations for an input size of $n$ would be considered intractable, as the number of operations grows exponentially with the input size, not polynomially. Generally tractable problems are considered to be those which may practically be solved in the real world. (Shor) In complexity theory, problems are often restated in terms of a \emph{decision problem} - this means that the function of interest takes in its input and produces a yes or no answer. Computer Scientists have grouped problems into a variety of complexity classes, below are some of the more well known. \begin{description} \item[P]: The class of decision problems which can be answered in polynomial time. \item[NP]: Nondeterministic polynomial time, the class of decision problems where a candidate for an answer can be verified as a correct answer or not in polynomial time. \item[NP-complete]: The ``hardest'' problems in NP, this set has the interesting property that if any NP-complete problem can be solved in polynomial time, then P and NP are in fact the same. Whether P equals NP is an outstanding problem in computer science and complexity theory. \end{description} (Cormen, Leiserson, Rivest) Given these distinctions, to determine if an function may be \emph{efficiently} computed by a classical computer there are two questions that must be answered. First, can a Turing machine compute the function, if so then the function is called computable. Second the time complexity of the algorithm to be used must by considered. In general if an algorithm requires polynomial time or less to compute it is considered to be tractable, and if not it is considered to be intractable. Accordingly, interest in quantum computing is twofold. First it is of interest if a quantum computer can compute functions which are uncomputable on a classical computer, and second it is of interest whether an algorithm which is intractable on a Turing machine may be tractable on a quantum computer. \section{The Quantum Computer} \subsection{Quantum Physics} When considering the possible efficiency of a computer which makes use of the results of quantum physics, it is helpful to know a little of quantum physics itself. Quantum physics arose from the failure of classical physics to offer correct predictions on the behavior of photons and other elementary particles. Quantum physics has since been under intense scrutiny, as some of its predictions seem very strange indeed. Nevertheless experiments verify the same strange behavior which leads skeptics to challenge the veracity of Quantum physics. \subsection{The Classical Bit} To understand the ways in which a quantum computer is different from a classical computer you must first understand the rudiments of the classical computer. The most fundamental building block of a classical computer is the bit. A bit is capable of storing one piece of information, it can have a value of either 0 or 1. Any amount of information can be encoded into a list of bits. In a classical computer a bit is typically stored in a silicone chip, or on a metal hard drive platter, or on a magnetic tape. About $10^{10}$ atoms were used to represent one bit of information in 1999. The smallest conceivable storage for a bit involves a single elementary particle of some sort. For example, for any particle with a spin-1/2 characteristic (such as a proton, neutron, or electron), it's spin-1/2 characteristic on measurement will be either $+1/2$ or $-1/2$. We can thus encode a bit using a single particle by mapping 1 to be spin $+1/2$ and 0 to be $-1/2$, assuming we can measure and manipulate the spin of such a particle. If we were to try to use this spin-1/2 particle as a classical bit, one that is always in the 0 or 1 state, we would fail. We would be trying to apply classical physics on a scale where it simply is not applicable. This single spin-1/2 particle will instead act in a quantum manner. (Williams, Clearwater) This spin-1/2 particle which behaves in a quantum manner could be the fundamental building block of a Quantum computer. We could call it a qubit, to denote that it is analogous in some ways to a bit in a classical computer. Just as a memory register in a classical computer is an array of bits, a quantum memory register is composed of several spin-1/2 particles, or qubits. There is no particular need for the spin-1/2 particle, equally well we could use a Hydrogen atom, and designate its electron being measured in the ground state to be the 0 state, and it being in the first excited state to be the 1 state. There are a multitude of possible physical qubit representations that could work. For simplicity I will discus only the spin-1/2 particle from here on, but analogous arguments could be made for many things. \subsection{State Vectors and Dirac Notation} We wish to know exactly how the behavior of the spin-1/2 particle, our qubit, differs from a that of a classical bit. Recall that a classical bit can store either a 1 or a 0, and when measured the value observed will always be the value stored. Quantum physics states that when we measure the spin-1/2 particles state we will determine that it is in the $+1/2$ state, or the spin $-1/2$ state. In this manner our qubit is not different from a classical bit, for it can be measured to be in the $+1/2$, or 1 state, or the $-1/2$, or 0 state. The differences between the qubit and the bit come from what sort of information a qubit can store when it is not being measured. According to quantum physics we may describe the state of this spin-1/2 particle by a state vector in a Hilbert Space. A Hilbert Space is a linear, complex vector space. In a linear vector space you may add and multiply vectors that lie within the space together and the resulting vector will still lie within that space. For example you are probably familiar with the $x$, $y$, $z$ coordinate system where the $x$, $y$, and $z$ axes are mutually perpendicular real number lines, which coincide at the point $x = 0$, $y = 0$, $z = 0$ - this is an example of a linear vector space. In a complex vector space the lengths of the vectors within the space are described with complex numbers. Complex numbers are numbers which take the form $a + i*b$, where $a$ and $b$ are real numbers, and $i$ is defined to be the square root of negative one. (Williams, Clearwater) In the Hilbert Space for our state vector, which describes the state of our spin-1/2 particle, these perpendicular axes will correspond to each possible state that the system can be measured in. So our Hilbert Space for a single qubit will have two perpendicular axes, one corresponding to the spin-1/2 particle being in the $+1/2$ state, and the other to the particle being in the $-1/2$ state. These states which the vector can be measured to be are referred to as \emph{eigenstates}. The vector which exists somewhere in this space which represents the state of our spin-1/2 particle is called the \emph{state vector}. The projection of the state vector onto one of the axes shows the contribution of that axes' eigenstate to the whole state. This means that in general, the state of the spin-1/2 particle can be any combination of the base states. In this manner a qubit is totally unlike a bit, for a bit can exist in only the 0 or 1 state, but the qubit can exist, in principle, in any combination of the 0 and 1 state, and is only constrained to be in the 0 or 1 state when we measure the state. Now I will introduce the standard notation for state vectors in Quantum physics. The state vector is written with the following way and called a \emph{ket vector} $|\psi \rangle$. Where $\psi$ is a list of numbers which contain information about the projection of the state vector onto its base states. The term ket and this notation come from the physicist Paul Dirac who wanted a concise shorthand for writing formulas that occur in Quantum physics. These formulas frequently took the form of the product of a row vector with a column vector. Thus he referred to row vectors as \emph{bra vectors} represented as $\langle y|$. The product of a bra and a ket vector would be written $\langle y|x \rangle$, and would be referred to as a \emph{bracket}. (Williams, Clearwater) There is nothing special about this vector notation, you may think of any state vector being written as a letter with a line over it, or as a bold letter as vectors are normally denoted. If you do further reading into this area you will almost assuredly come across the bra and ket notation, which is why it is presented. \subsection{Superposition and Eigenstates} Earlier I said that the projection of the state vector onto one of the perpendicular axes of its Hilbert Space shows the contribution of that axes' eigenstate to the whole state. You may wonder what is meant by the \emph{whole state}. You would think (and rightly so, according to classical physics) that our spin-1/2 particle could only exist entirely in one of the possible $+1/2$ and $-1/2$ states, and accordingly that its state vector could only exist lying completely along one of its coordinate axes. If the particle's axes are called $x$ and $y$, and the state vector's x coordinate which denotes the contribution to the $-1/2$, or 0 state, and y coordinate which denotes the contribution to the $+1/2$, or 1 state, should only be (1,0), or (0,1). That seems perfectly reasonable, but it simply is not correct. According to quantum physics a quantum system can exist in a mix of all of its allowed states simultaneously. This is the Principle of Superposition, and it is key to the power of the quantum computer. While the physics of superposition is not simple at all, mathematically it is not difficult to characterize this kind of behavior. \subsection{The Quantum qubit} Back to our qubit, our spin-1/2 particle. Now that we know that while it can only be measured to have a spin of $+1/2$ or $-1/2$, it may in general be in a superposition of these states when we are not measuring it. We could refer to its state in the following manner. Let $x_{1}$ be the eigenstate corresponding to the spin $+1/2$ state, and let $x_{0}$ be the eigenstate corresponding to the spin $-1/2$ state. Let $X$ be the total state of our state vector, and let $w_{1}$ and $w_{0}$ be the complex numbers that weight the contribution of the base states to our total state, then in general: \[|X \rangle = w_{0} * |x_{0} \rangle + w_{1} * |x_{1} \rangle \equiv (w_{0},w_{1})\] At this point it should be remembered that $w_{0}$ and $w_{1}$, the weighting factors of the base states are complex numbers, and that when the state of $X$ is measured, we are guaranteed to find it to be in either the state: \[0 * |x_{0} \rangle + w_{1} * |x_{1} \rangle \equiv (0,w_{1})\] or the state \[w_{0} * |x_{0} \rangle + 0 * |x_{1} \rangle \equiv (w_{0},0)\] This is analogous to a system you may be more familiar with, a vector with real weighting factors in the a two dimensional plane. Let the base states for this two dimensional plane be the unit vectors $x$, and $y$. In this case we know that the state of any vector $V$ can be described in the following manner: \[V = x_{0} * x + y_{0} * y \equiv (x_{0},y_{0})\] Our state vector is a unit vector in a Hilbert space, which is similar to vector spaces you may be more familiar with, but it differs in that the lengths of the vectors are complex numbers. It is not necessary from a physics perspective for the state vector to be a unit vector (by which I mean it has a length of 1), but it makes for easier calculations further on, so I will assume from here on out that the state vector has length 1. This assumption does not invalidate any claims about the behavior of the state vector. To see how to convert a state vector of any length to length 1 see appendix A. With all this in mind we have fully defined the basic building block of our quantum computer, the qubit. It is fundamentally different from our classical bit in that it can exist in any superposition of the 0 and 1 states when it is not being measured. (Barenco, Ekert, Sanpera, Machiavello) \subsection{The Quantum Memory Register} Up till now we have considered a two state quantum system, specifically our spin-1/2 particle. However a quantum system is by no means constrained to be a two state system. Much of the above discussion for a 2 state quantum system is applicable to a general $n$ state quantum system. In an $n$ state system our Hilbert Space has $n$ perpendicular axes, or eigenstates, which represent the possible states that the system can be measured in. As with the two state system, when we measure our $n$ state quantum system, we will always find it to be in exactly one of the $n$ states, and not a superposition of the $n$ states. The system is still allowed to exist in any superposition of the $n$ states while it is not being measured. Mathematically if two state quantum system with coordinate axes $x_{0}$, $x_{1}$ can be fully described by: \[|X \rangle = w_{0} * |x_{0} \rangle + w_{1} * |x_{1} \rangle \equiv (w_{0},w_{1})\] Then an $n$ state quantum system with coordinate axes $x_{0}, x_{1},\ldots, x_{n-1}$ can be fully described by: \[|X \rangle = \sum_{k = 0}^{n-1} w_{k} * |x_{k} \rangle \] In general a quantum system with $n$ base states can be represented by the $n$ complex numbers $w_{0}$ to $w_{n-1}$. When this is done the state may be written as: \[ |X \rangle = \left( \begin{array}{c} w_{0} \\ w_{1} \\ \vdots \\ w_{n-1} \end{array} \right) \] Where it is understood that $w_{k}$ refers to the complex weighting factor for the $k$'th eigenstate. Using this information we can construct a quantum memory register out of the qubits described in the previous section. We may store any number $n$ in our quantum memory register as long as we have enough qubits, just as we may store any number $n$ in a classical register as long as we have enough classical bits to represent that number. The state of the quantum register with $n$ states is give by the formula above. Note that in general a quantum register composed of $n$ qubits requires $2^{n}$ complex numbers to completely describe its state. A $n$ qubit register can be measured to be in one of $2^{n}$ states, and each state requires one complex number to represent the projection of that total state onto that state. In contrast a classical register composed of $n$ bits requires only $n$ integers to fully describe its state. This means that one can store an exponential amount of information in a quantum register relative to the size of that register. Here we see some of the first hints that a quantum computer could be exponentially more efficient than a classical computer in some respects. Recall that from our discussion of complexity that problems which can be solved in polynomial time are generally thought of as being tractable, and that problems which can be solved in exponential time are thought of as intractable. If a quantum computer can be is exponentially more efficient than a classical computer for some algorithms, then some problems currently intractable may become tractable! This is a large part of the motivation for the study of quantum computing. \subsection{Probability Interpretation} Now that we know how to represent our state vector as a superposition of states, and we know that we can only measure the state vector to be in one of the base states, it would seem that there would be some sort of discrepancy. We must determine what happens when we measure the state vector. We know from quantum physics that given an initial condition the state vector will evolve in time in accordance with Schr\"{o}denger's equation: \[ih\frac{\partial |X(t) \rangle }{\partial t} = H(t)|X(t) \rangle\] Where $i$ is the square root of negative one, $h$ is $1.0545*10^{-34}$ Js, and $H$ is the Hamiltonian operator, which is determined by the physical characteristics of the system being evolved. In our notation this expression is: \[ih\frac{\partial w_{i}(t)}{\partial t} = \sum_{j}H(t)_{ij} * w_{j}(t)\] This evolution would appear to be continuous, but these equations only apply to a quantum mechanical system evolving in isolation from the environment. The only way to observe the state of the state vector is to in some way cause the quantum mechanical system to interact with the environment. When the state vector is observed it makes a sudden discontinuous jump to one of the eigenstates, but this is not a violation of Schr\"{o}denger's equation. When the quantum mechanical system interacts with the outside environment and measured the state vector is said to have collapsed. (Williams, Clearwater) Now understanding that a state vector will collapse when it interacts with the external environment, we still need to know in what manner this collapse happens. To perform any sort of useful calculation we must be able so say something about which base state a quantum mechanical system will collapse into. The probability that the state vector will collapse into the $j$'th eigenstate, is given by $|w_{j}|^{2}$ which is defined to be $a_{j}^{2} + b_{j}^{2}$ if $w_{j} = a_{j} + i*b_{j}$, where $w_{j}$ is the complex projection of the state vector onto the $j$'th eigenstate. In general the chance of choosing any given state is \[Prob(j) = \frac{|w_{j}|^{2}}{\sum_{k = 0}^{ n-1} |w_{k}|^{2}}\] but as mentioned earlier we will insist on having the state vector of length one, and in this case the probability expression simplifies to $Prob(j) = |w_{j}|^{2}$. So now we know how to construct an $n$ state quantum system, which can be placed in an arbitrary superposition of states. We also know how to measure the resultant superposition and get a base state with a certain probability. This is all that we need to understand about our quantum memory register to be able to simulate its behavior. \subsection{Quantum Parallelism} Given that our quantum memory register differs from a classical one in that it can store a superposition of the base states of the register, one might wonder what this implies as to the efficiency of quantum computing. The study of quantum computing is relatively new, most give credit to Richard Feynman for being the first to suggest that there were tasks that a quantum computer could perform exponentially better than a classical computer. Feynman observed that a classical computer could not simulate a quantum mechanical system without suffering from exponential slowdown. At the same time he hinted that perhaps by using a device whose behavior was inherently quantum in nature one could simulate such a system without this exponential slowdown. (Feynman) Several quantum algorithms rely on something called quantum parallelism. Quantum parallelism arises from the ability of a quantum memory register to exist in a superposition of base states. A quantum memory register can exist in a superposition of states, each component of this superposition may be thought of as a single argument to a function. A function performed on the register in a superposition of states is thus performed on each of the components of the superposition, but this function is only applied one time. Since the number of possible states is $2^{n}$ where $n$ is the number of qubits in the quantum register, you can perform in one operation on a quantum computer what would take an exponential number of operations on a classical computer. This is fantastic, but the more superposed states that exist in you register, the smaller the probability that you will measure any particular one will become. As an example suppose that you are using a quantum computer to calculate the function $\mathcal{F}(x) = 2*x \bmod 7$, for $x$ integers between 0 and 7 inclusive. You could prepare a quantum register that was in a equally weighted superposition of the states 0-7. Then you could perform the $2*x \bmod 7$ operation once, and the register would contain the equally weighted superposition of 1,2,4,6,1,3,5,0 states, these being the outputs of the function $2*x \bmod 7$ for inputs 0 - 7. When measuring the quantum register you would have a $2/8$ chance of measuring 1, and a $1/8$ chance of measuring any of the other outputs. It would seem that this sort of parallelism is not useful, as the more we benefit from parallelism the less likely we are to measure a value of a function for a particular input. Some clever algorithms have been devised, most notably by Peter Shor and L. K. Grover which succeed in using quantum parallelism on a function where they are interested in some property of all the inputs, not just a particular one. \subsection{The Efficiency of the Quantum Computer} Given the possible efficiency of quantum parallelism much work has been done to show formally with mathematical proofs how quantum computers differ from classical ones in the complexity classes and running times of various algorithms. Here is a short list of some of the landmarks in the study of the efficiency of quantum computers. In 1980 Paul Benioff offered a classical Turing machine which used quantum mechanics in its workings, thus showing that theoretically a quantum computer was at least as powerful as a classical computer. (Benioff) In 1982 Richard Feynman showed that a classical Turing Machine (and hence any classical computer if the Complexity-Theoretic Church-Turing Thesis holds) could not simulate a quantum mechanical system without suffering exponential slowdown. (Feynman) David Deutsch and Richard Jozsa showed in a paper in 1992 that there was an algorithm that could be run in poly-log time on a quantum computer, but required linear time on a deterministic Turing machine. This may have been the first example of a quantum computer being shown to be exponentially faster than a deterministic Turing machine. However, the problem could also be solved in poly-log time in a probabilistic Turing machine, a Turing machine which is capable of making a random choice. (Deutsch, Jozsa) Andre Berthiaume and Giles Brassard proved in a 1992 paper that $P \subset QP$, where $P$ is a complexity class as mentioned earlier and $QP$ (also denoted as $EQP$) corresponds to problems which can be solved in worst case polynomial time by a quantum computer, so there are indeed problems which can be solved in polynomial time on a quantum computer that can not be solved in a polynomial time with a deterministic Turing machine. (Berthiaume, Brassard) \section{Quantum Algorithms} \subsection{Introduction to Shor's Algorithm} By the early nineties it was known that a quantum computer could be more efficient than any classical computer for certain tasks of the Complexity-Theoretic Church-Turing thesis holds . Nonetheless investigations into these areas were largely driven by academic curiosity. There was not much economic motive for people to spend lots of money or time trying to build a quantum computer. That changed in 1994 when Peter Shor, a scientist working for Bell Labs, devised a polynomial time algorithm for factoring large numbers on a quantum computer. This discovery drew great attention to the field of quantum computing. \subsection{Motivation for Shor's Algorithm} The algorithm was viewed as important because the difficulty of factoring large numbers is relied upon for many cryptography systems. If an efficient method of factoring large numbers is implemented most of the many encryption schemes would be next to worthless to protect their data from prying eyes. While it has not been proven that factoring large numbers can not be archived on a classical computer in polynomial time, as of 2015 the fastest algorithm publicly available for factoring large number runs in $O(exp(\frac{64}{9}n^{1/3} (log\ n)^{2/3})$, operations where $n$ is the number of bits used to represent the number: this runtime exceeds polynomial time. In contrast Shor's algorithm runs in $O((\log n)^{2} * \log \log n)$ on a quantum computer, and then must perform $O(\log n)$ steps of post processing on a classical computer. Overall then this time is polynomial. This discovery propelled the study of quantum computing forward, as such an algorithm is much sought after. (Shor) \subsection{Overview of Shor's Algorithm} Shor's algorithm hinges on a result from number theory. This result is: \begin{quote} The function $\mathcal{F}(a) = x^{a} \bmod n$ is a periodic function, where $x$ is an integer coprime to $n$. In the context of Shor's algorithm $n$ will be the number we wish to factor. When two numbers are coprime it means that their greatest common divisor is 1. \end{quote} Calculating this function for an exponential number of $a$'s would take exponential time on a classical computer. Shor's algorithm utilizes quantum parallelism to perform the exponential number of operations in one step. The reason why this function is of utility in factoring large numbers is this: \begin{quote} Since $\mathcal{F}(a)$ is a periodic function, it has some period $r$. We know that $x^{0} \bmod n = 1$, so $x^{r} \bmod n = 1$, and $x^{2r} \bmod n = 1$ and so on since the function is periodic. \end{quote} Given this information and through the following algebraic manipulation: \[x^{r} \equiv 1 \bmod n\] \[(x^{r/2})^{2} = x^{r} \equiv 1 \bmod n\] \[(x^{r/2})^{2} - 1 \equiv 0 \bmod n\] and if $r$ is an even number \[(x^{r/2} - 1)(x^{r/2} + 1) \equiv 0 \bmod n\] We can see that the product $(x^{r/2} - 1)(x^{r/2} + 1)$ is an integer multiple of $n$, the number to be factored. So long as $x^{r/2}$ is not equal to $\pm 1$, then at least one of $(x^{r/2} - 1)$ or $(x^{r/2} + 1)$ must have a nontrivial factor in common with $n$. So by computing $\gcd(x^{r/2} - 1, n)$, and $\gcd(x^{r/2} + 1, n)$, we will obtain a factor of $n$, where $\gcd$ is the greatest common denominator function. Here is a brief overview of Shor's algorithm, which is explained in detail in the next section. Shor's algorithm tries to find $r$, the period of $x^{a} \bmod n$, where $n$ is the number to be factored and $x$ is an integer coprime to $n$. To do this Shor's algorithm creates a quantum memory register with two parts. In the first part the algorithm places a superposition of the integers which are to be $a$'s in the $x^{a} \bmod n$ function. We will choose our $a$'s to be the integers 0 through $q - 1$, where $q$ is the power of two such that $n^{2} \leq q < 2n^{2}$. Then the algorithm calculates $x^{a} \bmod n$, where $a$ is the superposition of the states, and places the result in the second part of the quantum memory register. Next the algorithm measures the state of the second register, the one that contains the superposition of all possible outcomes for $x^{a} \bmod n$. Measuring this register has the effect of collapsing the state into some observed value, say $k$. It also has the side effect of projecting the first part of the quantum register into a state consistent with the value measured in the second part. Since we have partitioned our quantum register into two parts measurement of the second part collapses that part into exactly one value, while the other partition collapses into a state consistent with the observed value in the other portion of the register. It is still possible for the non measured part of the register to exist in a superposition of base states, as long as each of those base states are consistent with the measured value in the other part of the register. What this means in this instance is that after this measurement the second part of the register contains the value $k$, and the first part of the register contains a superposition of the base states which when plugged into $x^{a} \bmod n$ produce $k$. Since we know $x^{a} \bmod n$ is a periodic function, we know that the first part of the register will contain the values $c$, $c+r$, $c+2r\ldots$ and so on, where $c$ is the lowest integer such that $x^{c} \bmod n = k$. The next step is to perform a discrete Fourier transform on the contents of first part of the register (the one containing all integers such that $x^{a} \bmod n = k$) and to put the result back into register one. The application of the discrete Fourier transformation has the effect of peaking the probability amplitudes of the first part of the register at integer multiples of the quantity $q / r$. Now measuring the first part of the quantum register will yield an integer multiple of the inverse period. Once this number is retrieved from the quantum memory register, a classical computer can do some analysis of this number, make a guess as to the actual value of $r$, and from that compute the possible factors of $n$. This post processing will be covered in more detail later. (Shor) \subsection{Steps to Shor's Algorithm} Shor's algorithm for factoring a given integer $n$ can be broken into some simple steps. \begin{description} \item[Step 1] Determine if the number $n$ is a prime, a even number, or an integer power of a prime number. If it is we will not use Shor's algorithm. There are efficient classical methods for determining if a integer $n$ belongs to one of the above groups, and providing factors for it if it is. This step would be performed on a classical computer. \item[Step 2] Pick a integer $q$ that is the power of 2 such that $n^{2} \leq q < 2n^{2}$. This step would be done on a classical computer. \item[Step 3] Pick a random integer $x$ that is coprime to $n$. When two numbers are coprime it means that their greatest common divisor is 1. There are efficient classical methods for picking such an $x$. This step would be done on a classical computer. \item[Step 4] Create a quantum register and partition it into two parts, register 1 and register 2. Thus the state of our quantum computer can be given by: $|$reg1, reg2$ \rangle$. Register 1 must have enough qubits to represent integers as large as $q - 1$. Register 2 must have enough qubits to represent integers as large as $n - 1$. The calculations for how many qubits are needed would be done on a classical computer. \item[Step 5] Load register 1 with an equally weighted superposition of all integers from 0 to $q - 1$. Load register 2 with all zeros. This operation would be performed by our quantum computer. The total state of the quantum memory register at this point is: \[\frac{1}{\sqrt{q}}\sum_{a = 0}^{q - 1} |a, 0 \rangle\] \item[Step 6] Now apply the transformation $x^{a} \bmod n$ to for each number stored in register 1 and store the result in register 2. Due to quantum parallelism this will take only one step, as the quantum computer will only calculate $x^{|a \rangle } \bmod n$, where $|a \rangle $ is the superposition of states created in step 5. This step is performed on the quantum computer. The state of the quantum memory register at this point is: \[\frac{1}{\sqrt{q}}\sum_{a=0}^{q-1}|a,x^{a} \bmod n \rangle \] \item[Step 7] Measure the second register, and observe some value $k$. This has the side effect of collapsing register one into a equal superposition of each value a between 0 and $q-1$ such that \[x^{a} \bmod n = k\] This operation is performed by the quantum computer. The state of the quantum memory register after this step is: \[\frac{1}{\sqrt{||A||}}\sum_{a'=a'\in A}|a', k \rangle \] Where $A$ is the set of $a$'s such that $x^{a} \bmod n = k$, and $||A||$ is the number of elements in that set. \item[Step 8] Now compute the discrete Fourier transform on register one. The discrete Fourier transform when applied to a state $|a \rangle $ changes it in the following manner: \[|a \rangle = \frac{1}{\sqrt{q}}\sum_{c=0}^{q-1}|c \rangle *e^{2\pi iac/q}\] This step is performed by the quantum computer in one step through quantum parallelism. After the discrete Fourier transform our register is in the state: \[\frac{1}{\sqrt{||A||}}\sum_{a'\in A}\frac{1}{\sqrt{q}}\sum_{c=0}^{q-1} |c,k \rangle *e^{2\pi ia'c/q}\] \item[Step 9] Measure the state of register one, call this value $m$, this integer $m$ has a very high probability of being a multiple of $q/r$, where $r$ is the desired period. This step is performed by the quantum computer. \item[Step 10] Take the value $m$, and on a classical computer do some post processing which calculates $r$ based on knowledge of $m$ and $q$. In particular: \begin{itemize} \item $m$ has a high probability of being = $\lambda * (q / r)$ where $\lambda$ is an integer \item If we perform floating point division of $m / q$, and then calculate the best rational approximation to $m / q$ whose denominator is less than or equal to $q$ \item We take this denominator to be a candidate for $r$. \item If our candidate $r$ is odd we either double it if doing so leads to a value less than $q$ \end{itemize} There are efficient ways to do this post processing on a classical computer. \item[Step 11] Once you have attained $r$, a factor of $n$ can be determined by taking $\gcd(x^{r/2} + 1, n)$ and $\gcd(x^{r/2} - 1, n)$. If you have found a factor of $n$, then stop, if not go to step 4. This final step is done on a classical computer. \end{description} Step 11 contains a provision for what to do if Shor's algorithm failed to produce factors of $n$. There are a few reasons why Shor's algorithm can fail, for example the discrete Fourier transform could be measured to be 0 in step 9, making the post processing in step 10 impossible. The algorithm will sometimes find factors of 1 and $n$, which is not useful either. For these reasons step 11 must be able to jump back to step 4 to start over. (Williams, Clearwater) \subsection{Other Quantum Algorithms} Shor's algorithm is not the only algorithm that seems to be better on a quantum computer than any classical computer for a problem which is considered to be useful. In 1994 L. K. Grover, also of Bell Labs, devised an algorithm to find an item in an unsorted list of $N$ items in $.758\sqrt{N}$ operations. No classical algorithm can guarantee finding the item in less than $N$ operations, and in the average case it would take $N/2$ operations. (Grover) \section{A Simulation of Shor's Algorithm on a Classical Computer} In order to make the steps of Shor's algorithm more concrete to me I made a program which simulates the operation of a quantum computer performing Shor's algorithm. As noted by Feynman and others it is impossible to simulate a quantum mechanical system on a classical computing device without incurring exponential slowdown, thus the simulation is only tolerable to observe for relatively small numbers. The simulation implements all of the steps in Shor's algorithm, and successfully factors numbers. The code is written in C++ and can be found in Appendix D. \subsection{Introduction to the Code for the Simulation} The code base for the simulation of Shor's algorithm consists of four files. \begin{description} \item[complex.C]: This file contains the simple complex number class that I wrote for storing state information about my quantum memory register in the simulation. \item[qureg.C]: This file contains the quantum register class that simulates the behavior of the quantum memory register in Shor's algorithm. It is generic enough that it may be made to simulate any quantum memory register, although what happens in the event of a collapse not due to the direct measurement of the register is left to the programmer to perform. An example of such a collapse that is the collapse of the first partition of Shor's quantum memory register in step 7 of his algorithm. \item[shor.C]: This is the main body of the simulation, which contains the Steps of Shor's algorithm as described above. \item[util.C]: This is a library of useful functions used by shor.C \item[Makefile]: A simple makefile to produce the shor-simulator executable. \end{description} \subsection{The Complex Number Class} The complex class is needed because the probability amplitudes, or the lengths, of a vector in a Hilbert Space are in general complex. The state of our quantum memory register may be thought of as a vector in a Hilbert Space. Each complex number can be set to a complex value, which is stored internally as a real part, and an imaginary part, both of which are double precision floating point numbers in this simulation. \subsection{The Quantum Memory Register Class} The qureg class, my quantum memory register, is defined in qureg.C. When a quantum memory register is constructed it takes as input the number of qubits which it will contain. In my implementation quantum memory register of $n$ bits requires $2^{n}$ complex numbers to represent. This is because a register of $n$ bits can exist in any one of $2^{n}$ base states, and a quantum register may exist in any superposition of those base states. I use one complex number per base state to describe the probability of that state being measured. For example, if a quantum register has 3 bits it can exist in any one of the following states: \[|0,0,0 \rangle \] \[|0,0,1 \rangle \] \[|0,1,0 \rangle \] \[|0,1,1 \rangle \] \[|1,0,0 \rangle \] \[|1,0,1 \rangle \] \[|1,1,0 \rangle \] \[|1,1,1 \rangle \] If I take these values to be binary numbers where the most significant bit is the left most bit then these states correspond to the numbers 0,1,2,3,4,5,6,7. There are $2^{3} = 8 $possible base states. I follow the convention that the probability of measuring the $j$'th state whose complex amplitude is $w_{j}$ is $|w_{j}|^{2}/\sum_{j}|w_{j}|^{2}$. I also ensure throughout that algorithm that the $\sum_{j}|w_{j}|^{2}$ are equal to 1, so the probability of measuring the $j$'th state is $|w_{j}|^{2}$. Since my qureg class is only being used in Shor's algorithm, any given base state that is measured will correspond to an integer value, for example, in steps 7 and 9 we will measure only integer values. These the base states of our $n$ bit quantum register can be thought of as being the integers 0 through $2^{n - 1}$. In my code a quantum register of size $n$ will be represented by an array of complex numbers of size $2^{n}$. The value stored at the $j$'th array position is the complex probability amplitude associated with measuring the number $j$. For example, in a quantum memory register of size 2 that contains an equal superposition of the numbers 0-3 would be represented by: \[State[0] = 1/2 + i * 0 \] \[State[1] = 1/2 + i * 0 \] \[State[2] = 1/2 + i * 0\] \[State[3] = 1/2 + i * 0\] Observe that the probability of measuring any given state is $|w_{j}|^{2} = \frac{1}{2}^{2} + 0^{2} = 1/4$. If we attempt to measure this register, we will with equal probability get the number 0,1,2,3. Let us assume we measure the state and get the value 2, this has the effect of collapsing all probabilities not measured to 0, so the new state of our quantum register is: \[State[0] = 0 + i * 0\] \[State[1] = 0 + i * 0\] \[State[2] = 1 + i * 0\] \[State[3] = 0 + i * 0\] If we were to measure this register again we would measure the value 2 without fail. Sometimes during the course of operation we cause a state to exist in the quantum memory register where the sum over the squares of the probability amplitudes is not 1, for example a quantum register could be in the state: \[State[0] = 1 + i * 0\] \[State[1] = 1 + i * 0\] We would expect that when measured 0 and 1 would be equally likely outcomes, however we would also like for simplicity to have $|w_{j}|^{2}$ be the probability of measuring the $j$'th state. To accommodate this desire the qureg class contains a subroutine which will normalize the probability amplitudes, such that sum of the $|w_{j}|^{2}$ over all $j$'s is one. The qureg class also allows the state of the register to be set to any arbitrary state. This functionality is used in steps 5, 6, 7, and 8 of Shor's algorithm as described above. The quantum memory register class has two member functions which do things that a real quantum memory register could not do. The first such function is the Dump() function which dumps the entire state information of the register without collapsing it, and is included for debugging. The reason why a real quantum memory register could not do such a thing is that to detect the state information of a real quantum memory register without applying some measurement, which would collapse the state vector of the quantum memory register into one of the base states. The second such function is the GetProb(state) which gets the probability amplitude of any given state. This function suffers from the same problems as the Dump function, it is impossible to attain information about the state vector without causing it to collapse into one of the base states. The GetProb function is used in the main program to calculate the discrete Fourier transform in step 8. While a real quantum computer could not use this information in step 8, a real quantum computer would not need to. In a real quantum computer the Fourier transformation would take place in one step, and the transform function would take as its argument a state which could in general be any superposition of base states. This function is used to simulate this use of a function which takes a superposition of states as its argument. \subsection{The Simulation of Shor's Algorithm} The implementation of Shor's algorithm found in shor.C follows the steps outlined in the description of Shor's algorithm found above. There are some significant differences in the behavior of the simulator and the behavior of a actual quantum computer. The simulator uses $2^{n+1}$ double precision floating point numbers to represent the state of the quantum register. It uses one Complex object to represent the probability amplitude of each of the $2^{n}$ eigenstates of a $n$ bit quantum memory register, and each Complex object uses two double precision floating point numbers. On a machine in which a double precision floating point number is represented with 64 bits the simulator uses approximately $2^{n+7}$ classical bits used to represent the state of the quantum memory register, where $n$ is the number of bits required to represent the number the simulator is trying to factor. A real quantum computer would use exactly $n$ qubits as its memory register. This exponential space usage is apparently unavoidable, as a classical computer cannot simply cannot simulate a quantum computer without suffering from exponential time or space overhead if the Complexity-Theoretic Church-Turing Thesis holds. A second large difference is that during the modular exponentiation step of Shor's algorithm (Step 6 above) a quantum computer would perform one operation of $x^{a} \bmod n$, where a is the superposition of states in register 1. The simulation must calculate the superposition of values caused by calculating $x^{a} \bmod n$ for $a = 0$ through $q - 1$ iteratively. The simulation also stores the result of each modular exponentiation, and uses that information to collapse register 1 in step 7 in Shor's algorithm. A quantum computer would not be capable of performing this book keeping, as examining any particular result would collapse the existing superposition to be placed in register 2 at the end of step 6 of Shor's algorithm. A quantum computer would have no need whatsoever to do such bookkeeping, as when when register 2 is measured in step 7, the collapse of register 1 is an automatic and unavoidable consequence of the measurement. A analogous argument follows for the use of the get probability function in step 8 of Shor's algorithm for calculating the discrete Fourier transform. Aside from these differences, which are necessitated by the inability of a classical computer to accurately depict the behavior of a quantum mechanical system, the operations are performed by the shor.C program are identical to those called for in the description of Shor's algorithm. As the algorithm runs the state of the quantum memory register changes in the manner laid out in the description of Shor's algorithm. After the final measurement of register 1 in step 9 we obtain some integer $m$, which has a high probability of being an integer multiple $\lambda$ of $q / r$, where $\lambda$ is some integer, and $r$ is the period of $x^{a} \bmod n$ that we are trying to find, so that we may calculate $x^{r/2} - 1 \bmod n$ and $x^{r/2} + 1 \bmod n$ in an effort to find numbers which share factors with $n$. In the post processing step shor.C takes the integer $m$ and divides it by $q$, thus yielding some number $c$ which is approximately equal to $\lambda / r$, where $\lambda$ is an integer, and $r$ is the desired period. Then using a helper function it calculates the best rational approximation to $c$ which has a denominator that is less than or equal to $q$ (recall that $q$ is the power of two such that $n^{2} \leq q < 2n^{2}$, where $n$ is the number to be factored by Shor's algorithm). We take this denominator to be our period. Shor's algorithm can only use even periods in determining factors of $n$, and so we check to see if $r$ is even, if not we check to see if doubling the period would still yield a period less than $q$, if so we double our guessed period. Taking the period, which is now guaranteed to be even, we proceed to calculate $x^{r/2} - 1 \bmod n$ and $x^{r/2} + 1 \bmod n$, calling these values $a$ and $b$, we compute the greatest common denominator of $a$ and $n$, and $b$ and $n$, to see if we have attained a nontrivial factor of $n$. \subsection{Utility Functions for the Simulation} The util.C file contains functions which are called from shor.C. Many of the functions are written to avoid overflow when calculating functions in which there are intermediate steps which may have large values. For example while factoring 15, $x^{255} \bmod 15$ is calculated, but the modular exponentiation function in util.C calculates it in such a manner that $x^{255}$ need not be explicitly calculated. I have commented these helper functions in the code itself. \section{Conclusion} From the dawn of computer science the field has benefited from abstractions which simulate actual computing devices. Between the Turing machine and the Church-Turing Thesis a strong foundation was made for study of the computable and uncomputable. Complexity analysis provides a way to distinguish classes of problems based on their runtime characteristics, and the rough grouping of problems of polynomial runtime as tractable and others as intractable combined with the presumed correctness of the Complexity-Theoretic Church-Turing hypothesis suggest whether or not a problem is tractable is not depended on the model of computation used. Through the principle of superposition in quantum systems we can create useful memory components that are on the scale of an atom or smaller. These quantum memory registers may be able to facilitate exponential computational speed increases in algorithms that can take advantage of quantum parallelism. Peter Shor has shown an algorithm which makes factoring large numbers tractable for a quantum computer, where no such algorithm is published for a classical computer. In doing so has drawn great attention to the field of quantum computing. Due to Shor's algorithm, we may someday have to turn to other means of encrypting data than are today typically employed. L. K. Grover's database search algorithm shows another noteworthy task that a quantum computer can perform faster than any classical computer.(Brassard) The efforts to build a real quantum memory register that functions are in the most preliminary stages. As of the original time of writing this paper, 3-qubit registers had been built. In 2001 Shor's Algorithm was applied to the number 15 at IBM's Almaden Research Center and Stanford University. In 2005 the first qubyte (g. 8-bit quantum register) was created at The Institute of Quantum Optics and Quantum Information at the University of Innsbruck in Austria. In 2009 NIST reads and writes individual qubits, and demonstrates some computing operations on qubits. (Timeline of quantum computing) Operational quantum computers are by no means an inevitable consequence of this research. It may be that the problems surrounding keeping a quantum memory register isolated from any disturbance long enough for a calculation to take place will be insurmountable. In any case, quantum computing will remain an exciting topic for experimentalists and theorists alike for years to come. Hopefully this paper, and the simulation of Shor's algorithm have been as enlightening and fun for the reader as they were for the author. \section{Bibliography} Barenco, Ekert, Sanpera and Machiavello. "Un saut d'echelle pour les calculateurs," \emph{La Recherche}, November 1996. Benioff, p. "The Computer as a Physical System: A Microscopic Quantum Mechanical Hamiltonian Model of Computers as Represented by Turing Machines," \emph{Journal of Statistical Physics}, Vol. 22 (1980), pp. 563-591. Berthiaume, Andre and Brassard, Gilles. "The quantum Challenge to Complexity Theory," \emph{Proceedings of the 7th IEEE Conference on Structure in Complexity Theory} (1992), pp. 132-137. Brassard, Gilles. "Searching a Quantum Phone Book," \emph{Science}, 31 January 1997. Cormen, Thomas H., Leiserson, Charles E., and Rivest, Ronald L. "Introduction to Algorithms," St. Louis: McGraw-Hill, 1994. Deutsch, David and Jozsa, Richard. "Rapid Solution of Problems by Quantum Computation," \emph{Proceedings Royal Society London}, Vol. 439A (1992), pp. 553-558. Feynman, Richard. "Simulating Physics with Computers," \emph{Optics News} Vol. 11 (1982), pp. 467-488. Grover, L. K. "A Fast Quantum Mechanical Algorithm for Database Search," \emph{Proceedings of the 28'th Annual ACM Symposium on the Theory of Computing} (1996), pp. 212-219. Shor, Peter. "Algorithms for Quantum Computation: Discrete Logarithms and Factoring," \emph{Proceedings 35th Annual Symposium on Foundations of Computer Science} (1994), pp. 124-134. Steane, Andrew. "Quantum Computing," \emph{Reports on Progress in Physics}, vol 61 (1998), pp 117-173. Timeline of quantum computing (2015). \emph{Timeline of quantum computing}. Available at: http://en.wikipedia.org/wiki/Timeline\_of\_quantum\_computing [Accessed 14 Jan. 2015]. Williams, Colin P. and Clearwater, Scott H. "Explorations in Quantum Computing," New York: Springer-Verlag, 1998. \section{Glossary} This is a glossary of terms and variables used throughout this paper. \begin{description} \item[$\lambda$]: In the context of Shor's algorithm in integer such that $m = \frac{\lambda q}{r}$. \item[$a$]: A variable that when used in the context of Shor's algorithm is an argument to the function $\mathcal{F}(a) = x^{a} \bmod n$. It may be a single integer, or it may denote a superposition of states. \item[Algorithm]: A method for producing some answer given some inputs. \item[Binary]: The base two number system. For more information see Appendix A. \item[Bit]: A thing which can store one item of information, either a 1 or a 0. \item[Church-Turing Thesis]: The Church-Turing Thesis can be thought of as defining the class of functions which are effectively calculable as those which can be computed by a Turing machine. \item[Classical computer]: A computer whose internal workings behave in manner consistent with classical physics. Data registers in a classical computer can not exist in a superposition of states. By definition (e.g. the Church-Turing Thesis) a classical computer can only compute functions that a Turing machine can compute. \item[Classical physics]: The model that was used to describe physical phenomenon before the advent of quantum physics. The predictions of classical physics with regard to the behavior of fundamental particles are incorrect. \item[Collapse]: A collapse in the context of this paper is what happens to the state vector of a quantum mechanical system when that system is observed or measured. Since the system can only be measured to be in one of its base states, the state vector will collapse from some superposition of base states into the measured state only. \item[Complex number]: A number of the form $a + i * b$, where $a$ and $b$ are real numbers and $i$ is defined to be the square root of negative one. \item[Complex vector space]: A vector space in which the coordinates of a vector are complex numbers. \item[Complexity class]: A grouping of algorithms based on how their memory usage and number of operations scale with the size of the input. \item[Complexity-Theoretic Church-Turing Thesis]: The hypothesis that a Turing machine can efficiently simulate any realistic model of computation. Here efficiently means with polynomial time overhead. It should be noted that this is a hypothesis only, however the advent of a physically realizable computing device which violates this thesis would be a major breakthrough, and there are arguments as to why such a breakthrough may never occur. \item[Coprime]: Integers $a$ and $b$ are coprime if their greatest common denominator is one. \item[Discrete Fourier Transform]: A transformation converts a finite list of equally spaced samples of a function into a list of coefficients of finite combinations of circles, ordered by their frequencies, that have the same values. In Shor's algorithm it is used to calculate and multiple of the inverse period, where the period is the quantity which enables Shor's algorithm to find factors of $a$ number $n$. \item[Exponential]: A function which goes as: $\mathcal{F}(x) = a^{x}$, where $a$ is some constant greater than 1. \item[Exponential time]: This is an attribute of an algorithm which means the number of operations required to compute the answer grows exponentially with the size of the input. \item[gcd]: This is an abbreviation for the mathematical function which calculates the greatest common denominator of two integers. The greatest common denominator of two integers $a$ and $b$ is the largest integer $c$ such that $a / c$ and $b / c$ are integers. \item[Grover's Search Algorithm]: An algorithm designed by L. K. Grover of Bell Labs which finds a element in an unsorted database of size $n$ in $O(\sqrt{n})$ operations on a quantum computer. The lowest known running time for this problem on a classical computer is $O(n)$. \item[Hilbert Space]: A complex linear vector space. The complete state of a $n$ state quantum mechanical system can be represented by a vector in an $n$ dimensional Hilbert Space. \item[$i$]: The square root of -1. \item[Linear vector space]: A vector space such that vectors within the space which are added or multiplied together result in vectors that also lie within the same space. \item[Memory register]: A array of bits on a classical computer. A memory register of size $n$ may store one of $2^{n}$ values. \item[Mutually perpendicular]: In the context of vector spaces mutually perpendicular items are items such that no one can be decomposed into components of the others. \item[$n$]: In the context of Shor's algorithm, a number to be factored. \item[Periodic function]: A function with a period $r$ such that $\mathcal{F}(x) = \mathcal{F}(x + r) = \mathcal{F}(x + 2r)$ and so on. Sine and Cosine are typical examples of periodic functions. \item[Polynomial time]: This is an attribute of an algorithm which means the number of operations required to compute the answer grows polynomially with the size of the input $n$ (E.g. like $n^{c}$ for some constant $c$). \item[$q$]: In the context of Shor's algorithm the power of 2 such that $n^{2} \leq q < 2n^{2}$. \item[Quantum memory register]: A array of $n$ qubits which can exist in any superposition of its $2^{n}$ base states. \item[Quantum parallelism]: The ability of a quantum computer to perform an operation on a quantum memory register which results in the simultaneous calculation of a function function of $2^{n}$ different values where $n$ is the size of the quantum memory register. \item[Quantum physics]: Currently the most complete model for describing the behavior of small physical systems. \item[Qubit]: A two state quantum mechanical system, which can exist in any superposition of the 0 and 1 state. In this paper I have considered a spin-1/2 particle as a possible candidate for a qubit in a physical implementation of a qubit. \item[$r$]: In the context of Shor's algorithm the period of the periodic function $x^{a} \bmod n$. \item[Shor's Algorithm]: A algorithm designed by Peter Shor of Bell Labs which finds factors of a number $n$ in polynomial time relative to the number of bits in $n$'s binary representation on a quantum computer. The fastest published algorithm on a classical computer is slower than polynomial time, and the presumed intractability of this problem is the basis for many cryptographic systems. \item[Spin-1/2 particle]: A particle which can be characterized as having a spin of +1/2 or -1/2. Examples include the proton, neutron, and electron. \item[State vector]: In the context of this paper, the state vector in a Hilbert Space which completely describes a quantum mechanical state vector - such as the state of quantum memory register. \item[Superposition of states]: A mixture of base states. The state vector for a quantum mechanical systems which can be measured in one of $n$ base states can exist as any combination of components of the base states. \item[Turing machine]: A theoretical computing device consisting of an infinite tape divided into cells which can hold a 1, a 0, or a blank and a head which can move around the tape, read and write bits, and change its own internal state. The Church-Turing Thesis defines effectively computable functions as those which can be computed by a Turing machine. \item[Unit vector]: A vector whose length is 1. \item[$x$]: In the context of Shor's algorithm a integer which is coprime to $n$ and used in the function $\mathcal{F}(a) = x^{a} \bmod n$. \end{description} \appendix \section{Mathematics Used in this Paper} This appendix will review some of the mathematics used in the paper. The sections are not intended to be comprehensive for their topic, they only cover what is needed to understand this paper. \subsection{Binary Representation of Numbers} We commonly represent number is base 10, there are 10 elements in our base 10 numbering system, 0,1,2,3,4,5,6,7,8, and 9. In a base $n$ counting system there are $n$ distinct elements, 0 through $n - 1$. When a number which is greater than $n - 1$ needs to be displayed in base $n$ it is represented by a string of the $n - 1$ elements. The value of any given symbol in the string is found by multiplying that symbol by $n^{x}$, where $x$ is the number of symbols in the string that are to the right of the symbol in question. For example, in base 10 the number 982 is equal to $9 * 10^{2} + 8*10^{1} + 2*10^{0}$. Likewise, in base two the number 10101001 is equal to $1*2^{7}+0*2^{6}+1*2^{5}+0*2^{4}+1*2^{3}+0*2^{2}+0*2^{1}+1*2^{0} = 169$ in base 10. \subsection{Complex Numbers} A complex number is a number of the form $a + i*b$, where $a$ and $b$ are real numbers, and $i$ is defined to be the square root of negative one. Addition of two complex numbers $c_{1}$ and $c_{2}$ is defined to be: \[c_{1} = a_{1} + i * b_{1}\] \[c_{2} = a_{2} + i * b_{2}\] \[c_{1} + c_{2} = a_{1} + a_{2} + i * (b_{1} + b_{2})\] The complex conjugate of a complex number $c$, denoted $c^{*}$ is defined to be: \[c = a + i * b\] \[c^{*} = a - i * b\] Multiplication of two complex numbers $c_{1}$ and $c_{2}$ is defined to be: \[c_{1} = a_{1} + i * b_{1}\] \[c_{2} = a_{2} + i * b_{2}\] \[c_{1}*c_{2}=a_{1}*a_{2}-b_{1}*b_{2}+i*(a_{1}*b_{2}+a_{2}*b_{1})\] Euler's Formula for complex numbers states that $e^{ix} = \cos x + i * \sin x$, this relationship is used in the discrete Fourier transform of Shor's algorithm. \subsection{Vector Mathematics} The only operations that are used in Shor's algorithm on vectors are addition, length determination, and scaling. The vector in question is that state vector of a quantum mechanical system, it is a vector in a Hilbert Space. Being a vector in a Hilbert space means that the vector can be represented by projections of the vector onto each of the perpendicular base vectors which define the Hilbert Space. For example, a $n$ state quantum system requires a $n$ dimensional Hilbert Space to represent its state vector. The quantum system can be measured in any of the $n$ states, and to represent this we imagine each of the $n$ states as mutually perpendicular axes within our Hilbert space. Thus the state vector for a system in the $j$'th state is equal to: \[\left( \begin{array}{c} 0 \\ 0 \\ \vdots \\ 1 \\ \vdots \\ 0 \end{array} \right)\] For the $n$ states, where the number at the top of the column is the length of the state vector projected onto the 1st state, and the 1 appears in the $j$'th row. To add two vectors we simply add their components. \[\left( \begin{array}{l} a_{1} \\ a_{2} \\ a_{3} \end{array} \right) + \left( \begin{array}{c} b_{1} \\ b_{2} \\ b_{3} \end{array} \right) = \left( \begin{array}{r} a_{1}+b_{1} \\ a_{2}+b_{2} \\ a_{3}+b_{3} \end{array} \right)\] Since this vector lies in a Hilbert Space the projections of the state vector onto the coordinate axes are allowed to be complex numbers, thus the definition of length is slightly different from what is expected. The length of a vector in a Hilbert space with $n$ components is defined to be: $\sqrt{\sum_{j = 1}^{n}|w_{j}|^{2}}$ where $w_{j}$ is the value of the $j$'th component of the vector, and $|w_{j}|^{2}$ is defined to be $w_{j}$ times its complex conjugate, or when $w_{j} = a + i * b$, $|w_{j}|^{2} = a^{2} + b^{2}$ . To scale a vector by any length $l$ you simply multiply each component of the vector by the value $l$. In particular to scale a vector to length 1 you multiply each component by the inverse length of the vector. \section{Actual Quantum Computers and Further Research} There is great interest in building an actual quantum computer. Many strategies for building a quantum memory register are currently being evaluated. The study of quantum computing is a relatively new field, and the experimental implementation is still in its infancy. Consequently the world wide web is the best place to search for information regarding current research. Note: This section has not been updated since 1999. One possible way to implement a quantum memory register is by using nuclear magnetic resonance techniques to manipulate classical fluids to perform quantum calculations. To read more on this see: http://www.sciam.com/1998/0698issue/0698gershenfeld.html http://squint.stanford.edu/qc/overview.html A second approach is to use ion traps and lasers to hold atomic ions in place. To read more about this method check: http://physics.colorado.edu/faculty/monroe.c.html http://www.bldrdoc.gov/timefreq/ion/index.htm A good place to look for original research papers is the Los Alamos National Laboratories E-print archive at: http://quickreviews.org/ This is by no means a comprehensive list of approaches to and research regarding quantum computing. Given the changing topology of the web and the youth of the subject you will probably be able to search for and find a wealth of information on quantum computing from your favorite web search engine or directory service. \section{Other Quantum Computer Simulators} On a more theoretical note there is some people are writing other, more robust and generic simulators of a quantum computer. The following web sites have information about quantum computing, and have code for simulating the action of a quantum computer. I found these sites to be great resources during my study of quantum computing. Note: This section has not been updated since 1999. http://tph.tuwien.ac.at/~oemer/qc/qcl/qcl.html http://tph.tuwien.ac.at/~oemer/prakt/prakt.html http://www.openqubit.org/ \section{Code for my Simulation of Shor's Algorithm} Here is included the code for my simulation of Shor's algorithm. I have tried to document the code as well as possible, buy if you find something difficult please write me at mjhayward@gmail.com and I will attempt to explain it as well as I can. The code was originally compiled in 1999 with g++ version 2.8.1 on a Pentium Pro 200 running Linux 2.0.33, and had also been compiled it on a UltraSparc running SunOS 5.6 also using g++ 2.8.1, with one modification. Since then, the code has been updated to compile with gcc/g++ 4.8.2 and GNU Make 3.81. I hope that you will find it easy to compile the code with g++. If you have trouble compiling please send me information about your compiler, OS, and error messages, and I will see what I can do. \subsection{complex.h} \begin{verbatim} #ifndef COMPLEX_H #define COMPLEX_H class Complex { public: Complex(); //Default constructor Complex( double, double ); // Argument constructor virtual ~Complex(); //Default destructor. void Set(double new_real, double new_imaginary); //Set data members. double Real() const; //Return the real part. double Imaginary() const; //Return the imaginary part. Complex & operator+(const Complex&); //Overloaded + operator Complex & operator*(const Complex&); //Overloaded * operator Complex & operator=(const Complex&); //Overloaded = operator bool operator==(const Complex&) const; //Overloaded == operator private: double real; double imaginary; }; #endif \end{verbatim} \subsection{complex.C} \begin{verbatim} #include #include "complex.h" //Complex constructor, initializes to 0 + i0. Complex::Complex(): real( 0 ), imaginary( 0 ) {} //Argument constructor. Complex::Complex( double r, double i ): real( r ), imaginary( i ) {} //Complex destructor. Complex::~Complex() {} //Overloaded = operator. Complex & Complex::operator=(const Complex &c) { if (&c != this) { real = c.Real(); imaginary = c.Imaginary(); return *this; } } //Overloaded + operator. Complex & Complex::operator+(const Complex &c) { real += c.Real(); imaginary += c.Imaginary(); return *this; } //Overloaded * operator. Complex & Complex::operator*(const Complex &c) { real = real * c.Real() - imaginary * c.Imaginary(); imaginary = real * c.Imaginary() + imaginary * c.Real(); return *this; } //Overloaded == operator. Small error tolerances. bool Complex::operator==(const Complex &c) const { //This is to take care of round off errors. if (fabs(c.Real() - real) > pow(10,-14)) { return false; } if (fabs(c.Imaginary()- imaginary) > pow(10,-14)) { return false; } return true; } //Sets private data members. void Complex::Set(double new_real, double new_imaginary) { real = new_real; imaginary = new_imaginary; return; } //Returns the real part of the complex number. double Complex::Real() const { return real; } //Returns the imaginary part of the complex number. double Complex::Imaginary() const { return imaginary; } \end{verbatim} \subsection{qureg.h} \begin{verbatim} #include "complex.h" #ifndef QUREG_H #define QUREG_H using namespace std; class QuReg { public: //Default constructor. Size is the size in bits of our register. //In our implementation of Shor's algorithm we will need size bits //to represent our value for "q" which is a number we have chosen //with small prime factors which is between 2n^2 and 3n^2 inclusive //where n is the number we are trying to factor. We envision our //the description of our register of size "S" as 2^S complex //numbers, representing the probability of finding the register on //one of or 2^S base states. Thus we use an array of size 2^S, of //Complex numbers. Thus if the size of our register is 3 bits //array[7] is the probability amplitude of the state |1,1,1>, and //array[7] * Complex Conjugate(array[7]) = probability of choosing //that state. We use normalized state vectors thought the //simulation, thus the sum of all possible states times their //complex conjugates is = 1. QuReg(unsigned long long int size); QuReg(); //Default Constructor QuReg(const QuReg &); //Copy constructor virtual ~QuReg(); //Default destructor. //Measures our quantum register, and returns the decimal //interpretation of the bit string measured. unsigned long long int DecMeasure(); //Dumps all the information about the quantum register. This //function would not be possible for an actual quantum register, it //is only there for debugging. When verbose != 0 we return every //value, when verbose = 0 we return only probability amplitudes //which differ from 0. void Dump(int verbose) const; //Sets state of the qubits using the arrays of complex amplitudes. void SetState(Complex *new_state); //Sets the state to an equal superposition of all possible states //between 0 and number inclusive. void SetAverage(unsigned long long int number); //Normalize the state amplitudes. void Norm(); //Get the probability of a given state. This is used in the //discrete Fourier transformation. In a real quantum computer such //an operation would not be possible, on the flip side, it would //also not be necessary as you could simply build a DFT gate, and //run your superposition through it to get the right answer. Complex GetProb(unsigned long long int state) const; //Return the size of the register. int Size() const; private: unsigned long long int reg_size; Complex *State; }; #endif \end{verbatim} \subsection{qureg.C} \begin{verbatim} #include #include #include #include #include "complex.h" #include "qureg.h" using namespace std; QuReg::QuReg(): reg_size( 0 ), State( NULL ) {} //Constructor. QuReg::QuReg(unsigned long long int size) { reg_size = size; State = new Complex[(unsigned long long int)pow(2, reg_size)]; srand(time(NULL)); } //Copy Constructor QuReg::QuReg(const QuReg & old) { reg_size = old.reg_size; unsigned long long int reg_length = (unsigned long long int) pow(2, reg_size); State = new Complex[reg_length]; for (unsigned int i = 0 ; i < reg_length ; i++) { State[i] = old.State[i]; } } //Destructor. QuReg::~QuReg() { if ( State ) { delete [] State; } } //Return the probability amplitude of the state'th state. Complex QuReg::GetProb(unsigned long long int state) const { if (state >= pow(2, reg_size)) { cout << "Invalid state index " << state << " requested!" << endl << flush; throw -1; } else { return(State[state]); } } //Normalize the probability amplitude, this ensures that the sum of //the sum of the squares of all the real and imaginary components is //equal to one. void QuReg::Norm() { double b = 0; double f, g; for (unsigned long long int i = 0; i < pow(2, reg_size) ; i++) { b += pow(State[i].Real(), 2) + pow(State[i].Imaginary(), 2); } b = pow(b, -.5); for (unsigned long long int i = 0; i < pow(2, reg_size) ; i++) { f = State[i].Real() * b; g = State[i].Imaginary() * b; State[i].Set(f, g); } } //Returns the size of the register. int QuReg::Size() const { return reg_size; } //Measure a state, and return the decimal value measured. Collapse //the state so that the probability of measuring the measured value in //the future is 1, and the probability of measuring any other state is //0. unsigned long long int QuReg::DecMeasure() { int done = 0; int DecVal = -1; //-1 is an error, we did not measure anything. double rand1, a, b; rand1 = rand()/(double)RAND_MAX; a = b = 0; for (unsigned long long int i = 0 ; i < pow(2, reg_size) ;i++) { if (!done ){ b += pow(State[i].Real(), 2) + pow(State[i].Imaginary(), 2); if (b > rand1 && rand1 > a) { //We have just measured the i state. for (unsigned long long int j = 0; j < pow(2, reg_size) ; j++) { State[j].Set(0,0); } State[i].Set(1,0); DecVal = i; done = 1; } a += pow(State[i].Real(), 2) + pow(State[i].Imaginary(), 2); } } return DecVal; } //For debugging, output information about the register. void QuReg::Dump(int verbose) const { for (unsigned long long int i = 0 ; i < pow(2, reg_size) ; i++) { if (verbose || fabs(State[i].Real()) > pow(10,-14) || fabs(State[i].Imaginary()) > pow(10,-14)) { cout << "State " << i << " has probability amplitude " << State[i].Real() << " + i" << State[i].Imaginary() << endl << flush; } } } //Set the states to those given in the new_state array. void QuReg::SetState(Complex *new_state) { //Beware, this function will cause segfaults if new_state is too //small. for (unsigned long long int i = 0 ; i < pow(2, reg_size) ; i++) { State[i].Set(new_state[i].Real(), new_state[i].Imaginary()); } } //Set the State to an equal superposition of the integers 0 -> number //- 1 void QuReg::SetAverage(unsigned long long int number) { if (number >= pow(2, reg_size)) { cout << "Error, initializing past end of array in qureg::SetAverage.\n"; } else { double prob = pow(number, -.5); for (unsigned long long int i = 0 ; i <= number ; i++) { State[i].Set(prob, 0); } } } \end{verbatim} \subsection{shor.C} \begin{verbatim} #include #include #include #include #include "complex.h" #include "qureg.h" #include "util.h" using namespace std; int main() { //Establish a random seed. srand(time(NULL)); //Output standard greeting. cout << "Welcome to the simulation of Shor's algorithm." << endl << "There are four restrictions for Shor's algorithm:" << endl << "1) The number to be factored (n) must be >= 15." << endl << "2) The number to be factored must be odd." << endl << "3) The number must not be prime." << endl << "4) The number must not be a prime power." << endl << endl << "There are efficient classical methods of factoring " << "any of the above numbers, or determining that they are prime." << endl << endl << "Input the number you wish to factor." << endl << flush; //n is the number we are going to factor, get n. unsigned long long int n; cin >> n; cout << "Step 1 starting." << endl << flush; //Test to see if n is factorable by Shor's algorithm. //Exit if the number is even. if (n%2 == 0) { cout << "\tError, the number must be odd!" << endl << flush; exit(0); } //Exit if the number is prime. if (TestPrime(n)) { cout << "\tError, the number must not be prime!" << endl << flush; exit(0); } //Prime powers are prime numbers raised to integral powers. //Exit if the number is a prime power. if (TestPrimePower(n)) { cout << "\tError, the number must not be a prime power!" << endl << flush; exit(0); } cout << "Step 1 complete." << endl << flush; cout << "Step 2 starting." << endl << flush; //Now we must figure out how big a quantum register we need for our //input, n. We must establish a quantum register big enough to hold //an equal superposition of all integers 0 through q - 1 where q is //the power of two such that n^2 <= q < 2n^2. unsigned long long int q; cout << "\tSearching for q, the smallest power of 2 greater than or equal to n^2." << endl << flush; q = GetQ(n); cout << "\tFound q to be " << q << "." << endl << flush; cout << "Step 2 complete." << endl << flush; cout << "Step 3 starting." << endl << flush; //Now we must pick a random integer x, coprime to n. Numbers are //coprime when their greatest common denominator is one. One is not //a useful number for the algorithm. unsigned long long int x = 0; cout << "\tSearching for x, a random integer coprime to n." << endl << flush; x = 1+ (unsigned long long int)((n-1)*(double)rand()/(double)RAND_MAX); while (GCD(n,x) != 1 || x == 1) { x = 1 + (unsigned long long int)((n-1)*(double)rand()/(double)RAND_MAX); } cout << "\tFound x to be " << x << "." << endl << flush; cout << "Step 3 complete." << endl << flush; //Create the register. cout << "Step 4 starting." << endl << flush; QuReg * reg1 = new QuReg(RegSize(q) - 1); cout << "\tMade register 1 with register size = " << RegSize(q) << endl << flush; //This array will remember what values of q produced for x^q mod n. //It is necessary to retain these values for use when we collapse //register one after measuring register two. In a real quantum //computer these registers would be entangled, and thus this extra //bookkeeping would not be needed at all. The laws of quantum //mechanics dictate that register one would collapse as well, and //into a state consistent with the measured value in resister two. unsigned long long int * modex = new unsigned long long int[q]; //This array holds the probability amplitudes of the collapsed state //of register one, after register two has been measured it is used //to put register one in a state consistent with that measured in //register two. Complex * collapse = new Complex[q]; //This is a temporary value. Complex tmp; //This is a new array of probability amplitudes for our second //quantum register, that populated by the results of x^a mod n. Complex * mdx = new Complex[(unsigned long long int)pow(2,RegSize(n))]; // This is the second register. It needs to be big enough to hold // the superposition of numbers ranging from 0 -> n - 1. QuReg *reg2 = new QuReg(RegSize(n)); cout << "\tCreated register 2 of size " << RegSize(n) << endl << flush; cout << "Step 4 complete." << endl << flush; //This is a temporary value. unsigned long long int tmpval; //This is a temporary value. unsigned long long int value; //c is some multiple lambda of q/r, where q is q in this program, //and r is the period we are trying to find to factor n. m is the //value we measure from register one after the Fourier //transformation. double c, m; //This is used to store the denominator of the fraction p / den where //p / den is the best approximation to c with den <= q. unsigned long long int den; //This is used to store the numerator of the fraction p / den where //p / den is the best approximation to c with den <= q. unsigned long long int p; //The integers e, a, and b are used in the end of the program when //we attempts to calculate the factors of n given the period it //measured. //Factor is the factor that we find. unsigned long long int e, a, b, factor; //Shor's algorithm can sometimes fail, in which case you do it //again. The done variable is set to 0 when the algorithm has //failed. Only try a maximum number of tries. unsigned int done = 0; unsigned int tries = 0; while (!done) { if (tries >= 5) { cout << "\tThere have been five failures, giving up." << endl << flush; exit(0); } cout << "Step 5 starting attempt: " << tries+1 << endl << flush; //Now populate register one in an even superposition of the //integers 0 -> q - 1. reg1->SetAverage(q - 1); cout << "Step 5 complete." << endl << flush; cout << "Step 6 starting attempt: " << tries+1 << endl << flush; //Now we preform a modular exponentiation on the superposed //elements of reg 1. That is, perform x^a mod n, but exploiting //quantum parallelism a quantum computer could do this in one //step, whereas we must calculate it once for each possible //measurable value in register one. We store the result in a new //register, reg2, which is entangled with the first register. //This means that when one is measured, and collapses into a base //state, the other register must collapse into a superposition of //states consistent with the measured value in the other.. The //size of the result modular exponentiation will be at most n, so //the number of bits we will need is therefore less than or equal //to log2 of n. At this point we also maintain a array of what //each state produced when modularly exponised, this is because //these registers would actually be entangled in a real quantum //computer, this information is needed when collapsing the first //register later. //This counter variable is used to increase our probability amplitude. tmp.Set(1,0); //This for loop ranges over q, and puts the value of x^a mod n in //modex[a]. It also increases the probability amplitude of the value //of mdx[x^a mod n] in our array of complex probabilities. for (unsigned long long int i = 0 ; i < q ; i++) { //We must use this version of modexp instead of c++ builtins as //they overflow when x^i is large. tmpval = modexp(x,i,n); modex[i] = tmpval; mdx[tmpval] = mdx[tmpval] + tmp; } //Set the state of register two to what we calculated it should be. reg2->SetState(mdx); //Normalize register two, so that the probability of measuring a //state is given by summing the squares of its probability //amplitude. reg2->Norm(); cout << "Step 6 complete." << endl << flush; cout << "Step 7 starting attempt: " << tries+1 << endl << flush; //Now we measure reg2. value = reg2->DecMeasure(); //Now we must using the information in the array modex collapse //the state of register one into a state consistent with the value //we measured in register two. for (unsigned long long int i = 0 ; i < q ; i++) { if (modex[i] == value) { collapse[i].Set(1,0); } else { collapse[i].Set(0,0); } } //Now we set the state of register one to be consistent with what //we measured in state two, and normalize the probability //amplitudes. reg1->SetState(collapse); reg1->Norm(); cout << "Step 7 complete." << endl << flush; //Here we do our Fourier transformation. cout << "Step 8 starting attempt: " << tries+1 << endl << flush; DFT(reg1, q); cout << "Step 8 complete." << endl << flush; cout << "Step 9 starting attempt: " << tries+1 << endl << flush; //Next we measure register one, due to the Fourier transform the //number we measure, m will be some multiple of lambda/r, where //lambda is an integer and r is the desired period. m = reg1->DecMeasure(); cout << "\tValue of m measured as: " << m << endl << flush; cout << "Step 9 complete." << endl << flush; //If nothing goes wrong from here on out we are done. done = 1; //If we measured zero, we have gained no new information about the //period, we must try again. if (m == 0) { cout << "\tMeasured, 0 this trial a failure!" << endl << flush; done = 0; } //The DecMeasure subroutine will return -1 as an error code, due //to rounding errors it will occasionally fail to measure a state. if (m == -1) { cout << "\tWe failed to measure anything, this trial a failure!" << endl << flush; done = 0; } //If nothing has gone wrong, try to determine the period of our //function, and get factors of n. if (done) { //Now c =~ lambda / r for some integer lambda. Borrowed with //modifications from Berhnard Ohpner. c = (double)m / (double)q; cout << "Steps 10 and 11 starting attempt: " << tries+1 << endl << flush; //Calculate the denominator of the best rational approximation //to c with den < q. Since c is lambda / r for some integer //lambda, this will provide us with our guess for r, our period. den = denominator(c, q); //Calculate the numerator from the denominator. p = (unsigned long long int)floor(den * c + 0.5); //Give user information. cout << "\tMeasured m: " << m << ", rational approximation for m/q=" << c << " is: " << p << " / " << den << endl << flush; //The denominator is our period, and an odd period is not //useful as a result of Shor's algorithm. If the denominator //times two is still less than q we can use that. if (den % 2 == 1 && 2 * den < q ){ cout << "\tOdd candidate for r found, expanding by 2\n"; p = 2 * p; den = 2 * den; } //Initialize helper variables. e = a = b = factor = 0; // Failed if odd denominator. if (den % 2 == 1) { cout << "\tOdd period found. This trial failed." << " \t Trying again." << endl << flush; done = 0; } else { //Calculate candidates for possible common factors with n. cout << "\tCandidate period is " << den << endl << flush; e = modexp(x, den / 2, n); a = (e + 1) % n; b = (e + n - 1) % n; cout << "\t" << x << "^" << den / 2 << " + 1 mod " << n << " = " << a << "," << endl << "\t" << x << "^" << den / 2 << " - 1 mod " << n << " = " << b << endl << flush; factor = max(GCD(n,a),GCD(n,b)); } } //GCD will return a -1 if it tried to calculate the GCD of two //numbers where at some point it tries to take the modulus of a //number and 0. if (factor == -1) { cout << "\tError, tried to calculate n mod 0 for some n. Trying again." << endl << flush; done = 0; } if ((factor == n || factor == 1) && done == 1) { cout << "\tFound trivial factors 1 and " << n << ". Trying again." << endl << flush; done = 0; } //If nothing else has gone wrong, and we got a factor we are //finished. Otherwise start over. if (factor != 0 && done == 1) { cout << "\t" << n << " = " << factor << " * " << n / factor << endl << flush; } else if (done == 1) { cout << "\tFound factor to be 0, error. Trying again." << endl << flush; done = 0; } cout << "Steps 10 and 11 complete." << endl << flush; tries++; } delete reg1; delete reg2; delete [] modex; delete [] collapse; delete [] mdx; return 1; } \end{verbatim} \subsection{util.h} \begin{verbatim} #include "qureg.h" #ifndef UTIL_H #define UTIL_H //This function takes an integer input and returns 1 if it is a prime //number, and 0 otherwise. unsigned long long int TestPrime(unsigned long long int n); //This function takes an integer input and returns 1 if it is equal to //a prime number raised to an integer power, and 0 otherwise. unsigned long long int TestPrimePower(unsigned long long int n); //This function computes the greatest common denominator of two integers. //Since the modulus of a number mod 0 is not defined, we return a -1 as //an error code if we ever would try to take the modulus of something and //zero. unsigned long long int GCD(unsigned long long int a, unsigned long long int b); //This function takes and integer argument, and returns the size in bits //needed to represent that integer. unsigned long long int RegSize(unsigned long long int a); //q is the power of two such that n^2 <= q < 2n^2. unsigned long long int GetQ(unsigned long long int n); //This function takes three integers, x, a, and n, and returns x^a mod //n. This algorithm is known as the "Russian peasant method," I //believe, and avoids overflow by never calculating x^a directly. unsigned long long int modexp(unsigned long long int x, unsigned long long int a, unsigned long long int n); // This function finds the denominator q of the best rational // denominator q for approximating p / q for c with q < qmax. unsigned long long int denominator(double c, unsigned long long int qmax); //This function takes two integer arguments and returns the greater of //the two. unsigned long long int max(unsigned long long int a, unsigned long long int b); //This function computes the discrete Fourier transformation on a register's //0 -> q - 1 entries. void DFT(QuReg * reg, unsigned long long int q); #endif \end{verbatim} \subsection{util.C} \begin{verbatim} #include #include #include "qureg.h" using namespace std; //This function takes an integer input and returns 1 if it is a prime //number, and 0 otherwise. // // Not optimized at all. unsigned long long int TestPrime(unsigned long long int n) { unsigned long long int i; for (i = 2 ; i <= floor(sqrt(n)) ; i++) { if (n % i == 0) { return(0); } } return(1); } //This function takes an integer input and returns 1 if it is equal to //a prime number raised to an integer power, and 0 otherwise. unsigned long long int TestPrimePower(unsigned long long int n) { unsigned long long int i,j; j = 0; i = 2; while ((i <= floor(pow(n, .5))) && (j == 0)) { if((n % i) == 0) { j = i; } i++; } for (unsigned long long int i = 2 ; i <= (floor(log(n) / log(j)) + 1) ; i++) { if(pow(j , i) == n) { return(1); } } return(0); } //This function computes the greatest common denominator of two integers. //Since the modulus of a number mod 0 is not defined, we return a -1 as //an error code if we ever would try to take the modulus of something and //zero. unsigned long long int GCD(unsigned long long int a, unsigned long long int b) { unsigned long long int d; if (b != 0) { while (a % b != 0) { d = a % b; a = b; b = d; } } else { return -1; } return(b); } //This function takes and integer argument, and returns the size in bits //needed to represent that integer. unsigned long long int RegSize(unsigned long long int a) { unsigned long long int size = 0; while(a != 0) { a = a>>1; size++; } return(size); } //q is the power of two such that n^2 <= q < 2n^2. unsigned long long int GetQ(unsigned long long int n) { unsigned long long int power = 8; //256 is the smallest q ever is. while (pow(2,power) < pow(n,2)) { power = power + 1; } return 1 << power; } //This function takes three integers, x, a, and n, and returns x^a mod //n. This algorithm is known as the "Russian peasant method," I //believe, and avoids overflow by never calculating x^a directly. unsigned long long int modexp(unsigned long long int x, unsigned long long int a, unsigned long long int n) { unsigned long long int value = 1; unsigned long long int tmp; tmp = x % n; while (a > 0) { if (a & 1) { value = (value * tmp) % n; } tmp = tmp * tmp % n; a = a>>1; } return value; } // This function finds the denominator q of the best rational // denominator q for approximating p / q for c with q < qmax. unsigned long long int denominator(double c, unsigned long long int qmax) { double y = c; double z; unsigned long long int q0 = 0; unsigned long long int q1 = 1; unsigned long long int q2 = 0; while (1) { z = y - floor(y); if (z < 0.5 / pow(qmax,2)) { return(q1); } if (z != 0) { //Can't divide by 0. y = 1 / z; } else { //Warning this is broken if q1 == 0, but that should never happen. return(q1); } q2 = (unsigned long long int)floor(y) * q1 + q0; if (q2 >= qmax) { return(q1); } q0 = q1; q1 = q2; } } //This function takes two integer arguments and returns the greater of //the two. unsigned long long int max(unsigned long long int a, unsigned long long int b) { if (a > b) { return(a); } return(b); } //This function computes the discrete Fourier transformation on a register's //0 -> q - 1 entries. void DFT(QuReg * reg, unsigned long long int q) { //The Fourier transform maps functions in the time domain to //functions in the frequency domain. Frequency is 1/period, thus //this Fourier transform will take our periodic register, and peak it //at multiples of the inverse period. Our Fourier transformation on //the state a takes it to the state: q^(-.5) * Sum[c = 0 -> c = q - 1, //c * e^(2*Pi*i*a*c / q)]. Remember, e^ix = cos x + i*sin x. Complex * init = new Complex[q]; Complex tmpcomp( 0, 0 ); //Here we do things that a real quantum computer couldn't do, such //as look as individual values without collapsing state. The good //news is that in a real quantum computer you could build a gate //which would what this out all in one step. unsigned long long int count = 0; double tmpreal = 0; double tmpimag = 0; double tmpprob = 0; for (unsigned long long int a = 0 ; a < q ; a++) { //This if statement helps prevent previous round off errors from //propagating further. if ((pow(reg->GetProb(a).Real(),2) + pow(reg->GetProb(a).Imaginary(),2)) > pow(10,-14)) { for (unsigned long long int c = 0 ; c < q ; c++) { tmpcomp.Set(pow(q,-.5) * cos(2*M_PI*a*c/q), pow(q,-.5) * sin(2*M_PI*a*c/q)); init[c] = init[c] + (reg->GetProb(a) * tmpcomp); } } count++; if (count == 100) { cout << "Making progress in Fourier transform, " << 100*((double)a / (double)(q - 1)) << "% done!" << endl << flush; count = 0; } } reg->SetState(init); reg->Norm(); delete [] init; } \end{verbatim} \subsection{qubit.C} Below is the much simpler code for a single qubit. You may find it entertaining and educational to write a simple driver program for this class to toy with a single qubit. \begin{verbatim} #include #include #include #include "complex.h" using namespace std; class Qubit { public: Qubit(); //Default constructor. virtual ~Qubit(); //Default destructor. int Measure(); //Returns zero_state = 0 or one_state = 1 in accordance //with the probabilities of zero_state and one_state. void Dump() const; //Prints our zero_state, and one_state, without //disturbing anything, this operation has no physical //realization, it is only for information and debugging. //It should never be used in an algorithm for //information. void SetState(const Complex& zero_prob, const Complex& one_prob); // Takes two complex numbers and sets the states to //those values. void SetAverage(); //Sets the state to 2^(1/2) zero_state, 2^(1/2) //one_state. No imaginary/phase component. double MCC(int state) const; //Multiply the zero or one state by its complex //conjugate, and return the value. This value //will always be a real number, with no imaginary //component. private: Complex zero_state; //The probability of finding the Qubit in the zero or all imaginary //state. I currently use only the real portion. Complex one_state; //The probability of finding the Qubit in the one or all real state. //I currently use only the real portion. //|zero_state|^2 + |one_state|^2 should always be 1. //This notation means z_s * ComplexConjugate(z_s) + o_s * //ComplexConjugate(o_s) = 1. }; //Qubit constructor, initially sets things in the zero state. Qubit::Qubit() { zero_state = Complex(1,0); one_state = Complex(0,0); srand(time(NULL)); } //Returns _state * ComplexConjugate(_state). double Qubit::MCC(int state) const { if (state == 0) { return (pow(zero_state.Real(), 2) + pow(zero_state.Imaginary(), 2)); } else { return (pow(one_state.Real(), 2) + pow(one_state.Imaginary(), 2)); } } //Measurement operator. Destructively collapses superpositions. int Qubit::Measure() { double rand1 = rand()/(double)RAND_MAX; //This assumes that the sum of the squares of the amplitudes are = 1 if (MCC(0) > rand1) { zero_state.Set(1,0); one_state.Set(0,0); return 0; } else { zero_state.Set(0,0); one_state.Set(1,0); return 1; } } //Outputs state info for our qubit. For debugging purposes. void Qubit::Dump() const{ cout << "Amplitude of the zero state is " << zero_state.Real() << " +i" << zero_state.Imaginary() << endl << flush; cout << "Amplitude of the one state is " << one_state.Real() << " +i" << one_state.Imaginary() << endl << flush; } //Sets the zero and one states to arbitrary amplitudes. Outputs //an error message if the two values MCC'ed != 1 + 0i. void Qubit::SetState(const Complex& zero_prob, const Complex& one_prob) { zero_state = zero_prob; one_state = one_prob; //Determine if |zero_state|^2 + |one_state|^2 == 1, if not we //are not in a valid state. double probab; probab = MCC(0) + MCC(1); if (fabs(probab - 1) > pow(10, -10)) { //This funny expression //allows us some rounding errors. cout << "Warning, total probability for in SetState is different from 1." << endl << flush; } } //Sets the qubit 1/2 way between the 0 state and the 1 state. No phase. void Qubit::SetAverage() { zero_state.Set(pow(2, -.5), 0); one_state.Set(pow(2, -.5), 0); } \end{verbatim} \pagebreak \section{Sample Output} \subsection{Sample Output for $n = 17$} \begin{verbatim} Welcome to the simulation of Shor's algorithm. There are four restrictions for Shor's algorithm: 1) The number to be factored must be >= 15. 2) The number to be factored must be odd. 3) The number must not be prime. 4) The number must not be a prime power. There are efficient classical methods of factoring any of the above numbers, or determining that they are prime. Input the number you wish to factor. 17 Step 1 starting. Error, the number must not be prime! \end{verbatim} \pagebreak \subsection{Sample Output for $n = 15$} \begin{verbatim} Welcome to the simulation of Shor's algorithm. There are four restrictions for Shor's algorithm: 1) The number to be factored (n) must be >= 15. 2) The number to be factored must be odd. 3) The number must not be prime. 4) The number must not be a prime power. There are efficient classical methods of factoring any of the above numbers, or determining that they are prime. Input the number you wish to factor. 15 Step 1 starting. Step 1 complete. Step 2 starting. Searching for q, the smallest power of 2 greater than or equal to n^2. Found q to be 256. Step 2 complete. Step 3 starting. Searching for x, a random integer coprime to n. Found x to be 13. Step 3 complete. Step 4 starting. Made register 1 with register size = 9 Created register 2 of size 4 Step 4 complete. Step 5 starting attempt: 1 Step 5 complete. Step 6 starting attempt: 1 Step 6 complete. Step 7 starting attempt: 1 Step 7 complete. Step 8 starting attempt: 1 Making progress in Fourier transform, 38.8235% done! Making progress in Fourier transform, 78.0392% done! Step 8 complete. Step 9 starting attempt: 1 Value of m measured as: 0 Step 9 complete. Measured, 0 this trial a failure! Steps 10 and 11 complete. Step 5 starting attempt: 2 Step 5 complete. Step 6 starting attempt: 2 Step 6 complete. Step 7 starting attempt: 2 Step 7 complete. Step 8 starting attempt: 2 Making progress in Fourier transform, 38.8235% done! Making progress in Fourier transform, 78.0392% done! Step 8 complete. Step 9 starting attempt: 2 Value of m measured as: 0 Step 9 complete. Measured, 0 this trial a failure! Steps 10 and 11 complete. Step 5 starting attempt: 3 Step 5 complete. Step 6 starting attempt: 3 Step 6 complete. Step 7 starting attempt: 3 Step 7 complete. Step 8 starting attempt: 3 Making progress in Fourier transform, 38.8235% done! Making progress in Fourier transform, 78.0392% done! Step 8 complete. Step 9 starting attempt: 3 Value of m measured as: 128 Step 9 complete. Steps 10 and 11 starting attempt: 3 Measured m: 128, rational approximation for m/q=0.5 is: 1 / 2 Candidate period is 2 13^1 + 1 mod 15 = 14, 13^1 - 1 mod 15 = 12 15 = 3 * 5 Steps 10 and 11 complete. \end{verbatim} \pagebreak \subsection{Sample Output for $n = 33$} \begin{verbatim} Welcome to the simulation of Shor's algorithm. There are four restrictions for Shor's algorithm: 1) The number to be factored (n) must be >= 15. 2) The number to be factored must be odd. 3) The number must not be prime. 4) The number must not be a prime power. There are efficient classical methods of factoring any of the above numbers, or determining that they are prime. Input the number you wish to factor. 33 Step 1 starting. Step 1 complete. Step 2 starting. Searching for q, the smallest power of 2 greater than or equal to n^2. Found q to be 2048. Step 2 complete. Step 3 starting. Searching for x, a random integer coprime to n. Found x to be 8. Step 3 complete. Step 4 starting. Made register 1 with register size = 12 Created register 2 of size 6 Step 4 complete. Step 5 starting attempt: 1 Step 5 complete. Step 6 starting attempt: 1 Step 6 complete. Step 7 starting attempt: 1 Step 7 complete. Step 8 starting attempt: 1 Making progress in Fourier transform, 4.83635% done! Making progress in Fourier transform, 9.72154% done! Making progress in Fourier transform, 14.6067% done! Making progress in Fourier transform, 19.4919% done! Making progress in Fourier transform, 24.3771% done! Making progress in Fourier transform, 29.2623% done! Making progress in Fourier transform, 34.1475% done! Making progress in Fourier transform, 39.0327% done! Making progress in Fourier transform, 43.9179% done! Making progress in Fourier transform, 48.8031% done! Making progress in Fourier transform, 53.6883% done! Making progress in Fourier transform, 58.5735% done! Making progress in Fourier transform, 63.4587% done! Making progress in Fourier transform, 68.3439% done! Making progress in Fourier transform, 73.2291% done! Making progress in Fourier transform, 78.1143% done! Making progress in Fourier transform, 82.9995% done! Making progress in Fourier transform, 87.8847% done! Making progress in Fourier transform, 92.7699% done! Making progress in Fourier transform, 97.6551% done! Step 8 complete. Step 9 starting attempt: 1 Value of m measured as: 409 Step 9 complete. Steps 10 and 11 starting attempt: 1 Measured m: 409, rational approximation for m/q=0.199707 is: 273 / 1367 Odd period found. This trial failed. Trying again. Steps 10 and 11 complete. Step 5 starting attempt: 2 Step 5 complete. Step 6 starting attempt: 2 Step 6 complete. Step 7 starting attempt: 2 Step 7 complete. Step 8 starting attempt: 2 Making progress in Fourier transform, 4.83635% done! Making progress in Fourier transform, 9.72154% done! Making progress in Fourier transform, 14.6067% done! Making progress in Fourier transform, 19.4919% done! Making progress in Fourier transform, 24.3771% done! Making progress in Fourier transform, 29.2623% done! Making progress in Fourier transform, 34.1475% done! Making progress in Fourier transform, 39.0327% done! Making progress in Fourier transform, 43.9179% done! Making progress in Fourier transform, 48.8031% done! Making progress in Fourier transform, 53.6883% done! Making progress in Fourier transform, 58.5735% done! Making progress in Fourier transform, 63.4587% done! Making progress in Fourier transform, 68.3439% done! Making progress in Fourier transform, 73.2291% done! Making progress in Fourier transform, 78.1143% done! Making progress in Fourier transform, 82.9995% done! Making progress in Fourier transform, 87.8847% done! Making progress in Fourier transform, 92.7699% done! Making progress in Fourier transform, 97.6551% done! Step 8 complete. Step 9 starting attempt: 2 Value of m measured as: 1438 Step 9 complete. Steps 10 and 11 starting attempt: 2 Measured m: 1438, rational approximation for m/q=0.702148 is: 719 / 1024 Candidate period is 1024 8^512 + 1 mod 33 = 32, 8^512 - 1 mod 33 = 30 33 = 3 * 11 Steps 10 and 11 complete. \end{verbatim} \end{document}