Example 5: Shor's Algorithm
Purpose
This notebook uses TinyQsim to demonstrate the principle of Shor's Algorithm for integer factorisation. It makes use of the Quantum Fourier Transform (QFT) and Quantum Phase Estimation (QPE) algorithms that were explored in previous examples, so it would be useful to become familiar with those examples first.
Introduction
Factoring large composite numbers is an interesting mathematical problem in its own right, but it is has particular practical implications for the widely-used RSA encryption system. RSA relies for its security on the fact that it is easy to multiply two large prime numbers but very difficult to find the factors if you are just given their product.
In 1994, Peter Shor published a quantum algorithm, now known as Shor's algorithm, that can theoretically factor large composite numbers exponentially faster than a classical computer. Daniel Simon had previously shown the possibility of exponential speedup on a contrived problem, but Shor's algorithm was the first to offer the promise of exponential speedup on a real-world problem. To have any practical advantage over classical factoring methods, Shor's algorithm would need to be run on a quantum computer with thousands of logical (error free) qubits, which is still some way in the future.
Shor's algorithm is a hybrid algorithm that combines both classical and quantum computation. The top-level routine is classical and calls a quantum order-finding subroutine that runs on a quantum computer.
We will start by looking at how the factoring process works and initially demonstrate it using a classical implementation of the order-finding part. We will then implement a simulation of quantum order-finding using TinyQsim and demonstate the factorisation of some small composite numbers using the quantum algorithm.
Imports and Definitions
The example will make use of the following Python imports and definitions. These include the Quantum Fourier Transform (QFT) and Quantum Phase Estimation (QPE) code from previous example notebooks.
import random
from fractions import Fraction
from functools import partial
from math import gcd, log2, ceil
import numpy as np
from examples_lib import create_qpe_circuit, continued_fraction, convergents
from tinyqsim.utils import bits_to_int
MAX_TRIES = 10 # Give up after 10 attempts
def format_bits(bitlist, nbits=0):
"""Format a list of bits (0 or 1) as a string, optionally padded with zeros to nbits.
Example: format_bits([0,1,0,1],6) -> '000101'
"""
return ''.join([str(i) for i in bitlist]).zfill(nbits)
The Top-Level Algorithm
The algorithm is based on number theory and makes use of modular arithmetic as well as Euclid's algorithm for efficiently finding the Greatest Common Divisor (GCD) of two integers.
The aim is to factor a composite number \(N\) into its factors \(p_1\) and \(p_2\):
\[N=p_1\ p_2 \]We will assume that these factors are odd primes. If \(N\) is even, then \(2\) is trivially a factor. If the values are not prime, we can recursively call the factorization routine, using a classical primality test to decide when to stop recursing.
We start by choosing a random integer \(a\) such that \(1<a<N\).
If \(\gcd(a,N)\ne 1\) then the GCD is a factor and we can stop.
Otherwise, we find the multiplicative order of \(a\) modulo \(N\). This is defined as the smallest positive integer \(r\) such that:
\[a^r\equiv 1\pmod{N} \]where \(a\) is an integer coprime to \(N\):
\[\gcd(a, N)=1 \]Finding the order is the part of the algorithm that can make use of a quantum processor to give an exponential speed improvement.
If \(r\) is odd we restart the algorithm with a new random value of \(a\). If \(r\) is even, then \(\frac{r}{2}\) is an integer and we have found a square root of \(1\) modulo \(N\):
\[(a^{r/2})^2 \equiv 1\pmod{N} \]By the difference of squares:
\[(a^{r/2})^2 - 1^2 \equiv 0\pmod{N} \] \[(a^{r/2}-1)(a^{r/2}+1)\equiv 0\pmod{N} \]Hence:
\[N\ |\ (a^{r/2}-1)(a^{r/2}+1) \]If \(N\nmid (a^{r/2}\pm 1)\) then the following must be factors of \(N\) and we have a solution:
\[\qquad\gcd(a^{r/2}+1, N)\quad\text{ and }\quad\gcd(a^{r/2}-1, N) \]Note: Modular exponentiation and Euclid's algorithm can be used to calculate expressions of this form efficiently.
We require that neither of the terms \((a^{r/2}+1)\) and \((a^{r/2}-1)\) is zero. The latter can't be zero because it would imply that \(r/2\) was the order of \(a\) modulo \(N\); but it is smaller than \(r\) which is defined as the smallest positive integer \(r\) such that \(a^r\equiv 1\pmod{N}\), which would lead to a contradiction. To guard against the first one being zero, we restart the algorithm with a new random \(a\) if:
\[\quad a^{r/2}\equiv -1\pmod{N} \]It can be seen from the above that the algorithm is probabilistic in nature and may fail at at several points if certain conditions are not met. In these cases, the algorithm is restarted with a new random value of \(a\). The condition that \(r\) must be even has a probability of about 50%. The conditional probability that \(a^{r/2}\not\equiv -1\pmod{N}\), given that \(r\) is even, is at least 50%. From this simplistic analysis, the probability of success for a particular initial guess is \(\approx\)25%, which implies that about 4 tries are needed on average.
When the order-finding is run on a quantum processor, the speed is exponentially greater than factoring on a classical computer. Even if a few runs are needed because of the restarts, the overall speedup is still exponentially better.
Simple Numerical Example
To demonstrate how the algorithm works, we will factor \(N=21\).
Let us assume \(a=2\) as the random starting value, chosen so that \(1<a<N\).
First, calculate \(\gcd(2,21)=1\)
As the GCD is 1, it follows that 2 and 21 are coprime, so there isn't a trivial factor and we can proceed with the "difference of squares" approach.
Next, we need to find the order \(r\), which is the smallest positive integer such that \(2^r\equiv 1\pmod{21}\). To do this classically, we could try succesive integers for 'r':
\[\begin{align*} 2^0&\equiv 1&\pmod{21}&\quad\text{(start of cycle)}\\ 2^1&\equiv 2&\pmod{21}\\ 2^2&\equiv 4&\pmod{21}\\ 2^3&\equiv 8&\pmod{21}\\ 2^4&\equiv 16&\pmod{21}\\ 2^5&\equiv 11&\pmod{21}\\ 2^6&\equiv 1&\pmod{21}\\ \end{align*} \]Hence the order \(r=6\). The pattern repeats with a cycle length of 6, with the first cycle starting at \(r=0\).
Since \(r\) is even, we can proceed as follows, using the difference of squares:
\[2^6-1\equiv 0 \pmod{21} \] \[(2^\frac{6}{2}+1)(2^\frac{6}{2}-1)\equiv 0\pmod{21} \] \[9\times 7 \equiv 0\pmod{21} \]\(21\) doesn't divide \(7\) or \(9\), so the GCDs of \(7\) and \(9\) with \(21\) must divide \(21\).
The required factors of 21 are therefore:
\[\gcd(9,21)=3\quad\text{ and }\quad\gcd(7,21)=7 \]Classical Order Finding
The top-level algorithm described above can be coded in Python as follows. To demonstrate that this works as intended, a simple classical implemention of the order-finding routine will be used initially.
def order_classical(a, N):
"""A non-quantum order-finding routine
:param a: An initial random value 1<a<N
:param N: The number to be factored
:return: the multiplicative order
"""
for r in range(1, N):
if pow(a, r, N) == 1:
return r
raise ValueError(f'Failed to find order of {a} modulo {N}')
def factorize_classical(N: int):
"""Factorize N using classical order-finding.
:param N: The number to be factored
:return: Tuple of factors (or None if not found)
"""
for _ in range(MAX_TRIES):
a = random.randint(2, N - 1)
g = gcd(a, N)
if g > 1:
return g, N // g # Trivial factors if a & N are not coprime
r = order_classical(a, N)
if (r % 2 == 0) and (a ** (r // 2)) % N != N - 1:
return gcd(a ** (r // 2) + 1, N), gcd(a ** (r // 2) - 1, N)
return None
The 'for' loop limits the number of restarts to avoid an infinite loop if the factorize function is called with a prime number.
To demonstrate that this top-level routine works as intended, we can try factoring a few numbers:
print(factorize_classical(21))
print(factorize_classical(119))
print(factorize_classical(247))
(3, 7)
(7, 17)
(13, 19)
Quantum Order Finding
Shor's algorithm achieves its performance gain on a quantum computer by using a quantum implementation of order-finding. This actually has both quantum and classical parts:
- A quantum circuit that performs Quantum Phase Estimation (QPE)
- A classical algorithm that extracts the order \(r\) from the phase
We need to find the multiplicative order of \(a\) modulo \(N\), which is defined as the smallest integer such that:
\[a^r\equiv 1\pmod N \]where \(a\) is coprime to \(N\).
Quantum order finding works by defining a unitary \(U\) such that:
\[U: \ket{x}\mapsto \ket{ax\bmod N} \]The QPE returns an approximation to an eigenvalue of \(U\), when the input of the target register is the corresponding eigenvector.
If \(U_{a,N}\ket{x} = \ket{ax \bmod N}\), then it satisfies \(U^r_{a,N}=I\) when \(r\) is the order.
Hence, its eigenvalues will be \(e^{2\pi\frac{k}{r}}\) for \(k=0,\dots,r-1\).
Applying \(U\) to the eigenstate \(k\)-times multiplies the phase \(\theta\) by \(k\):
\[U^k\ket{u_e}=e^{2\pi ik\theta}\ket{u_e} \] \[U^k: \ket{x}\mapsto \ket{a^kx\bmod N} \]where \(\ket{u_e}\) is the eigenvector of \(U\) for which we want to find the phase of the eigenvalue \(\lambda\).
We then use Quantum Phase Estimation (QPE) to find the phase as described in the earler notebook on QPE.
Using QPE with \(m\) qubits for the control register, we will measure some value \(x\), such that \(\frac{x}{2^m}\) is an approximation to the true phase \(\frac{k}{r}\):
\[\frac{x}{2^m}\approx\frac{k}{r} \]We need to deduce the order \(r\) by finding the exact rational \(\frac{k}{r}\) from the approximation \(\frac{x}{2^m}\).
Implementing the Unitary
The unitary operator \(U\) is specific to given values of \(a\) and \(N\) and consequently needs to be constructed for each run of the circuit. There are several ways to build such an operation out of gates, but they are quite complicated. Since the purpose of this notebook is just to demonstrate the principle of Shor's algorithm, we will implement the unitary as a single \(2^t\times 2^t\) unitary matrix, where \(t\) is the number of qubits necessary for a binary represententation of the number to be factored. Modular exponentiation can be viewed as a permutation operation, so the matrix can be constructed by starting with an identity matrix and permuting the columns, as shown in the following simple implementation.
def modexp_unitary(a, N, t, k):
"""Create unitary gate for modular exponentiation: a^k mod N.
:param a: An initial random value 1<a<N
:param N: The number to be factored
:param t: Number of bits required to represent N in binary
:param k: The exponent
:return: the unitary matrix
"""
assert gcd(a, N) == 1 # Pre-condition: 'a' is coprime with N
m = np.eye(2 ** t, dtype=int)
m[[(a ** k * i) % N for i in range(N)]] = m[range(N)]
return m
Quantum Phase Estimation
The following Python function constructs and runs the quantum circuit for performing Quantum Phase Estimation using the modular exponentiation unitary. See the previous notebook on Quantum Phase Estimation for more background and the definition of the function 'create_qpe_circuit'.
def run_quantum_circuit(a, N, m=None):
"""Create and run the quantum phase estimation circuit.
:param a: An initial random value 1<a<N
:param N: The number to be factored
:param m: Overrides the default number of control qubits
:return: (circuit, bits)
"""
# Calculate register sizes
t = int.bit_length(N) # Size of state register (to hold N)
if not m:
m = ceil(2 * log2(N)) + 1 # Size of control register
print(f'N={N}, m={m}, t={t}')
# Create the QPE circuit
u = partial(modexp_unitary, a, N, t)
qc = create_qpe_circuit(u, m, t, eigenvec=1)
# Measure the state
qc.barrier('3')
bits = qc.measure(*range(m))
return qc, bits
We can now build the quantum circuit for specific values of \(a\) and \(N\).
The following diagram shows the general form of the Quantum Phase Estimation circuit for the trivial case of factoring 15. The number of qubits used for the 'C' register has been reduced from 9 to 4 to make the diagram easier to understand.
The circuit is the same as the example in the QPE notebook except that it is using the modular exponentiation unitary and measurements are performed on the output of the control register.
qc, bits = run_quantum_circuit(a=11, N=15, m=4)
qc.draw()
N=15, m=4, t=4
The qubits are split into two registers, C and T. Register C is the control register, with \(m\) qubits, and register T is the target register, with \(t\) qubits, that is initialized with the eigenvector.
The operation of this circuit was discussed in detail in the notebook on Quantum Phase Estimation (QPE), so it will not be repeated here.
For Shor's algorithm, register T requires enough qubits to hold the binary representation of the number \(N\) that is to be factored:
\[t=\lceil\log_2 N\rceil \]Register C normally requires more qubits than register T to give the required phase resolution. The following value is commonly used:
\[m = \lceil 2\log_2 N\rceil + 1 \]Note: This value can be derived using Legendre's theorem on continued fractions. It guarantees that the correct partial-fraction convergent will appear in the list of convergents (see above). However, the convergents are reduced fractions, so the required order may be a multiple of the convergent denominator.
The circuit diagram above is just to illustrate the generate form of the QPE circuit, so a smaller value of \(m\) has been used to make the diagram easier to understand.
The target register is initialized to the eigenvector \(\ket{u_e}=\ket{0\dots 01}\) of the unitary \(U\). Then applying QPE gives an approximation \(\frac{x}{2^m}\) to the phase, as discussed above.
Computing the Order
Using Quantum Phase Estimation (QPE), measurement of the \(m\)-qubit control register will give some value \(x\), such that \(\frac{x}{2^m}\) is an approximation to the true phase \(\theta=\frac{k}{r}\):
\[\frac{x}{2^m}\approx\frac{k}{r} \]where \(k=0,1,\dots,r-1\)
The next step is to deduce the multiplicative order \(r\) by finding the exact rational \(\frac{k}{r}\) from the approximation \(\frac{x}{2^m}\). This can be done by computing the sequence of continued-fraction convergents, which are a sequence of progressively better approximations to a fraction.
We just need to check the denominator \(q\) of each convergent until we find the first one that satisfies:
\[a^q\equiv 1\pmod{N} \]For example, for \(a=40, N=119, m=15\), the QPE might return \(x=27989\). Then the continued-fraction convergents would be:
\[\small\frac{1}{1}, \frac{5}{6}, \frac{6}{7}, \frac{35}{41}, \frac{41}{48} \]Checking each of the denominators \(q\) as a candidate \(r\), the first to succeed is \(r=48\):
\[a^r = 40^{48} \equiv 1 \pmod{119} \]Hence, we have found that the order \(r\) is \(48\).
The approximation error in the value returned by QPE is:
\[\left|\frac{x}{2^m}-\frac{k}{r}\right| \]Legendre's theorem on continued fractions states that if a rational number \(\frac{p}{q}\) is a sufficently good approximation of a real number \(\alpha\), such that:
\[\left|\alpha-\frac{p}{q}\right| < \frac{1}{2q^2} \]then \(\frac{p}{q}\) must be one of the convergents of the continued fraction of \(\alpha\).
Rather than test each convergent denominator to see whether it is the order 'r' of \(a\) mod \(N\), it can be faster to only test ones that satisfy the approximation error condition:
\[\left|\frac{x}{2^m}-\frac{k}{r}\right| < \frac{1}{2r^2} \]The procedure may still sometimes fail to find the order because the convergents are reduced fractions in their lowest terms. The actual order may therefore be a multiple of the convergent denominator. The simplest solution is to restart the algorithm with a new random value \(a\), as is done for other failure conditions.
Note: Another solution to the reduced fraction problem is to try small multiples of the candidate order. Alternatively, the quantum circuit may be run two or more times and the Least Common Multiple (LCM) of the resulting denominators used as the candidate \(r\) value.
The following function is used to compute the sequence of convergents from the phase returned by the Quantum Phase Estimation circuit.
def build_convergents(x):
"""Return the continued-fraction convergents of the fraction 'x'.
:param x: The fraction
:return: list of tuples (p,q) for the convergents
"""
f = Fraction(x)
fraction = continued_fraction(f)
return list(convergents(fraction))
The order is then found by checking each convergent to find the first one whose denominator \(q\) satisfies:
\[a^q \equiv 1 \pmod{N} \]The check based on Legendre's theorem on continued fraction, which was mentioned earlier, is used to skip over convergents that are not sufficiently good approximations.
This may fail as the convergents are reduced fractions, in which case small multiples of the candidate convergent are tried.
def order_from_measured_bits(bits, a, N, multiples=8):
"""Compute the multiplicative order from bits returned by QPE.
:param bits: QPE result as an array of bit values (0 or 1)
:param a: An initial random value 1<a<N
:param N: The number to be factored
:param multiples: Largest multiplier to be tried
:return: the multiplicative order (or None)
"""
m = len(bits)
x = bits_to_int(bits)
if x == 0:
return None
M = 2 ** m
# Find convergents
f = Fraction(x, M)
convergents = build_convergents(f)
print('Convergents:', [f'{p}/{q}' for p, q in convergents if q < N])
for p, q in convergents:
if q > N:
return None
if q == 1:
continue
# Check whether convergent is a good enough approximation
ok = abs(x / M - p / q) <= 1 / (2 * 2 ** m)
print(f'Legendre test on convergent: {p}/{q} -> {ok}')
if ok:
# Try small multiples as convergents are reduced fractions
for i in range(1, multiples + 1):
print(f'Check: {a}^({q}*{i}) mod {N} = {pow(a, q * i, N)}')
if pow(a, q * i, N) == 1:
return q * i
return None
Putting the above functions together, we can define the quantum order-finding routine. This involves generating and running the quantum circuit and then extracting the order from the measured phase.
def order_quantum(a, N):
"""Quantum order-finding function.
:param a: An initial random value 1<a<N
:param N: The number to be factored
:return: The multiplicative order 'r'
"""
# Run the quantum circuit
qc, bits = run_quantum_circuit(a, N)
bits = qc.results().values() # FIXME: Assumes dict is ordered
print(f'bits: {format_bits(bits, len(bits))} = {bits_to_int(bits)}')
# Convert measurement result to phase
return order_from_measured_bits(bits, a, N)
Finally, we can put everything together into a top-level factoring function as follows. The output includes the order found by a classical algorithm for comparison.
def shor_factorize(N, n_tries=MAX_TRIES):
"""Factor N using Shor's Algorithm.
:param N: The number to be factored
:param n_tries: Maximum number of attempts
:return: Tuple of factors (or None if not found)
"""
# Attempt to find the order r of a^x mod N.
for trial in range(1, n_tries + 1):
print(f'\nTrial # {trial}')
a = random.randint(2, N - 1)
print(f'Random a={a}')
g = gcd(a, N)
if g > 1:
print('Success: a and N are not coprime')
return g, N // g # Trivial factors if a & N are not coprime
r = order_quantum(a, N)
print(f'Order r (classical) = {order_classical(a, N)}') # For comparison
print(f'Order r (quantum) = {r}')
if not r:
print('Retry as convergent not found')
continue
if r % 2 != 0:
print('Retry as r is odd')
continue
if (a ** (r // 2)) % N == N - 1:
print(f'Retry as a**r//2 mod N = N-1')
continue
print(f'Success with difference of squares')
return gcd(a ** (r // 2) + 1, N), gcd(a ** (r // 2) - 1, N)
print(f'Failed to find a solution after {n_tries} attempts')
return None
Demonstration
The following demonstration shows \(119\) being factored using the quantum algorithm.
shor_factorize(119)
Trial # 1
Random a=19
N=119, m=15, t=7
bits: 001010101010101 = 5461
Fraction = 5461/32768
Convergents: ['0/1', '1/6']
Legendre test on convergent: 1/6 -> True
Check: 19^(6*1) mod 119 = 64
Check: 19^(6*2) mod 119 = 50
Check: 19^(6*3) mod 119 = 106
Check: 19^(6*4) mod 119 = 1
Order r (classical) = 24
Order r (quantum) = 24
Success with difference of squares
(17, 7)
Shor's algorithm can, in principle, give an exponential speed-up over a classical computer. However, when the quantum algorithm is simulated on a classical computer, it is very much slower than a classical factorization algorithm.
You can try factoring other numbers but the algorithm will become very slow for larger values of \(N\). Factoring 119, which is a 7-bit number, requires \(m+t=15+7=22\) qubits and may take a few seconds. Factoring 253, which is an 8-bit number, requires \(m+t=17+8=25\) qubits and will consequently take longer to run.
Example Output
Because of the probabilistic nature of the algorithm, there may be a number of retries and different paths to the solution. If the random variable \(a\) is not coprime with \(N\), it can return the result without even using the quantum algorithm.
The following snapshots show some example outputs that may occur when factoring 119.
Example A
Trial # 1
Random a=40
N=119, m=15, t=7
bits: 110110101010101 = 27989
Fraction = 27989/32768
Convergents: ['0/1', '1/1', '5/6', '6/7', '35/41', '41/48']
Legendre test on convergent: 5/6 -> False
Legendre test on convergent: 6/7 -> False
Legendre test on convergent: 35/41 -> False
Legendre test on convergent: 41/48 -> True
Check: 40^(48*1) mod 119 = 1
Order r (classical) = 48
Order r (quantum) = 48
Success with difference of squares
(17, 7)
This is the nominal situation in which the algorithm finds a convergent whose denominator is the required order.
\[a^r \bmod{N} = 40^{48} \bmod{119} = 1 \]Example B
shor_factorize(119)
Trial # 1
Random a=9
N=119, m=15, t=7
bits: 111010101010101 = 30037
Fraction = 30037/32768
Convergents: ['0/1', '1/1', '10/11', '11/12']
Legendre test on convergent: 10/11 -> False
Legendre test on convergent: 11/12 -> True
Check: 9^(12*1) mod 119 = 50
Check: 9^(12*2) mod 119 = 1
Order r (classical) = 24
Order r (quantum) = 24
Success with difference of squares
(17, 7)
In this example, the convergent \(11/12\) satisfies the appproximation criterion, but 12 is not the order because the convergent is a reduced fraction.
\[9^{12\times 1} \bmod{119} = 50 \]Small multiples of 12 are then tried, with a multiple of 2 giving the required order 24:
\[9^{12\times 2} \bmod{119} = 1 \]Example C
shor_factorize(119)
Trial # 1
Random a=55
N=119, m=15, t=7
bits: 000000000000000 = 0
Order r (classical) = 4
Order r (quantum) = None
Retry as convergent not found
Trial # 2
Random a=68
Success: a and N are not coprime
(17, 7)
In this particular example, the cause of the failure was that the measured phase was zero. This value has a non-zero probability in the QPE, which leads to a failure.
The algorithm starts again with a different random value \(a=68\) which happens not to be coprime to \(N\), so one factor is \(\gcd(119,68)=17\) and the other is \(119/17=7\).