Implementing a Simulation of Shor S Quantum Factoring Algorithm

Implementing A Simulation of Shor’s Quantum Factoring Algorithm

Introduction

In 1994, AT&T researcher Peter W. Shor presented two algorithms that used quantum mechanics to solve problems regarded as being exponentially hard on conventional computers (1): factoring integers and finding discrete logarithms. His quantum algorithms accomplish the task in merely polynomial time, and that threatens the security of cryptographic systems that depend upon the problems being exponentially difficult. The loss of exponential complexity, aside from the cryptographic applications, makes the algorithms intensely interesting and adjures us to obtain a thorough understanding of how and where they accomplish their reduction of time complexity. This paper is a narrative of its author’s attempt to understand the quantum algorithm for factoring by implementing a simulation of it on a classical computer.

The Algorithm

First of all, a compact summary of the factoring algorithm should be presented. The core of the factoring algorithm (2), the part that uses quantum circuits, is a quantum order finding algorithm (3). That algorithm in turn is built around quantum phase estimation (4), which is entirely quantum mechanical. It consists of a pair of multi-qubit registers, the first of which is set to all zero (false qubits), and the second of which is put into some other special prepared state. The qubits of the first register are then independently modified by Hadamard gates (5) (which changes the qubits from a pure state to a mixed state) and then the entire assemblage of qubits is run through a set of quantum circuits that implements some function of interest (a “black box” operation). This is where the search process becomes de-exponentiated. The qubits form an entangled state due to the black box having gates that use one or more of the qubits at time to control what is done to another qubit. That ends up forming an exponentially large state in terms of the number of basis elements in it. The entanglement free-expands to an exponentially large size with only a polynomial quantity of operations. In contrast, a classical process must use an exponential number of operations to search among the exponential quantity of states. Note that the black box must have qubit-controlled qubit gates to grow any entanglement. What the entanglement state’s expansion ends up searching for depends upon what the arrangement of circuits is inside the black box.

After the black box operation does its search work, the second register’s qubits are measured. This will choose one of its sets of states out of all the states of the entanglement. The equivalent new entanglement may be regarded as having shrunk by some orders of magnitude. The reduced entanglement is then subjected to a quantum inverse Fourier transform. That creates a new entanglement of the same size as the untransformed one, but hopefully with an exponentially sparser set of states in it, each having exponentially higher probabilities of being the outcome when measured. The first register’s qubits are finally measured, collapsing this final entangled state to one outcome, usable as the output result of the quantum circuit.

The order finding algorithm adapts the phase estimation circuit described above by modifying two things: the second register is prepared into a state of all true qubits (notated as |1>), and by having the black box circuits perform a function of xkmod N, where k is the first register of qubits, and the value of the function ends up in the second (“work”) register (6). The order finding algorithm then tacks on a classical numerical computation after running the quantum circuits by doing a continued fractions expansion upon the measurement obtained from the first register (7).

Finally getting to the quantum factoring algorithm, the usage of the quantum circuit is preceded by some classical low-cost numerical tests for prime factors. If they don’t find anything then the quantum circuit is employed. For the quantum circuit, the second (“work”) register is prepared into a state of all false qubits (|0>) instead of true ones, and the rest of the quantum circuit is the same as for order finding. The continued fractions algorithm is also in place after the quantum circuit, followed by another classical numerical test of a pair of greatest common divisors for the number being factored to find one or two non-trivial factors.

The Simulation

Implementing a simulation of a quantum system involves implementing a system of linear algebra. The data objects of a very limited symbolic algebra system have been coded from scratch in C++ and undergone several revisions. The end product is a collection of six notable classes:

· A class for representing a complex number (the elements of quantum state vectors). The class has particular features for easy recognition:

o Separate integer numerators and denominators to represent the squared magnitude of the complex number.

o Separate integer numerators and denominators to represent a fraction that 2πi is multiplied by for the angle of the complex number, expressed in exponential form.

o The two features above are not necessary for numerical purposes, but are there for the “symbolic algebra system” aspect of making the numbers more recognizable by human users. Stating that a number is square root of an integer divided by another integer is more recognizable and preferable to displaying a cryptic floating point number.

o Both polar and rectangular floating point members are maintained in parallel with the discrete members for handling cases such as adding numbers with different phases that are not diametrically opposed.

· A class for representing a single qubit, comprising a two element state vector.

· A class for an arbitrarily long quantum state vector used to represent entangled states.

· A class for a Hadamard gate for the single qubit class.

· A class for a controlled rotation gate for the single qubit class (returning a single qubit turns out to be an incorrect implementation, a four-element long entanglement must be returned instead).

· A class for a quantum Fourier transform (and its inverse).

Each of the data classes (complex number, qubit, long quantum state) comes equipped with the ability to print itself out neatly, multiply by each other, and add to each other.

The effort for this project focused on just the quantum circuit portion at the core of factoring, neglecting the numerical methods fore and aft of the quantum circuit. The major guidance for setting milestones and testing the project was Box 5.4 on page 235 (the entire page) of Nielsen & Chuang, the subject of the next section.

Simulating The Factoring of 15

The first action for setting up a simulation of the factoring according to Box 5.4 is to establish the number of qubits for the first register (used for exponentiation) and the number of qubits for the second register (work). The number of work qubits is set to 4 by the example, so the number of qubits for the exponent register is set to 11 in accordance with t=2L+1+log⁡(2+12ϵ (8), where L is the number of work qubits (4), and the value for ϵ is 14 for the maximum error (9).

Setting the exponent and work qubits to all |0>s establishes the state |0>|0> as recommended in box 5.4, and the next step of Hadamard gating all the exponent qubits to 12(|0>+ |1>) each is further taken to create the state 12tk=02t-1|k>|0>. In neither case, though, are the registers of qubits actually tensored with each other as the notation indicates. That step is saved for creating the result of sending the qubits through the “black box” of xkmod N.

At this point the parameterization of the black box must be mentioned for the sake of determining what value for x is exponentiated and what value for N is divided by for the number to factor. Rather obviously, N is a parameter given to the entire algorithm at the beginning, and stated to be 15. The value for x is generated by step 3 of the algorithm (page 233), the last classical numerical step before invoking the quantum process. Step 3 has to fail for the x it generated to be employed with the quantum process. Box 5.4 uses 7 for x as a number that does not satisfy step 3.

The proper way to simulate the next step is to perform all the operations of the quantum circuits of the black box upon all of the involved qubits. Such a complex circuit was not already given in the text, and the time available to try constructing and testing cycles was considered to be too short against the need to develop the operations after the black box.

The method developed to create the black box state of 12tk=02t-1|k>|7kmod 15> is a completely literal implementation of that summation of a tensored product of the exponent register with the black box result. The 11 qubits (Hadamard gated) were made into a tensored long quantum state object of 211=2048 state elements. State |k> for each summation iteration would be created from the full tensored state by selecting the kth element of the full state and surrounding it with all other elements set to zero. The only remaining task would be to generate the modular exponential value, load the four work qubits with its binary representation, and form the 16 element tensored quantum state from that. The two quantum states of |k> and |7kmod 15> could then be tensor multiplied to form a 32768 element state to add to the 32768 element total sum state.

There was still one little hitch to that plan: ordinary arithmetic with floating point numbers could not be used to generate all the modular exponent values. The lowest order mantissa bit of the double type falls above 20 beyond 718, and the lowest order mantissa bit of the long double type falls above 20 beyond 722. This will clearly not make it anywhere near 72047.

Fortunately, the very reason that modular exponentiation is used to search for factors is that it generates numbers in repeating periods. One merely needs to generate one period of the numbers and just re-apply them to the larger and larger |k> vectors beyond the first period (although the periods had better be shorter than 24: the 24th number would have to be generated by 723, which is out of range!).

7kmod 15 has a period of 4: { 1, 7, 4, 13 }, so those numbers were made into tensored qubits states loaded into a four element array of them. The index to the array was treated in a modular fashion as |k> went up and the indexed element was tensor multiplied with each |k> vector to generate each summed vector.

The summation finally generated the first intermediate result to check against Box 5.4’s example:

sqrt(1/2048)|1> + sqrt(1/2048)|23> + sqrt(1/2048)|36> + sqrt(1/2048)|61> + …

… + sqrt(1/2048)|32705> + sqrt(1/2048)|32727> + sqrt(1/2048)|32740> + sqrt(1/2048)|32765>

…which compares with:

12t[|0>|1> + |1>|7> + |2>|4> + |3>|13> + |4>|1> + |5>|7> + |6>|4> + …] (10)

The printout actually does match the formal expression above from Box 5.4. One gets the printed basis numbers if the |k> basis numbers from the text are multiplied by 16 and the |7kmod 15> basis numbers from the text are added to them.

The next step was to recreate the change to the summed state due to measurement of the work register. Box 5.4 invokes implicit measurement (11) to accomplish this at this point in the circuit. Following the choice in the text to declare |4> as the |7kmod 15> measurement result, another special method of the long quantum state class condenses the 32768 element sum by a factor of 16, and a further method adds the sizes of the non-selected elements to the selected elements to effect the choice of |4> for the work register.

The printed result is:

… + sqrt(1/512)|2034> + sqrt(1/512)|2038> + sqrt(1/512)|2042> + sqrt(1/512)|2046>

...which compares with:

42t[|2> + |6> + | 10> |14> + …] (12) Perfect agreement.

The next step leaps completely out of the software to create a set of qubits to feed to the inverse Fourier transform, as a method to de-tensor the above state back into qubits would be problematic in terms of time to explore its development. The time was better spent writing and testing the inverse transform based on the circuit diagram for the Fourier transform (13).

The set of qubits equivalent to the measured state was worked out on paper and found to be

H(|0>)x9|1>|0> . This was created in software and fed to the inverse Fourier transform. The result was tensored to a new long quantum state and the result turned out to be:

sqrt(1/4)|0> - sqrt(1/4)|512> + sqrt(1/4)|1024> - sqrt(1/4)|1536>

This is exactly the graph of the probability distribution presented for the inverse Fourier transform in Box 5.4. The 1536 basis is obtainable as a measurement outcome, which Box 5.4 runs through the classical numeric processes of continued fractions and GCD testing, resulting in finding 3 and 5 as the factors of 15.

The software has been coaxed to a state of success.

Conclusion

The original ideal goal was to have a complete system of data objects that could simply run straight through the results of quantum circuits. Two major issues thwarted this. A very major one was the lack of a quantum circuit design to perform the modular exponentiation with 11 exponent qubits and four work qubits. No circuit design for modular exponentiation with any number of qubits was evident in the Nielsen & Chuang text. Admittedly, further research to find a sample design in literature outside the textbook was not undertaken. The struggle with basic representational issues for the software system and the availability of summing products of state vectors for a solution were substituted for any outreach effort. The uncertainty of how expensive it would be against deadline to translate any found literature into working code and data structure was weighted against it.