Shor’s Algorithm [1]

Modular Exponentation

Consider a large composite integer number \(N\) with two distinct primes \(p\) and \(q\), such that \(N = p * q\). Then, according to the number theory, there exists a periodic function \(\mathcal{F}(x)\) with period \(r\), such that, \(\mathcal{F}(x) = a^x\hspace{0.2 em}mod\hspace{0.2 em} N\) ; where \(a\) is an integer, \(a ≠ N\), not prime with \(p\) or \(q\) and coprime with \(N\). The period \(r = lcm \left( \left( p-1 \right),\left( q-1 \right) \right)\) is defined by Carmichael’s Totient function. Here \(lcm\) stands for least common multiple.

Let us illustrate this with an example.

Let \(N = 21\), \(p = 7\), \(q= 3\)

Then, \(N = p×q = 7×3 = 21\)

\(r = lcm((7-1),(3-1)) = lcm(6,2) = 6\)

Let, \(a = 5\). Consider the following table:

\(\mathcal{F}(x_0)\)	\(\mathcal{F}(x_1)\)	\(\mathcal{F}(x_2)\)
\(\mathcal{F}(1)=5^1\hspace{0.2 em} mod\hspace{0.2 em} 21=5\)	\(\mathcal{F}(7)=5^7\hspace{0.2 em} mod\hspace{0.2 em} 21=5\)	\(\mathcal{F}(13)=5^{13}\hspace{0.2 em} mod\hspace{0.2 em} 21=5\)
\(\mathcal{F}(2)=5^2\hspace{0.2 em} mod\hspace{0.2 em} 21=4\)	\(\mathcal{F}(8)=5^8\hspace{0.2 em} mod\hspace{0.2 em} 21=4\)	\(\mathcal{F}(14)=5^{14}\hspace{0.2 em} mod\hspace{0.2 em} 21=4\)
\(\mathcal{F}(3)=5^3\hspace{0.2 em} mod\hspace{0.2 em} 21=20\)	\(\mathcal{F}(9)=5^9\hspace{0.2 em} mod\hspace{0.2 em} 21=20\)	\(\mathcal{F}(15)=5^{15}\hspace{0.2 em} mod\hspace{0.2 em} 21=20\)
\(\mathcal{F}(4)=5^4\hspace{0.2 em} mod\hspace{0.2 em} 21=16\)	\(\mathcal{F}(10)=5^{10}\hspace{0.2 em} mod\hspace{0.2 em} 21=16\)	\(\mathcal{F}(16)=5^{16}\hspace{0.2 em} mod\hspace{0.2 em} 21=16\)
\(\mathcal{F}(5)=5^5\hspace{0.2 em} mod\hspace{0.2 em} 21=17\)	\(\mathcal{F}(11)=5^{11}\hspace{0.2 em} mod\hspace{0.2 em} 21=17\) `	\(\mathcal{F}(17)=5^{17}\hspace{0.2 em} mod\hspace{0.2 em} 21=17\)
\(\mathcal{F}(6)=5^6 mod\hspace{0.2 em} 21=1\)	\(\mathcal{F}(12)=5^{12}\hspace{0.2 em} mod\hspace{0.2 em} 21=1\)	\(\mathcal{F}(18)=5^{18}\hspace{0.2 em} mod\hspace{0.2 em} 21=1\)

From the table it is clear that the period \(r = 6\), as found earlier.

Modular Exponentation of 15

To recap, the modular exponentiation calculates the remainder, when an integer \(a\) is raised to the \(x\) th power (i.e., \(a^x\)) is divided by a positive integer \(N\), called the modulus. Using symbols, we can write this as \(c= a^x mod\hspace{0.2 em} N\). Assuming \(a = 2\) and \(N = 15\); the following table illustrates this operation:

\(a^x \hspace{0.2 em}mod\hspace{0.2 em} N\)	\(c\)
\(2^0 \hspace{0.2 em}mod\hspace{0.2 em} 15\)	\(1\)
\(2^1 \hspace{0.2 em}mod\hspace{0.2 em} 15\)	\(2\)
\(2^2 \hspace{0.2 em}mod\hspace{0.2 em} 15\)	\(4\)
\(2^3 \hspace{0.2 em}mod\hspace{0.2 em} 15\)	\(8\)

This table can be verified by calculating \(a^x\) first, and then by performing \(mod\hspace{0.2 em} N\), using a calculator.

If we can state \(x\) in its binary form \((x_{n-1} x_{n-2} x_{n-3} \cdots x_0 )\), we can expand the modular exponentiation as follows.

\[a^x\hspace{0.2 em}mod\hspace{0.2 em}N= (a^{2^{0}}\hspace{0.2 em}mod\hspace{0.2 em} N )^{x_0} \cdot (a^{2^{1}}\hspace{0.2 em}mod\hspace{0.2 em} N )^{x_1}\cdots (a^{2^{n-1}} \hspace{0.2 em}mod \hspace{0.2 em} N )^{x_{n-1}} \hspace{0.2 em}mod\hspace{0.2 em} N\]

Modular exponentiation can be constructed as a circuit using controlled multiplier blocks, as shown in the following figure.

In this circuit, each of the controlled multiplier blocks performs a multiplication of the input REG 2, by a constant function \(a^{2{^i}}\hspace{0.2 em} mod\hspace{0.2 em} N\).Each of these controlled multiplier blocks is controlled by the respective qubit \(x_i\) in REG 1. If we assume \(|x \rangle\) as the control qubits (REG 1) and \(|y\rangle\) as the target qubits (REG 2), this operation produces:

\[CU_{a^{2{^i }}} |x\rangle |y\rangle = |x\rangle (a^{2{^i}} )^x y \hspace{0.2 em}mod\hspace{0.2 em} N\]

We can explain the inner workings of this circuit with an example that performs \(2^x \hspace{0.2 em}mod\hspace{0.2 em} 15\). Consider the following circuit:

In this circuit, REG 1 is formed by qubits q[0] to q[2]. REG 2 is formed by qubits q[3] to q[6]. We can write the state of the system after application of the Walsh-Hadamard transform and the X gate as follows.

\[\begin{split}|ψ \rangle = \frac{1}{\sqrt{8}} ( |000\rangle|0001\rangle + |001\rangle|0001\rangle + |010\rangle|0001\rangle + \\ |011\rangle|0001\rangle + |100\rangle|0001\rangle + |101\rangle|0001\rangle + \\ |110\rangle|0001\rangle + |111\rangle|0001\rangle)\end{split}\]

The application of the first block of multiplier converts this state as follows.

\[\begin{split}|ψ \rangle = \frac{1}{\sqrt{8}} ( |000\rangle|0001\rangle + |001\rangle|0010\rangle + |010\rangle|0001\rangle + \\ |011\rangle|0010\rangle +|100\rangle|0001\rangle + |101\rangle|0010\rangle + \\ |110\rangle|0100\rangle + |111\rangle|0010\rangle )\end{split}\]

After the application of the second block:

\[\begin{split}|ψ \rangle = \frac{1}{\sqrt{8}} (|000\rangle|0001\rangle + |001\rangle|0010\rangle + |010\rangle|0100\rangle + \\ |011\rangle|1000\rangle +|100\rangle|0001\rangle + |101\rangle|0010\rangle + \\ |110\rangle|0100\rangle + |111\rangle|1000\rangle)\end{split}\]

Application of the final block produces the following state.

\[\begin{split}|ψ \rangle = \frac{1}{\sqrt{8}} (|000\rangle|1000\rangle + |001\rangle|0010\rangle +|010\rangle|0100\rangle + \\ |011\rangle|0001\rangle + |100\rangle|1000\rangle + |101\rangle|0010\rangle + \\ |110\rangle|0100\rangle + |111\rangle|0001\rangle )\end{split}\]

At this stage, if we measure REG 2 (i.e., the qubits q[3] to q[6] in the diagram) , we can find that the system is projected to the following states: 0001, 0010, 0100, 1000 with equal probability. By comparing these values with the previous table, we can conclude that this circuit has performed modular exponentiation.

RSA Security

The following table outlines the innerworkings of RSA as a sequence of steps.

Now, we have the public exponent \(e=19\), and the private exponent \(d=139\). With this, we can try to perform an encryption and decryption operation of the alphabet \(A\). The required steps are outlined in the following table.

Period Estimation

If \(p\) and \(q\) are known, it is easy to find \(r\). But if only \(N\) is known, it is a hard problem. The hard problem of finding the period \(r\) of \(\mathcal{F}(x) = a^x\hspace{0.2 em} mod\hspace{0.2 em} N\) ensures the security of RSA. Once \(r\) is known, it is easy to calculate the Private Exponent \(d\) using the formula \(d = e^{-1}\hspace{0.2 em} mod\hspace{0.2 em} r\), and the RSA security is broken (where \(e\) is the Public Exponent. It is part of the Public Key \(\{ e,N \}\), and usually \(65537\).) Finding the value of \(r\) is called order finding or period estimation. Shor’s algorithm (c.1994) is a Quantum version of period estimation. With this background, we are ready to jump into Shor’s algorithm.

Shor’s Algorithm

In 1994, the American mathematician Peter Shor proposed a novel quantum algorithm for factorizing large integers, which runs on a polynomial-time. In 2001, a group of scientists at IBM factorized the integer 15, using Nuclear Magnetic Resonance (NMR) [2]. In 2012, the integer 21 was factorized, as of 2020, this is the largest number reliably factored using Shor’s algorithm [3]. There are several other claims of large numbers being factorized using classical-quantum hybrid methods.

Shor solved the problem of factorizing integers by shifting the problem domain to finding the period of the modular exponentiation function. Shor’s problem is stated as follows:

Given an integer \(N\), and two prime numbers \(p\) and \(q\) such that \(N = p \times q\). Let \(a \in \mathbb{Z}_n\). Let a modular exponentiation function \(f_a: \mathbb{Z}_+→ \mathbb{Z}_n\) defined as \(f_a (x)= a^x \hspace{0.2 em} mod\hspace{0.2 em} N\) satisfy the condition, \(f(x)⟺f(x_0 ),x=x_0+kr\), where \(r\) is the period of \(f_a\), and \(k\) is a constant. Shor’s problem is determining the period \(r\) by querying \(f_a\).

The block diagram of Shor’s algorithm is given below:

The block diagram uses modular exponentiation and inverse QFT. Inverse QFT is used to find the period r of the equation \(\mathcal{F}(x)≡ a^x \hspace{0.2 em} mod\hspace{0.2 em} N\). The period \(r\) can be used to find the factors of \(N\). The algorithm uses two banks of qubit registers. The upper bank is called the source register, and the lower bank is called the target register. The oracle function performs modular exponentiation on the source register and stores the values in the target register. This process is done in one step using quantum parallelism. The number of qubits \(s\) in the source, and target registers must be in the following range.

\[N^2 ≤ 2^s ≤ (2N)^2\]

From this range, \(\frac{2^s}{r} > N\) and it guarantees that when estimating the period \(r\), \(N\) terms are contributing to the probability amplitude, even if \(r\) is close to \(\frac{N}{2}\). This assumption creates a sufficiently large state-space in the source register, which can contain even the longest periods.By having a large value of \(s\), we can explore a dense solution space, and the period \(r\) can be estimated as accurately as possible.The modular exponentiation circuits require additional qubits to store intermediate information. Hence, we have selected that the target register is also made of \(s\) qubits. Some preprocessing tasks are required before the algorithm is started. First, we need to check whether \(N\) is a prime number.There are some classical methods to do the primality check, and it can be done efficiently on a classical computer. The second task is to identify \(a\) and confirm that \(N\) and \(a\) are coprime.The third task is to determine the number of qubits required. After completing these steps, the circuit can be constructed for execution.

For a complete derivation of Shor’s algorithm, see: [15].

Quantum ciruit to factorize 15

In this section, we construct a quantum circuit and demonstrate that Shor’s algorithm works by factorizing 15. This circuit is based on the work done by Lieven M.K et al. in 2001 using Nuclear Magnetic Resonance [23]. The following circuit finds the period of the function \(\mathcal{F}(x)= 7^x \hspace{0.2 em} mod \hspace{0.2 em} 15\).

This circuit is an optimized version that uses 3-qubits for the source register. Since the period of this function is relatively smaller, the choice of 3-qubits holds good. We start this circuit by preparing the source register in an equal superposition state. Modular exponentiation is done for all the states of the source register and stored in the target register. After this stage, we perform the inverse QFT and measure the source register. The X‍ gate in qubit q[6] in the diagram provides for the phase kickback.

Let us now verify the modular exponentiation circuit in this implementation. The system state after the initialization stage is:

\[\begin{split}|ψ \rangle=\frac{1}{\sqrt{8}} ( |000\rangle|0001 \rangle+ |001\rangle|0001 \rangle+ \\ |010\rangle|0001 \rangle+ |011\rangle|0001 \rangle+ \\ |100\rangle|0001 \rangle+ |101\rangle|0001 \rangle+ \\ |110\rangle|0001 \rangle+ |111\rangle|0001 \rangle)\end{split}\]

After preparing the register to the initial state, we apply a series of CNOT and CCNOT gates to perform the modular exponentiation.

STEP 1: The first CNOT gate uses q[2] as the control qubit and q[4] as the target qubit

\[\begin{split}|ψ \rangle=\frac{1}{\sqrt{8}} ( |000\rangle|0001 \rangle+ |001\rangle|0101 \rangle+ \\ |010\rangle|0001 \rangle+ |011\rangle|0101 \rangle+ \\ |100\rangle|0001 \rangle+ |101\rangle|0101 \rangle+ \\ |110\rangle|0001 \rangle+ |111\rangle|0101 \rangle)\end{split}\]

STEP 2: The second CNOT gate uses q[2] as the control qubit and q[5] as the target qubit.

\[\begin{split}|ψ \rangle=\frac{1}{\sqrt{8}} ( |000\rangle|0001 \rangle+ |001\rangle|0111 \rangle+ \\ |010\rangle|0001 \rangle+ |011\rangle|0111 \rangle+ \\ |100\rangle|0001 \rangle+ |101\rangle|0111 \rangle+ \\ |110\rangle|0001 \rangle+ |111\rangle|0111 \rangle)\end{split}\]

STEP 3: The third CNOT gate with q[6] as the control and q[4] as the target.

\[\begin{split}|ψ \rangle=\frac{1}{\sqrt{8}} ( |000\rangle|0101 \rangle+ |001\rangle|0011 \rangle+ \\ |010\rangle|0101 \rangle+ |011\rangle|0011 \rangle+ \\ |100\rangle|0101 \rangle+ |101\rangle|0011 \rangle+ \\ |110\rangle|0101 \rangle+ |111\rangle|0011 \rangle)\end{split}\]

STEP 4: CCNOT gate with q[1] and q[5] as the control gates and q[3] as the target qubit.

\[\begin{split}|ψ \rangle=\frac{1}{\sqrt{8}} ( |000\rangle|0101 \rangle+ |001\rangle|0011 \rangle+ \\ |010\rangle|0101 \rangle+ |011\rangle|1011 \rangle+ \\ |100\rangle|0101 \rangle+ |101\rangle|0011 \rangle+ \\ |110\rangle|0101 \rangle+ |111\rangle|1011 \rangle)\end{split}\]

STEP 5: The fourth CNOT gate with q[3] as the control and q[5] as the target.

\[\begin{split}|ψ \rangle=\frac{1}{\sqrt{8}} ( |000\rangle|0101 \rangle+ |001\rangle|0011 \rangle+ \\ |010\rangle|0101 \rangle+ |011\rangle|1001 \rangle+ \\ |100\rangle|0101 \rangle+ |101\rangle|0011 \rangle+ \\ |110\rangle|0101 \rangle+ |111\rangle|1001 \rangle)\end{split}\]

STEP 6: CCNOT gate with q[1] and q[4] as the control and q[6] as the target.

\[\begin{split}|ψ \rangle=\frac{1}{\sqrt{8}} (|000\rangle|0101 \rangle+ |001\rangle|0011 \rangle+ \\ |010\rangle|0100 \rangle+ |011\rangle|1001 \rangle+ \\ |100\rangle|0101 \rangle+ |101\rangle|0011 \rangle+ \\ |110\rangle|0100 \rangle+ |111\rangle|1001 \rangle)\end{split}\]

STEP 7: The fifth CNOT gate with q[6] as the control and q[4] as the target.

\[\begin{split}|ψ \rangle=\frac{1}{\sqrt{8}} ( |000\rangle|0001 \rangle+ |001\rangle|0111 \rangle+ \\ |010\rangle|0100 \rangle+ |011\rangle|1101 \rangle+ \\ |100\rangle|0001 \rangle+ |101\rangle|0111 \rangle+ \\ |110\rangle|0100 \rangle+ |111\rangle|1101 \rangle)\end{split}\]

After performing the modular exponentiation, the target register is in equal probability with the states: \(|1 \rangle, |7 \rangle, |4 \rangle, |13 \rangle.\) This result very well aligns with the modular exponentiation shown earlier. We readily see that the period \(r\) is \(4\). Nevertheless, let us continue testing this circuit and validate this later.

\[\begin{split}7^0 \hspace{0.2 em}mod \hspace{0.2 em} 15=1 \\ 7^1 \hspace{0.2 em}mod\hspace{0.2 em} 15=7 \\ 7^2 \hspace{0.2 em}mod\hspace{0.2 em} 15=4 \\ 7^3 \hspace{0.2 em}modv 15=13 \\ 7^4 \hspace{0.2 em}mod\hspace{0.2 em} 15=1\end{split}\]

The final step of Shor’s algorithm is to perform the inverse QFT on the source register and measure it. When we do this experiment with 1024 shots, we get a distributed measurement of states \(|0 \rangle\), \(|2 \rangle\), \(|4 \rangle\), \(|6 \rangle\) with almost equal probability.The outcome of this experiment is shown in the histogram at the end of this section. By analyzing the distribution pattern, we can figure out that the periodicity in the amplitude of \(|y \rangle\) is \(2\). Therefore, we can calculate \(r\) as \(r= \frac{2^3}{2}=4\), which we saw when we performed the modular exponentiation. Now since \(r\) is known, we can calculate the factors of \(N\) using the brute force method.

\[\begin{split}gcd\hspace{0.2 em} ⁡(7^2+1,15 )=5\\ gcd\hspace{0.2 em} ⁡(7^2-1,15 )=3\end{split}\]

We conclude that the factors of 15 are 5 and 3. Shor’s algorithm indeed works!

Source code to factorize 15

Note

Be sure to use your API token and your account name.

Step 1. Import the required modules and obtain the backend

import QuantumRingsLib
from QuantumRingsLib import QuantumRegister, AncillaRegister, ClassicalRegister, QuantumCircuit
from QuantumRingsLib import QuantumRingsProvider
from QuantumRingsLib import job_monitor
from QuantumRingsLib import JobStatus
from matplotlib import pyplot as plt
import numpy as np

provider = QuantumRingsProvider(token =<YOUR_TOKEN_HERE>, name=<YOUR_ACCOUNT_NAME_HERE>)
backend = provider.get_backend("scarlet_quantum_rings")
shots = 1024

provider.active_account()

Step 2. Define the core methods

def iqft_cct(qc, b, n):
    """
    The inverse QFT circuit

    Args:

        qc (QuantumCircuit):
                The quantum circuit

        b (QuantumRegister):
                The target register

        n (int):
                The number of qubits in the registers to use

    Returns:
        None

    """

    for i in range (n):
        for j in range (1, i+1):
            # for inverse transform, we have to use negative angles
            qc.cu1(  -math.pi / 2** ( i -j + 1 ), b[j - 1], b[i])
        # the H transform should be done after the rotations
        qc.h(b[i])
    qc.barrier()
    return

def plot_histogram (counts, title=""):
    """
    Plots the histogram of the counts

    Args:

        counts (dict):
            The dictionary containing the counts of states

        titles (str):
            A title for the graph.

    Returns:
        None

    """
    fig, ax = plt.subplots(figsize =(10, 7))
    plt.xlabel("States")
    plt.ylabel("Counts")
    mylist = [key for key, val in counts.items() for _ in range(val)]

    unique, inverse = np.unique(mylist, return_inverse=True)
    bin_counts = np.bincount(inverse)

    plt.bar(unique, bin_counts)

    maxFreq = max(counts.values())
    plt.ylim(ymax=np.ceil(maxFreq / 10) * 10 if maxFreq % 10 else maxFreq + 10)
    # Show plot
    plt.title(title)
    plt.show()
    return

Step 3. Perform the algorithm

# Shor’s algorithm to factorize 15 using 7^x mod 15.
numberofqubits = 7
shots = 1024

InitializeDebugger()

q = QuantumRegister(numberofqubits , 'q')
c = ClassicalRegister(4 , 'c')
qc = QuantumCircuit(q, c)

# Initialize source and target registers
qc.h(0)
qc.h(1)
qc.h(2)
qc.x(6)
qc.barrier()

# Modular exponentiation 7^x mod 15
qc.cx(q[2],q[4] )
qc.cx(q[2],q[5] )
qc.cx(q[6],q[4] )
qc.ccx(q[1],q[5],q[3] )
qc.cx(q[3],q[5] )
qc.ccx(q[1],q[4],q[6] )
qc.cx(q[6],q[4] ) #
qc.barrier()

# IQFT. Refer to implementation from earlier examples
iqft_cct (qc, q, 3)

# Measure
qc.measure(q[0], c[0])
qc.measure(q[1], c[1])
qc.measure(q[2], c[2])

# Execute the circuit
job = backend.run(qc, shots)
job_monitor(job)
result = job.result()
counts = result.get_counts()

#clean up
del q, c, qc
del result
del job

#visualize
plot_histogram(counts)

A plot of the execution results is shown below. Compare this with the calculated values.

Footnotes