The Shor algorithm

Quantum
Proving the Shor algorithm.

Shor’s algorithm performs integer factorization on a quantum computer, which can break many asymmetric (public/private key) cryptosystems, such as RSA or Diffie–Hellman. Many secure protocols, including HTTPS, SSH, and TLS, rely on these cryptosystems to guarantee encryption and authenticity. In mathematical terms, Shor’s solves the hidden subgroup problem for finite Abelian groups. In layman’s terms, Shor’s algorithm could expose encrypted information, such as passwords, credit cards, or other confidential items, transmitted over the Internet.

Before explaining Shor’s algorithm we’ll describe the quantum Fourier transform used in the proof and highlight some elements of number theory (the order of an integer). Thereafter you can find an explanation of how Shor’s algorithm can crack RSA (and similar) encryption, provided Shor’s algorithm is (physically) implemented. Finally, we show how the algorithm is related to the hidden subgroup problem. A snippet of Python code is added as well which highlights how one can implement the quantum Fourier transform with nothing but Numpy.

Quantum Fourier transform

The Hilbert space operator \[U_F = \frac{1}{\sqrt{2}}\sum_{j=0}^{2^n-1}\sum_{k=0}^{2^n-1} e^{-i\,2\pi kj/2^n}|k\rangle\langle j|\]

is called the quantum Fourier transform. Taking the adjoint of this expression

\[U_F^\dagger = \frac{1}{\sqrt{2}}\sum_{j=0}^{2^n-1}\sum_{k=0}^{2^n-1} e^{i\,2\pi kj/2^n}|j\rangle\langle k|\]

and multiplying leads to

\[\begin{align}U_F . U_F^\dagger &= \frac{1}{2^n}\sum_{jklm}e^{i2\pi(kj-lm)/2^n}|j\rangle\langle k | \,l \rangle\langle m|\\ &= \frac{1}{2^n}\sum_{jklm}e^{i2\pi(kj-lm)/2^n}|j\rangle \delta_{kl}\langle m| \\ &= \frac{1}{2^n}\sum_{jklm}e^{i2\pi k(j-m)/2^n}|j\rangle \langle m|. \end{align}\]

The sum

\[\sum_{k}e^{i2\pi k(j-m)/2^n} = \sum_{k}\left(e^{i2\pi (j-m)/2^n}\right)^k \]

is equal to \(2^n\) if \(j=m\) and is zero if not:

\[\sum_{k}\left(e^{i2\pi (j-m)/2^n}\right)^k = \frac{1-e^{i2\pi (j-m)}}{1- e^{i2\pi (j-m)/2^n}} = 0.\]

Hence,

\[\sum_{k}\left(e^{i2\pi (j-m)/2^n}\right)^k =2^n\, \delta_{jm}\]

and

\[U_F . U_F^\dagger = \sum_j |j\rangle\langle j| = 1.\]

thus proving that the quantum Fourier operator is a unitary operation and can be implemented as a quantum gate in a circuit.

As an example, if we take the Hilbert space \(\mathbb{C}^8\) and the state

\[|\psi\rangle = \frac{1}{2}\sum_{j=0}^7 \cos(2\pi j/8)|j\rangle\] as a superposition of \(|000\rangle,\dots |7\rangle = |111\rangle = |1\rangle\otimes|1\rangle\otimes|1\rangle\) The quantum Fourier transformation is then

\[U_F = \frac{1}{2\sqrt{2}}\sum_{j=0}^7\sum_{k=0}^7 e^{-i 2\pi kj/8}|k\rangle\langle j|\]

and the coefficients are

\[\begin{align} \hat{\psi}_k &= \sum_0^7 e^{-i 2\pi kj/8} \cos(2\pi j/8) \\ &= \frac{1}{2} \sum_0^7 \left(e^{i 2\pi (1-k)j/8} + e^{i 2\pi (1+k)j/8} \right)\\ &= 4(\delta_{k1} + \delta_{k7}). \end{align} \]

and the Fourier transformed state becomes

\[|\hat{\psi}\rangle = \frac{1}{\sqrt{2}}(|1\rangle + |7\rangle).\]

The quantum Fourier transform can be performed efficiently on a quantum computer, with a particular decomposition into a product of simpler unitary matrices. Using a simple decomposition, the discrete Fourier transform on \(2^{n}\) amplitudes can be implemented as a quantum circuit consisting of only \(O(n^2\)) Hadamard gates and controlled phase shift gates, where \(n\) is the number of qubits. This can be compared with the classical discrete Fourier transform, which takes \(O(n2^{n})\) gates (where \(n\) is the number of bits), which is exponentially more than \(O(n^{2})\). However, the quantum Fourier transform acts on a quantum state, whereas the classical Fourier transform acts on a vector, so not every task that uses the classical Fourier transform can take advantage of this exponential speedup.

Euler’s totient function

By definition, the Euler \(\phi\)-function of an integer \(n\) counts the number of positive integers less than \(n\) that are coprime to \(n\). Equivalently, \(\phi(n) = m\) where \(1\leq m<n\) and \(\text{gcd}(m,n)=1.\) Below is a drawing representing the first hundred values of the totient function.

Related to this definitions is the following theorem:

if \(a\in\mathbb{Z}_n\) then \(a^{\phi(n)}=1 \,(\text{mod} \;n).\)

The proof of this can be found in any textbook on number theory. If you compute a bit the residues in, say mod 19, you get the following symmetric visualization

and one defines the order of an integer \(b\) the smallest value \(k\) such that

\[b^k = 1 \,(\text{mod} \;n) \]

and denotes this by \(\text{ord}_nb.\) There is a large body of terrific number theory I could add here but the above definitions are sufficient towards Shor’s algorithm.

Breaking RSA

Suppose Bob picks a number \(n = pq\) where \(p\) and \(q\) are two large primes. He then picks a random \(c\) such that \(gcd(c; (p - 1)(q -1)) = 1\). He also computes \(d = c^{-1} \,(\text{mod} \;n)\). The values \((n,c)\) are made available to the public. Alice wants to send the message a to Bob. She computes \(b= a^c \,(\text{mod} \;n)\), and sends it to Bob. Bob computes \(a' = bd = a^{cd} \,(\text{mod} \;n)\) . In this way Bob gets \(a\) without anyone else finding out. The reason is as follows. The obvious method of recovering a from \(b\) is to compute \(d\). However, to compute \(d\) one needs to know \((p - 1)(q- 1)\). This in turn requires factoring \(n\), for which there is no known polynomial time algorithm. Other attacks can also be made, some more and some less clever. Unfortunately, there is no known classical polynomial algorithm for recovering the message.

There is, however, a clever approach to breaking the code via periods.

Suppose that one can somehow compute the order \(r\) of \(b\) (the message that Alice sent) in the modular group \(n\). From the above discussion we know that \(r\) divides the order of the group \((p - 1)(q - 1)\).

One can also assert that the order of \(b\) in \(n\) is the same as the order of \(a\). This is because \(b^d= a\,(\text{mod} \;n)\) and \(a^c = a\,(\text{mod} \;n)\) . So, the subgroup generated by \(a\) contains \(b\) and the subgroup generated by \(b\) contains \(a\). This can only happen if the two subgroups are equivalent. If they are equivalent then they have the same number of elements and so the order of \(a\) and \(b\) are the same. We have already computed the order of \(b\). Now we know this is also the order of \(a\).

We also know that \(c\) does not have any common factors with \((p - 1)(q- 1)\). Also, \(r\) divides \((p-1)(q-1)\). This implies that \(c\) and \(r\) have no common factors. We can do the following computation:

\[c = c'\,(\text{mod} \;r) \]

hence \[c'd' = 1\,(\text{mod} \;r) \]

meaning that

\[cd' = 1\,(\text{mod} \;r) \]

and there exists an \(m\) such that \(cd' = 1 +mr\). The value of \(d'\) can easily be computed via the Euclidean algorithm. Finally

\[\begin{align}b^{d'} &= a^{cd'}\\ &= a^{1+mr} \\ &= a\,(a^r)^m\\ &= a. \end{align}\] In this way the value of \(a\) has been recovered and the encrypted message has been broken. This of course depends on being able to find the period of \(b\). However, there is no known polynomial time algorithm for doing this and this is where quantum computing and Shor’s algorithm comes in.

Shor’s algorithm

First, prepare two registers of qubits to work with and call them input and output registers. The output needs to have space for \(n_0=\log n\) qubits. The input register needs space for \(n = 2n_0\) qubits. The initial state of the system is

\[|\psi_0\rangle =|0\rangle\otimes |0\rangle. \]

Apply the quantum Fourier transform \(U_F\) to the first register and since all values are zero the exponent is always 1:

\[\begin{array}\,|\psi_1\rangle &= U_F\,|\psi_0\rangle \otimes |0\rangle \\ &= \left(\frac{1}{2^{n/2}}\sum_x|x\rangle\right)\otimes|0\rangle \end{array}\]

Define the function \(f(x) = b^x \,(\text{mod} \;n)\) and create a quantum gate implementing this function. That is, cook a unitary transformation \(\hat{f}\) such that \[\hat{f}\,|x\rangle\otimes |0\rangle = |x\rangle\otimes |f(x)\rangle\]

and apply it to our Fourier transformed state:

\[\begin{array}\,|\psi_2\rangle &= \hat{f}(x)\,|\psi_1\rangle \\ &= \frac{1}{2^{n/2}}\sum_x|x\rangle\otimes|f(x)\rangle. \end{array}\]

Now, measure the output register and get some value \(f_0(x_r)\) where \(x_r = x_0 + kr.\) The value \(x_0\) is the smallest value for which \(f(x_0) = f(x_r).\) Through the measurement the input register is altered to

\[ \,|\psi_3\rangle = \frac{1}{\sqrt{m}}\sum_k|x_0+kr\rangle\otimes|f(x_r)\rangle \] and \(m\) is the smallest integer such that \(x_0 +mr \geq 2^n.\) From here on, forget about the output register and focus on the input one \[ \,|\psi_4\rangle = \frac{1}{\sqrt{m}}\sum_k|x_0+kr\rangle. \] Apply the quantum Fourier transform once again to obtain \[\begin{array} \,|\psi_5\rangle &= U_F\,|\psi_4\rangle \\ &= \sum_k\frac{1}{2^{n/2}\sqrt{m}}\sum_y e^{2\pi i\,(x_0+kr)y/2^n}\,|y\rangle\\ &=\sum_y \left(e^{2\pi ix_0y/2^n}\frac{1}{\sqrt{m2^n}}\sum_k e^{2\pi ikry/2^n}\right)|y\rangle. \end{array} \] The phase factor in front of the terms have modulus one and since the wave function sits in \(P_1(\mathbb{C})\) we can ommit it: \[ |\psi_6\rangle = \sum_y \left(\frac{1}{\sqrt{m2^n}}\sum_k e^{2\pi ikry/2^n}\right)|y\rangle.\]

Upon measurement we have the probability \(p(y)\) of finding \(y\) equal to \[p(y) = \frac{1}{m2^n}\left|\sum_k e^{2\pi ikry/2^n}\right|^2.\]

The exponential has it’s maximum when it’s argument is of the form \(2\pi l\) where \(l\) is some integer. The exponential will have maximums when \(y\) is close to an integer multiple of \(2^n/r\). It can be shown that for integers \(j\) if \(|\,y - j^{2n}/r\,|<1/2\) then the probability of obtaining such a \(y\) is at least 40%. We can repeat the quantum part efficiently and with near certainity obtain such a \(y\). The expression can be rewritten as

\[\left|\frac{y}{2^n}-\frac{j}{r}\right|<\frac{1}{2^{n+1}}\]

and for \(n\) large the value of \(y/2^n\) is near \(j/r.\) There is a low chance that any two random numbers would have a common factor. If \(j\) and \(r\) don’t have a common factor than r is simply the numerator. If not we can either use a classical computer to compute multiples of the numerator and and out which one of them is \(r\). Or we can run the quantum part all over again and get a new \(j=r\).

To check if we have the right \(r\) we can simply compute \(b^r \,(\text{mod} \;r)\) . If it’s equal to 1, we have our period and thus cracked the encryption.

The hidden subgroup problem

If you have a group \(G\) the normalizer \(N\) of the group consits of the elements which commute with all elements of the group. That is

\[N= \{n\in G \, | n.G = G.n\}.\]

The normalizer is often denoted by \(N \triangleleft G\) and is a step towards defining the quotient \(G/N\). A function \(f\) on the group \(G\) is hiding \(N\) if it’s constant on the cosets of \(G/N\). That is, for \(f:G\mapsto X\) one has \(f(g.N) = f(g)\) where \(X\) is an arbitrary target space. One could say that the function is gauge invariant under the normalizer.

The hidden subgroup problem is inverting the situation. Given a group \(G\) and a function which supposedly hides a subgroup \(N\) can one find the generator of the group \(N\)? The generator of a (sub)group consists of a subset of elements which through group combinations generate the subgroup.

The period-finding problem tackled by the Shor algorithm can be seen as a hidden subgroup problem. Let \(f\) be a periodic function on \(\mathbb{Z}_n\) with period \(r\) that divides \(n\). The subgroup \(H\triangleleft \mathbb{Z}_n\) generated by \(r\) is the hidden subgroup. Once a generator \(h\) for \(H\) has been found, the period \(r\) can be found by taking the GCD of \(h\) and \(n\); \(r = \text{GCD}(h,n)\).

Sample Python code

SymPy as well as many other packages have an implementation of Shor’s algorithm. ProjectQ even has rendering of the circuit diagrams (like below)

On Github and elsewhere in various Jupyter notebooks you can also find fun snippets to experiment with the quantum Fourier transform, like the one below:

import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline
 
def Hadamard(n):
    def Hn(H=np.array([[1, 1], [1, -1]], dtype=np.complex64), n=n):
        if n > 1:
            return Hn(H=np.kron(np.array([[1, 1], [1, -1]], dtype=np.complex64), H), n=n-1)
        return H 
    return Hn(n=n)
 
def QFT(t):
    Q = np.zeros(shape=(2 ** t, 2 ** t), dtype=np.complex64)
    N = 2 ** t
    for i in range(N):
        for j in range(N):
            Q[i][j] = np.exp(np.pi * 2j * ((i * j) % N) / N)
 
    return Q
 
N = 21
t = 9
H = Hadamard(t)
 
reg1 = np.zeros(shape=(2 ** t), dtype=np.complex64)
reg2 = np.ones(shape=(2 ** t), dtype=np.complex64)
reg1[0] = 1
reg1 = H.dot(reg1)
 
for i in range(2 ** t):
    reg2[i] = 2 ** i % N
 
r = reg2[0]
 
for i in range(2 ** t):
    if reg2[i] != r:
        reg1[i] = 0
 
Q = QFT(9)
reg1 = np.linalg.inv(Q).dot(reg1)
 
fig, ax = plt.subplots( nrows=1, ncols=1 )
ax.plot(abs(reg1))
[[plt]].close(fig)
plt.show()

which produces the Fourier spectrum