The Algorithm

This section describes Shor's algorithm from a functional point of view which means that it doesn't deal with the implementation for a specific hardware architecture. A detailed implementation for the Cirac-Zoller device can be found in [8]. For a more rigid mathematical description, please refer to [9].

Modular Exponentiation

Let $N=n_1 n_2$ with the greatest common divisor ${\rm gcd}(n_1,n_2)=1$ be the number to be factorised, $x$ a randomly selected number relatively prime to $N$ (i.e. ${\rm gcd}(x,N)=1)$ and $F_N$ the following function with the period $r$:

$$F_N(k)=x^k \,{\rm mod}\,N, \; F_N(k+r)=F_N(k), \; x^r \equiv 1 \,{\rm mod}\,N$$

The function $F_N$ performs a modular exponentiation, its period $r$ is the order of $x \,{\rm mod}\,N$. If $r$ is even, we can define a $y=x^{r/2}$, which satisfies the condition $y^2 \equiv 1 \,{\rm mod}\,N$ and therefore is the solution of one of the following systems of equations:

$$y_1 \equiv 1 \,{\rm mod}\,n_1 \equiv 1 \,{\rm mod}\,n_2$$

$$y_2 \equiv -1 \,{\rm mod}\,n_1 \equiv -1 \,{\rm mod}\,n_2$$

$$y_3 \equiv 1 \,{\rm mod}\,n_1 \equiv -1 \,{\rm mod}\,n_2$$

$$y_4 \equiv -1 \,{\rm mod}\,n_1 \equiv 1 \,{\rm mod}\,n_2$$

The first two systems have the trivial solutions $y_1=1$ and $y_2=-1$ which don't differ from those of the quadratic equation $y^2=1$ in ${\bf Z}$ or a Galois field ${\rm GF}(p)$ (i.e. ${\bf Z}_p$ with prime $p$). The last two systems have the non-trivial solutions $y_3=a$, $y_4=-a$, as postulated by the Chinese remainder theorem stating that a system of $k$ simultaneous congruences (i.e. a system of equations of the form $y \equiv a_i \,{\rm mod}\,m_i$) with coprime moduli $m_1, \ldots , m_k$ (i.e. ${\rm gcd}(m_i,m_j)=1$ for all $i \neq j$) has a unique solution $y$ with $0 \leq x < m_1 m_2 \ldots m_k$.

Finding a Factor

If $r$ is even and $y = \pm a$ with $a \neq 1$ and $a \neq N-1$, then $(a+1)$ or $(a-1)$ must have a common divisor with $N$ because $a^2 \equiv 1 \,{\rm mod}\,N$ which means that $a^2 = c N+1$ with $c \in {\bf N}$ and therefore $a^2-1 = (a+1)(a-1)= c N$. A factor of $N$ can then be found by using Euclid's algorithm for determing ${\rm gcd}(N,a+1)$ and ${\rm gcd}(N,a-1)$ which is defined as

$${\rm gcd}(a,b)=\left\{\begin{array}{c l} b & {\rm if}\; a \,... ...od}\,b \neq 0 \\ \end{array}\right. \quad{\rm with}\quad a>b$$

It can be shown that a random $x$ matches the above mentioned conditions with a probability $p > {1 \over 2}$ if $N$ is not of the form $N=p^\alpha$ or $N=2 p^\alpha$. Since there are efficient classical algorithms to factorise pure prime powers (and of course to recognise a factor of 2), an efficient probabilistic algorithm for factorisation can be found if the period $r$ of the modular exponentiation can be determined in polynomial time.

Period of a Sequence

Let $F$ be an operator of the form $F {\vert x,0 \rangle} \to {\vert x,f(x) \rangle}$ and $f: {\bf Z} \to {\bf Z}_{2^m}$ a function with the unknown period $r<2^n$.

To determine $r$, we need two registers, with the sizes of $2n$ and $m$ qubits, which should be reset to ${\vert,0 \rangle}$.

As a first step we produce a homogenous superposition of all basevectors in the first register by applying an operator $U$ with

$$U {\vert,0 \rangle} = \sum_{i=0}^{2^{2n}-1} c_i {\vert i,0 \rangle} \quad{\rm with}\quad \vert c_i\vert={1 \over 2^n}$$

This can e.g. be achieved by transforming each qubit with the ${\it U2}({\pi \over 2})$ operator (see section 2.1.2). Another possibility is the discrete fast Fourier transform (${\it FFT}$) which is defined for $2n$ qubits as

$${\it FFT}{\vert i \rangle} = {1 \over 2^n} \sum{j=0}^{2^{2n}-... ...\rm e}^{{2 \pi {\rm i} \over 2^{2n}} \, i j} {\vert j \rangle}$$

Applying $F$ to the resulting state gives

$${\vert\psi \rangle}=F \cdot {\it FFT}{\vert,0 \rangle} = F \... ... = {1 \over 2^n} \sum_{i=0}^{2^{2n}-1} {\vert i,f(i) \rangle}$$

A measurement of the second register with the result $k=f(s)$ with $s<r$ reduces the state to

$${\vert\psi' \rangle}=\sum_{j=0}^{{\left\lceil q/r \right\rcei... ...and}\quad c_j' = \sqrt{{\left\lceil r \over q \right\rceil}}$$

The post-measurement state ${\vert\psi' \rangle}$ of the first register consists only of basevectors of the form ${\vert rj+s \rangle}$ since $f(rj+s)=f(s)$ for all $j$. It therefore has a discrete, homogenous spectrum.

It is not possible to directly extract the period $r$ or a multiple of it by measurement of the first register because of the random offset $s$. The result of a Fourier transform, however, is invariant (except for phase factors which don't effect the probability spectrum) to offsets of a periodic distribution.

$${\vert\tilde{\psi}' \rangle} = {\it FFT}{\vert\psi' \rangle} = \sum_{i=0}^{q-1} \tilde{c}_i' {\vert i,k \rangle}$$

$$\tilde{c}_i'= {\sqrt{r} \over q} \sum_{j=0}^{p-1} \exp\left(... ...m_{j=0}^{p-1} \exp\left(2 \pi {\rm i}\,{ijr \over q} \right)$$

$$\quad{\rm with}\quad \phi_i=2 \pi {\rm i}\,{is \over q} \quad{\rm and}\quad p={\left\lceil q \over r \right\rceil}$$

If $q=2^{2n}$ is a multiple of $r$ then $\tilde{c}_i'={\rm e}^{\phi_i} / \sqrt{r}$ if $i$ is a multiple of $q/r$ and $0$ otherwise. But even if $r$ is not a power of $2$, the spectrum of ${\vert\tilde{\psi}' \rangle}$ shows distinct peaks with a period of $q/r$ because

$$\lim_{n \to \infty} {1 \over n} \sum_{k=0}^{n-1} {\rm e}^{2... ...quad {\rm if}\; \alpha \not\in {\bf Z} \\ \end{array}\right.$$

This is also the reason why we use a first register of $2n$ qubits when $r<2^n$ because it guarantees at least $2^n$ elements in the above sum and thus a peak width of order $O(1)$.

If we now measure the first register we will get a value $c$ close to $\lambda q/r$ with $\lambda \in {\bf Z}_r$. This can be written as $c / q = c \cdot 2^{-2n} \approx \lambda / r$. We can think of this as finding a rational approximation $a/b$ with $a,b < 2^n$ for the fixed point binary number $c \cdot 2^{-2n}$. An an efficient classical algorithm for solving this problem using continued fractions is described in [10] and also implemented in the program shor (see section 4.3).

Since the form of a rational number is not unique, $\lambda$ and $r$ are only determined by $a/b=\lambda / r$ if ${\rm gcd}(\lambda,r)=1$. The probability that $\lambda$ and $r$ are coprime is greater then $1/ {\rm ln} r$, so only $O(n)$ tries are necessary for a constant probability of success as close at 1 as desired.