Factoring 15 with Shor's Algorithm¶
Introduction¶
The integer factorization problem [1] is a famous problem in number theory: given a number $N$ which is composite, find its prime factors. The importance of the problem stems from the fact that no efficient (polynomial-time, in the number of bits needed to represent $N$) classical algorithm is known for it to this day, 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], are as follows:
- 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 found a co-prime $a$, otherwise we have found 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.
In this demo, we will factor the number $N=15$ using Shor's algorithm, by applying the quantum subroutine (step 2) with $a=7$. This particular $a$ is chosen 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 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 will 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 picked $a=7$. $i$ is the index of the control qubit, located in the upper register.
It is quiet involved to implement these unitaries, so for this demo we will make a shortcut, and compute exactly the unitary matrix that implements the computation (which in the general case is not applicable as this pre-processing is exponential). We will do so by calculating the modular-multiplication by $a$ matrix, then using its powers.
The function unitary is used for decomposing the unitary matrix into quantum gates.
import sympy
from classiq import (
H,
Output,
QArray,
QBit,
QFunc,
QParam,
X,
allocate,
control,
create_model,
invert,
join,
qft,
repeat,
split,
uniform_superposition,
unitary,
)
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(
power: QParam[int], target: QArray[QBit, num_auxilliary_qubits]
) -> None:
unitary(elements=sympy.Matrix(MODULAR_MUL_UNITARY) ** (2**power), target=target)
Building the 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, which consists 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 uniform superposition in the counting register
allocate(num_counting_qubits, qv_counting)
uniform_superposition(num_counting_qubits, qv_counting)
# Prepare the |1> state on the lower register
allocate(num_auxilliary_qubits, qv_auxilliary)
X(qv_auxilliary[0])
# Apply the contolled modular-exponentiations using each of the counting qubits
qv_all = QArray("qv_all")
join(qv_counting, qv_auxilliary, qv_all)
repeat(
count=num_auxilliary_qubits,
iteration=lambda index: control(
operand=lambda target: modular_exponentiation(index, target),
ctrl=qv_all[index],
target=qv_all[
num_counting_qubits : num_counting_qubits + num_auxilliary_qubits
],
),
) # ! not working with qv[a:]
split(
num_counting_qubits,
num_auxilliary_qubits,
in_=qv_all,
out1=qv_counting,
out2=qv_auxilliary,
) #! without kwargs does not work
# Lastly, apply an inverse QFT
invert(lambda target: qft(num_qbits=num_counting_qubits, qbv=target), qv_counting)
Quantum entry point¶
In order to synthesize the circuit, we define a quantum main function. As are we 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 the 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)
with open("shor.qmod", "w") as f:
f.write(qmod)
We now send the model to the synthesis engine, taking a few seconds:
from classiq import synthesize
qprog = synthesize(qmod)
We can now view the circuit and its depth:
from classiq import show
show(qprog)
Executing the circuit¶
Now, we turn to executing the circuit above, using the simulator:
from classiq import execute
results = execute(qprog).result()
res = results[0].value
import collections
hist_counting_qubits = collections.defaultdict(int)
for key, value in res.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())
We obtained 4 results $y$ from the circuit, each with probability roughly $1/4$: $0, 64, 128$ and $192$. Dividing by $Q = 256$ we obtain 4 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)