Factoring 15 with Shor's Algorithm
The integer factorization problem [1] is a famous problem in number theory: given a composite number \(N\), find its prime factors. The importance of the problem stems from there being no known efficient (polynomial-time, in the number of bits needed to represent \(N\)) classical algorithm, and much of modern-day cryptography relies on this fact. In 1994, Peter Shor came up with an efficient quantum algorithm for the problem [2], providing one of the first concrete pieces of evidence for the power of quantum computers.
Shor's Algorithm
Shor's algorithm consists of a classical part and a quantum subroutine. The steps of the algorithm for factoring an input number \(N\), summarized from [3]:
- Pick a random number \(1 < a < N\) that is co-prime with \(N\). Co-primality can be checked by computing the GCD (greatest common divisor) of \(a\) and \(N\). If it is 1 then we have a co-prime \(a\), otherwise we have a non-trivial factor of \(N\) and we are done.
-
Find the period \(r\) of the following function, using the quantum period finding algorithm (described in [4]): \(\(f(x) = a^x \mod N\)\)
-
If \(r\) is odd or \(a^{r/2} = -1 \mod N\), return to step 1 (this event can be shown to happen with probability at most \(1/2\)).
- Otherwise, \(\gcd(a^{r/2} \pm 1, N)\) are both factors of \(N\), and computing one of them yields the required result.
This demo factors the number \(N=15\) using Shor's algorithm by applying the quantum subroutine (step 2) with \(a=7\). We choose this particular \(a\) since it is co-prime with 15 and satisfies the conditions of step 3, providing us with a high probability of finding a factor of \(N\).
Building the Quantum Period Finding Circuit
We begin by declaring the number of qubits in the upper (counting) register that the quantum subroutine uses. In our case, \(N = 15\), and according to the algorithm the upper register must contain \(q = \log(Q)\) qubits for \(Q\) such that \(N^2 \le Q < 2N^2\); namely, \(225 < Q < 450\), and therefore \(q = 8\). In addition, the second register should be large enough to encode 15, hence:
import numpy as np
N = 15
num_counting_qubits = int(np.ceil(np.log2(N**2)))
num_auxilliary_qubits = int(np.ceil(np.log2(N)))
We implement a Phase Estimation [5] circuit. Each element in the circuit is a controlled operation of \(\(|x\rangle \rightarrow |x\cdot a^{2^i}\mod N \rangle\)\) where \(a < N\) is a number such that \(\gcd(a, N)=1\). For this demonstration we pick \(a=7\). \(i\) is the index of the control qubit, located in the upper register.
It is quite involved to implement these unitaries, so for this demo we take a shortcut and compute the exact unitary matrix that implements the computation (which in general is not applicable as this preprocessing is exponential). We do so by calculating the modular-multiplication by \(a\) matrix and then using its powers.
The unitary function is used for decomposing the unitary matrix into quantum gates.
from classiq import *
def get_modular_multiplication_matrix():
# fmt: off
swap = np.array(
[
[1, 0, 0, 0],
[0, 0, 1, 0],
[0, 1, 0, 0],
[0, 0, 0, 1]
],
dtype=complex
)
# fmt: on
swap32 = np.kron(np.identity(4), swap)
swap21 = np.kron(np.kron(np.identity(2), swap), np.identity(2))
swap10 = np.kron(swap, np.identity(4))
x = np.array([[0, 1], [1, 0]])
x_all = np.kron(np.kron(x, x), np.kron(x, x))
u = x_all @ swap10 @ swap21 @ swap32
return u
MODULAR_MUL_UNITARY = get_modular_multiplication_matrix().real.tolist()
@qfunc
def modular_exponentiation(
exponent: CInt, target: QArray[QBit, num_auxilliary_qubits]
) -> None:
power(2**exponent, lambda: unitary(elements=MODULAR_MUL_UNITARY, target=target))
Complete Circuit
At the first layer of the quantum circuit, we prepare the equal superposition state in the top (counting) register, and prepare the \(|1\rangle\) state in the bottom (auxiliary) register.
We then apply the second layer of the circuit, consisting of the controlled \(U^{2^i}\) gates. Lastly, we apply an inverse QFT on the counting register to get the period.
@qfunc
def period_finding(
qv_counting: Output[QArray[QBit, num_counting_qubits]],
qv_auxilliary: Output[QArray[QBit, num_auxilliary_qubits]],
) -> None:
# start with a Hadamard transform in the counting register
allocate(num_counting_qubits, qv_counting)
hadamard_transform(qv_counting)
# Prepare the |1> state on the lower register
allocate(num_auxilliary_qubits, qv_auxilliary)
X(qv_auxilliary[0])
# Apply the controlled modular-exponentiations using each of the counting qubits
repeat(
count=num_auxilliary_qubits,
iteration=lambda index: control(
ctrl=qv_counting[index],
stmt_block=lambda: modular_exponentiation(index, qv_auxilliary),
),
) # ! not working with qv[a:]
# Lastly, apply an inverse QFT
invert(lambda: qft(qv_counting))
Quantum Entry Point
To synthesize the circuit, we define a quantum main function. As we are only interested in the output of the counting register, we only define it in the signature of the function.
Next, we translate it to Qmod using create_model.
@qfunc
def main(qv_counting: Output[QArray[QBit, num_counting_qubits]]) -> None:
qv_auxilliary = QArray("qv_auxilliary")
period_finding(qv_counting, qv_auxilliary)
qmod = create_model(main, out_file="shor")
We now send the model to the synthesis engine, which may take a few seconds:
qprog = synthesize(qmod)
We can now view the circuit and its depth:
show(qprog)
Executing the Circuit
Now, we execute the circuit above, using the simulator:
result = execute(qprog).result_value()
import collections
hist_counting_qubits = collections.defaultdict(int)
for key, value in result.counts.items():
hist_counting_qubits[key] += value
Plotting the result:
import matplotlib.pyplot as plt
plt.bar(hist_counting_qubits.keys(), hist_counting_qubits.values())
<BarContainer object of 4 artists>
We obtained four \(y\) results from the circuit, each with a rough probability of \(1/4\): \(0, 64, 128\), and \(192\). By dividing by \(Q = 256\) we obtain four reduced fractions: \(0, 1/4, 1/2\), and \(3/4\), with two of them having the correct period \(r=4\) in the denominator. With this period, we can compute the factors of \(N = 15\): \(\gcd(a^{r/2} \pm 1, N) = \gcd(7^2 \pm 1, 15) = 3, 5\).
References
[1]: Integer Factorization (Wikipedia)
[2]: Shor, Peter W. "Algorithms for quantum computation: Discrete logarithms and factoring." Proceedings 35th annual symposium on foundations of computer science. IEEE, 1994.
[3]: Shor's Algorithm Procedure (Wikipedia)