Shor's algorithm

Given the high error rates of contemporary quantum computers and too few qubits to use quantum error correction, laboratory demonstrations obtain correct results only in a fraction of attempts.

In 2001, Shor's algorithm was demonstrated by a group at IBM, who factored ${\displaystyle 15}$ into ${\displaystyle 3\times 5}$, using an NMR implementation of a quantum computer with ${\displaystyle 7}$ qubits.[11] After IBM's implementation, two independent groups implemented Shor's algorithm using photonic qubits, emphasizing that multi-qubit entanglement was observed when running the Shor's algorithm circuits.[12][13] In 2012, the factorization of ${\displaystyle 15}$ was performed with solid-state qubits.[14] Later, in 2012, the factorization of ${\displaystyle 21}$ was achieved.[15] In 2019, an attempt was made to factor the number ${\displaystyle 35}$ using Shor's algorithm on an IBM Q System One, but the algorithm failed because of accumulating errors.[16] Though larger numbers have been factored by quantum computers using other algorithms,[17] these algorithms are similar to classical brute-force checking of factors, so unlike Shor's algorithm, they are not expected to ever perform better than classical factoring algorithms.[18]

Theoretical analyses of Shor's algorithm assume a quantum computer free of noise and errors. However, near-term practical implementations will have to deal with such undesired phenomena (when more qubits are available, Quantum error correction can help). In 2023, Jin-Yi Cai studied the impact of noise and concluded that "Shor's Algorithm Does Not Factor Large Integers in the Presence of Noise."[5]

A complete factoring algorithm is possible if we're able to efficiently factor arbitrary ${\displaystyle N}$ into just two integers ${\displaystyle p}$ and ${\displaystyle q}$ greater than 1, since if either ${\displaystyle p}$ or ${\displaystyle q}$ are not prime then the factoring algorithm can in turn be run on those until only primes remain.

A basic observation is that, using Euclid's algorithm, we can always compute the GCD between two integers efficiently. In particular, this means we can check efficiently whether ${\displaystyle N}$ is even, in which case 2 is trivially a factor. Let us thus assume that ${\displaystyle N}$ is odd for the remainder of this discussion. Afterwards, we can use efficient classical algorithms to check if ${\displaystyle N}$ is a prime power.[19] Since whether ${\displaystyle N}$ is a prime power or not can be determined efficiently using efficient classical algorithms, in which case ${\displaystyle N}$ has been completely factored, the rest of the algorithm may assume that ${\displaystyle N}$ is not a prime power.

If those easy cases do not produce a nontrivial factor of ${\displaystyle N}$, the algorithm proceeds to handle the remaining case. We pick a random integer ${\displaystyle 2\leq a<N}$. A possible nontrivial divisor of ${\displaystyle N}$ can be found by computing ${\displaystyle \gcd(a,N)}$, which can be done classically and efficiently using the Euclidean algorithm. If this produces a nontrivial factor (meaning ${\displaystyle \gcd(a,N)\neq 1}$), the algorithm is finished, and the other nontrivial factor is ${\textstyle {\frac {N}{\gcd(a,N)}}}$. If a nontrivial factor was not identified, then that means that ${\displaystyle N}$ and the choice of ${\displaystyle a}$ are coprime. Here, the algorithm runs the quantum subroutine, which will return the order ${\displaystyle r}$ of ${\displaystyle a}$, meaning

${\displaystyle a^{r}\equiv 1{\bmod {N}}.}$

The quantum subroutine requires that ${\displaystyle a}$ and ${\displaystyle N}$ are coprime,[2] which is true since at this point in the algorithm, ${\displaystyle \gcd(a,N)}$ did not produce a nontrivial factor of ${\displaystyle N}$. It can be seen from the equivalence that ${\displaystyle N}$ divides ${\displaystyle a^{r}-1}$, written ${\displaystyle N\mid a^{r}-1}$. This can be factored using difference of squares:

$${\displaystyle N\mid (a^{r/2}-1)(a^{r/2}+1)}$$

Since we have factored the expression in this way, the algorithm doesn't work for odd ${\displaystyle r}$ (because ${\displaystyle a^{r/2}}$ must be an integer), meaning the algorithm would have to restart with a new ${\displaystyle a}$. Hereafter we can therefore assume ${\displaystyle r}$ is even. It cannot be the case that ${\displaystyle N\mid a^{r/2}-1}$, since this would imply ${\displaystyle a^{r/2}\equiv 1{\bmod {N}}}$, which would contradictorily imply that ${\textstyle {\frac {r}{2}}}$ would be the order of ${\displaystyle a}$, which was already ${\displaystyle r}$. At this point, it may or may not be the case that ${\displaystyle N\mid a^{r/2}+1}$. If it is not true that ${\displaystyle N\mid a^{r/2}+1}$, then that means we are able to find a nontrivial factor of ${\displaystyle N}$. We compute

$${\displaystyle d=\gcd(N,a^{r/2}-1)}$$

If ${\displaystyle d=1}$, then that means ${\displaystyle N\mid a^{r/2}+1}$ was true, and a nontrivial factor of ${\displaystyle N}$ cannot be achieved from ${\displaystyle a}$, and the algorithm must restart with a new ${\displaystyle a}$. Otherwise, we have found a nontrivial factor of ${\displaystyle N}$, with the other being ${\textstyle {\frac {N}{d}}}$, and the algorithm is finished. For this step, it is also equivalent to compute ${\displaystyle \gcd(N,a^{r/2}+1)}$; it will produce a nontrivial factor if ${\displaystyle \gcd(N,a^{r/2}-1)}$ is nontrivial, and will not if it's trivial (where ${\displaystyle N\mid a^{r/2}+1}$). The algorithm restated shortly follows: let ${\displaystyle N}$ be odd, and not a prime power. We want to output two nontrivial factors of ${\displaystyle N}$.

It has been shown that this will be likely to succeed after a few runs.[2] In practice, a single call to the quantum order-finding subroutine is enough to completely factor ${\displaystyle N}$ with very high probability of success if one uses a more advanced reduction.[20]

The goal of the quantum subroutine of Shor's algorithm is, given coprime integers ${\displaystyle N}$ and ${\displaystyle 1<a<N}$, to find the order ${\displaystyle r}$ of ${\displaystyle a}$ modulo ${\displaystyle N}$, which is the smallest positive integer such that ${\displaystyle a^{r}\equiv 1{\pmod {N}}}$. To achieve this, Shor's algorithm uses a quantum circuit involving two registers. The second register uses ${\displaystyle n}$ qubits, where ${\displaystyle n}$ is the smallest integer such that ${\displaystyle N\leq 2^{n}}$. The size of the first register determines how accurate of an approximation the circuit produces. It can be shown that using ${\displaystyle 2n+1}$ qubits gives sufficient accuracy to find ${\displaystyle r}$. The exact quantum circuit depends on the parameters ${\displaystyle a}$ and ${\displaystyle N}$, which define the problem.

The algorithm consists of two main steps:

1. Use quantum phase estimation with unitary ${\displaystyle U}$ representing the operation of multiplying by ${\displaystyle a}$ (modulo ${\displaystyle N}$), and input state ${\displaystyle |0\rangle ^{\otimes 2n+1}\otimes |1\rangle }$ (where the second register is ${\displaystyle |1\rangle }$ made from ${\displaystyle n}$ qubits). The eigenvalues of this ${\displaystyle U}$ encode information about the period, and ${\displaystyle |1\rangle }$ can be seen to be writable as a sum of its eigenvectors. Thanks to these properties, the quantum phase estimation stage gives as output a random integer of the form ${\displaystyle {\frac {j}{r}}2^{2n+1}}$ for random ${\displaystyle j=0,1,...,r-1}$.
2. Use the continued fractions algorithm to extract the period ${\displaystyle r}$ from the measurement outcomes obtained in the previous stage. This is a procedure to post-process (with a classical computer) the measurement data obtained from measuring the output quantum states, and retrieve the period.

The connection with quantum phase estimation was not discussed in the original formulation of Shor's algorithm,[2] but was later proposed by Kitaev.[21]

Quantum subroutine in Shor's algorithm

In general the quantum phase estimation algorithm, for any unitary ${\displaystyle U}$ and eigenstate ${\displaystyle |\psi \rangle }$ such that ${\displaystyle U|\psi \rangle =e^{2\pi i\theta }|\psi \rangle }$, sends inputs states ${\displaystyle |0\rangle |\psi \rangle }$ into output states close to ${\displaystyle |\phi \rangle |\psi \rangle }$, where ${\displaystyle \phi }$ is an integer close to ${\displaystyle 2^{2n+1}\theta }$. In other words, it sends each eigenstate ${\displaystyle |\psi _{j}\rangle }$ of ${\displaystyle U}$ into a state close to the associated eigenvalue. For the purposes of quantum order-finding, we employ this strategy using the unitary defined by the action

$${\displaystyle U|k\rangle ={\begin{cases}|ak{\pmod {N}}\rangle &0\leq k<N,\\|k\rangle &N\leq k<2^{n}.\end{cases}}}$$

The action of ${\displaystyle U}$ on states ${\displaystyle |k\rangle }$ with ${\displaystyle N\leq k<2^{n}}$ is not crucial to the functioning of the algorithm, but needs to be included to ensure the overall transformation is a well-defined quantum gate. Implementing the circuit for quantum phase estimation with ${\displaystyle U}$ requires being able to efficiently implement the gates ${\displaystyle U^{2^{j}}}$. This can be accomplished via modular exponentiation, which is the slowest part of the algorithm. The gate thus defined satisfies ${\displaystyle U^{r}=I}$, which immediately implies that its eigenvalues are the ${\displaystyle r}$-th roots of unity ${\displaystyle \omega _{r}^{k}=e^{2\pi ik/r}}$. Furthermore, each eigenvalue ${\displaystyle \omega _{r}^{k}}$ has an eigenvector of the form ${\textstyle |\psi _{j}\rangle =r^{-1/2}\sum _{k=0}^{r-1}\omega _{r}^{-kj}|a^{k}\rangle }$, and these eigenvectors are such that

$${\displaystyle {\begin{aligned}{\frac {1}{\sqrt {r}}}\sum _{j=0}^{r-1}|\psi _{j}\rangle &={\frac {1}{r}}\sum _{j=0}^{r-1}\sum _{k=0}^{r-1}\omega _{r}^{jk}|a^{k}\rangle \\&=|1\rangle +{\frac {1}{r}}\sum _{k=1}^{r-1}\left(\sum _{j=0}^{r-1}\omega _{r}^{jk}\right)|a^{k}\rangle =|1\rangle ,\end{aligned}}}$$

where the last identity follows from the geometric series formula, which implies ${\textstyle \sum _{j=0}^{r-1}\omega _{r}^{jk}=0}$.

Using quantum phase estimation on an input state ${\displaystyle |0\rangle ^{\otimes 2n+1}|\psi _{j}\rangle }$ would result in an output ${\displaystyle |\phi _{j}\rangle |\psi _{j}\rangle }$ with each ${\displaystyle \phi _{j}}$ representing a superposition of integers that approximate ${\displaystyle 2^{2n+1}j/r}$, with the most accurate measurement having a chance of ${\textstyle {\frac {4}{\pi ^{2}}}\approx 40.55}$ of being measured (which can be made arbitrarily high using extra qubits). Thus using as input ${\displaystyle |0\rangle ^{\otimes 2n+1}|1\rangle }$ instead, the output is a superposition of such states with ${\displaystyle j=0,...,r-1}$. In other words, using this input amounts to running quantum phase estimation on a superposition of eigenvectors of ${\displaystyle U}$. More explicitly, the quantum phase estimation circuit implements the transformation

$${\displaystyle |0\rangle ^{\otimes 2n}|1\rangle ={\frac {1}{\sqrt {r}}}\sum _{j=0}^{r-1}|0\rangle ^{\otimes 2n}|\psi _{j}\rangle \to {\frac {1}{\sqrt {r}}}\sum _{j=0}^{r-1}|\phi _{j}\rangle |\psi _{j}\rangle .}$$

Measuring the first register, we now have a balanced probability ${\displaystyle 1/r}$ to find each ${\displaystyle |\phi _{j}\rangle }$, each one giving an integer approximation to ${\displaystyle 2^{2n+1}j/r}$, which can be divided by ${\displaystyle 2^{2n+1}}$ to get a decimal approximation for ${\displaystyle j/r}$.

Then, we apply the continued fractions algorithm to find integers ${\textstyle b}$ and ${\textstyle c}$, where ${\textstyle {\frac {b}{c}}}$ gives the best fraction approximation for the approximation measured from the circuit, for ${\textstyle b,c<N}$ and coprime ${\textstyle b}$ and ${\textstyle c}$. The number of qubits in the first register, ${\displaystyle 2n+1}$, which determines the accuracy of the approximation, guarantees that

$${\displaystyle {\frac {b}{c}}={\frac {j}{r}}}$$

given the best approximation from the superposition of ${\textstyle |\phi _{j}\rangle }$ was measured (which can be made arbitrarily likely by using extra bits and truncating the output). However, while ${\textstyle b}$ and ${\textstyle c}$ are coprime, it may be the case that ${\textstyle j}$ and ${\textstyle r}$ are not coprime. Because of that, ${\textstyle b}$ and ${\textstyle c}$ may have lost some factors that were in ${\textstyle j}$ and ${\textstyle r}$. This can be remedied by rerunning the quantum subroutine an arbitrary number of times, to produce a list of fraction approximations

$${\displaystyle {\frac {b_{1}}{c_{1}}}_{\textstyle ,}\;{\frac {b_{2}}{c_{2}}}_{\textstyle ,}\ldots {\vphantom {\frac {b_{1}}{c_{1}}}}_{\textstyle ,}\;{\frac {b_{s}}{c_{s}}}}$$

where ${\textstyle s}$ is the number of times the algorithm was run. Each ${\textstyle c_{k}}$ will have different factors taken out of it because the circuit will (likely) have measured multiple different possible values of ${\textstyle j}$. To recover the actual ${\textstyle r}$ value, we can take the least common multiple of each ${\textstyle c_{k}}$:

$${\displaystyle \mathrm {lcm} (c_{1},c_{2},\ldots ,c_{s})}$$

The least common multiple will be the order ${\textstyle r}$ of the original integer ${\textstyle a}$ with high probability.

Phase estimation requires choosing the size of the first register to determine the accuracy of the algorithm, and for the quantum subroutine of Shor's algorithm, ${\displaystyle 2n+1}$ qubits is sufficient to guarantee that the optimal bitstring measured from phase estimation (meaning the ${\displaystyle |k\rangle }$ where ${\textstyle k/2^{2n+1}}$ is the most accurate approximation of the phase from phase estimation) will allow the actual value of ${\displaystyle r}$ to be recovered.

Each ${\displaystyle |\phi _{j}\rangle }$ before measurement in Shor's algorithm represents a superposition of integers approximating ${\displaystyle 2^{2n+1}j/r}$. Let ${\displaystyle |k\rangle }$ represent the most optimal integer in ${\displaystyle |\phi _{j}\rangle }$. The following theorem guarantees that the continued fractions algorithm will recover ${\displaystyle j/r}$ from ${\displaystyle k/2^{2{n}+1}}$:

Theorem — If ${\displaystyle j}$ and ${\displaystyle r}$ are ${\displaystyle n}$ bit integers, and

$${\displaystyle \left\vert {\frac {j}{r}}-\phi \right\vert \leq {\frac {1}{2r^{2}}}}$$

then the continued fractions algorithm run on ${\displaystyle \phi }$ will recover both ${\textstyle {\frac {j}{\gcd(j,\;r)}}}$. and ${\textstyle {\frac {r}{\gcd(j,\;r)}}}$.

[3] As ${\displaystyle k}$ is the optimal bitstring from phase estimation, ${\displaystyle k/2^{2{n}+1}}$ is accurate to ${\displaystyle j/r}$ by ${\displaystyle 2n+1}$ bits. Thus,

$${\displaystyle \left\vert {\frac {j}{r}}-{\frac {k}{2^{2n+1}}}\right\vert \leq {\frac {1}{2^{2{n}+1}}}\leq {\frac {1}{2N^{2}}}\leq {\frac {1}{2r^{2}}}}$$

which implys that the continued fractions algorithm will recover ${\displaystyle j}$ and ${\displaystyle r}$ (or with their greatest common divisor taken out).

The runtime bottleneck of Shor's algorithm is quantum modular exponentiation, which is by far slower than the quantum Fourier transform and classical pre-/post-processing. There are several approaches to constructing and optimizing circuits for modular exponentiation. The simplest and (currently) most practical approach is to mimic conventional arithmetic circuits with reversible gates, starting with ripple-carry adders. Knowing the base and the modulus of exponentiation facilitates further optimizations.[22][23] Reversible circuits typically use on the order of ${\displaystyle n^{3}}$ gates for ${\displaystyle n}$ qubits. Alternative techniques asymptotically improve gate counts by using quantum Fourier transforms, but are not competitive with fewer than 600 qubits owing to high constants.

Given a group ${\displaystyle G}$ with order ${\displaystyle p}$ and generator ${\displaystyle g\in G}$, suppose we know that ${\displaystyle x=g^{r}\in G}$, for some ${\displaystyle r\in \mathbb {Z} _{p}}$, and we wish to compute ${\displaystyle r}$, which is the discrete logarithm: ${\displaystyle r={\log _{g}}(x)}$. Consider the abelian group ${\displaystyle \mathbb {Z} _{p}\times \mathbb {Z} _{p}}$, where each factor corresponds to modular addition of values. Now, consider the function

${\displaystyle f\colon \mathbb {Z} _{p}\times \mathbb {Z} _{p}\to G\;;\;f(a,b)=g^{a}x^{-b}.}$

This gives us an abelian hidden subgroup problem, where ${\displaystyle f}$ corresponds to a group homomorphism. The kernel corresponds to the multiples of ${\displaystyle (r,1)}$. So, if we can find the kernel, we can find ${\displaystyle r}$. A quantum algorithm for solving this problem exists. This algorithm is, like the factor-finding algorithm, due to Peter Shor and both are implemented by creating a superposition through using Hadamard gates, followed by implementing ${\displaystyle f}$ as a quantum transform, followed finally by a quantum Fourier transform.[3] Due to this, the quantum algorithm for computing the discrete logarithm is also occasionally referred to as "Shor's Algorithm."

The order-finding problem can also be viewed as a hidden subgroup problem.[3] To see this, consider the group of integers under addition, and for a given ${\displaystyle a\in \mathbb {Z} }$ such that: ${\displaystyle a^{r}=1}$, the function

${\displaystyle f\colon \mathbb {Z} \to \mathbb {Z} \;;\;f(x)=a^{x},\;f(x+r)=f(x).}$

For any finite abelian group ${\displaystyle G}$, a quantum algorithm exists for solving the hidden subgroup for ${\displaystyle G}$ in polynomial time.[3]