Shor's Algorithm for Quantum Factorisation

Motivation

In contrast to finding and multiplying of large prime numbers, no efficient classical algorithm for the factorisation of large number is known. An algorithm is called efficient if its execution time i.e. the number of elementary operations is assymtotically polynomial in the length of its input measured in bits. The best known (or at least published) classical algorithm (the quadratic sieve) needs $O\left(\exp\left(({64 \over 9})^{1/3} N^{1/3}(\ln N)^{2/3}\right)\right)$ operations for factoring a binary number of $N$ bits [12] i.e. scales exponentially with the input size.

The multiplication of large prime numbers is therefore a one-way function i.e. a function which can easily be evaluated in one direction, while its inversion is practically impossible. One-way functions play a major roll in cryptography and are essential to public key crypto-systems where the key for encoding is public and only the key for decoding remains secret.

In 1978, Rivest, Shamir and Adleman developed a cryptographic algorithm based on the one-way character of multiplying two large (typically above 100 decimal digits) prime numbers. The RSA method (named after the initials of their inventors) became the most popular public key system and is implemented in many communication programs.

While it is generally believed (although not formally proved) that efficient prime factorisation on a classical computer is impossible, an efficient algorithm for quantum computers has been proposed in 1994 by P.W. Shor [11].

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 gate can be found in [13]. For a more rigid mathematical description, please refer to [14].

Modular Exponentiation

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

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

The 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$$ (3.22)

$$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. $\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 to determine $\gcd(N,a+1)$ and $\gcd(N,a-1)$ which is defined as

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

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 quantum function $F : {\vert x,0 \rangle} \to {\vert x,f(x) \rangle}$ of the integer function $f: {\bf Z} \to {\bf Z}_{2^m}$ 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 base-vectors 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}$$ (3.24)

This can e.g. be achieved by the Hadamard transform $H$. Applying $F$ to the resulting state gives

$${\vert\psi \rangle}=F \, H \,{\vert,0 \rangle} = F \, {1 \ov... ... = {1 \over 2^n} \sum_{i=0}^{2^{2n}-1} {\vert i,f(i) \rangle}$$ (3.25)

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\rc... ...nd}\quad c_j' = \sqrt{{\left\lceil r \over q \right\rceil}}$$ (3.26)

The post-measurement state ${\vert\psi' \rangle}$ of the first register consists only of base-vectors 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} = \mathit{DFT} {\vert\psi' \rangle} = \sum_{i=0}^{q-1} \tilde{c}_i' {\vert i,k \rangle}$$ (3.27)

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

$$\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}^... ...uad {\rm if}\; \alpha \not\in {\bf Z} \\ \end{array}\right.$$ (3.29)

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 efficient classical algorithm for solving this problem using continued fractions is described in [15] and is implemented in the denominator function (appendix A.2).

Since the form of a rational number is not unique, $\lambda$ and $r$ are only determined by $a/b=\lambda / r$ if $\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.^3.3

QCL Implementation

Auxiliary Functions

The implementation of the Shor algorithm uses the following functions:

• boolean testprime(int n)
  Tests whether $n$ is a prime number ^3.4
• boolean testprimepower(int n)
  Tests whether $n$ is a prime power
• int powmod(int x,int a,int n)
  Calculates $x^a\,\,{\rm mod}\,n$
• int denominator(real x,int qmax)
  Returns the denominator $q$ of the best rational approximation $\frac{p}{q}\approx x$ with $p,q<q_\mathit{max}$

For the actual implementations of these functions, please refer to appendix A.2.

The Procedure shor

The procedure shor checks whether the integer number is suitable for quantum factorisation, and then repeats Shor's algorithm until a factor has been found.

procedure shor(int number) {
  int width=ceil(log(number,2));   // size of number in bits
  qureg reg1[2*width];             // first register
  qureg reg2[width];               // second register
  int qmax=2^width;
  int factor;                      // found factor
  int m; real c;                   // measured value
  int x;                           // base of exponentiation
  int p; int q;                    // rational approximation p/q
  int a; int b;                    // possible factors of number
  int e;                           // e=x^(q/2) mod number

  if number mod 2 == 0 { exit "number must be odd"; }
  if testprime(number) { exit "prime number"; }
  if testprimepower(number) { exit "prime power"; };

  {
    {                              // generate random base
      x=floor(random()*(number-3))+2;
    } until gcd(x,number)==1;
    print "chosen random x =",x;
    Mix(reg1);                     // Hadamard transform
    expn(x,number,reg1,reg2);      // modular exponentiation
    measure reg2;                  // measure 2nd register
    dft(reg1);                     // Fourier transform
    measure reg1,m;                // measure 2st register
    reset;                         // clear local registers
    if m==0 {                      // failed if measured 0
      print "measured zero in 1st register. trying again ...";
    } else {
      c=m*0.5^(2*width);           // fixed point form of m
      q=denominator(c,qmax);       // find rational approximation
      p=floor(q*m*c+0.5);
      print "measured",m,", approximation for",c,"is",p,"/",q;
      if q mod 2==1 and 2*q<qmax { // odd q ? try expanding p/q
        print "odd denominator, expanding by 2";
        p=2*p; q=2*q;
      }
      if q mod 2==1 {              // failed if odd q
        print "odd period. trying again ...";
      } else {
        print "possible period is",q;
        e=powmod(x,q/2,number);    // calculate candidates for
        a=(e+1) mod number;        // possible common factors
        b=(e+number-1) mod number; // with number
        print x,"^",q/2,"+ 1 mod",number,"=",a,",",
              x,"^",q/2,"- 1 mod",number,"=",b;
        factor=max(gcd(number,a),gcd(number,b));
      }
    }
  } until factor>1 and factor<number;   
  print number,"=",factor,"*",number/factor;
}

Factoring 15

15 is the smallest number that can be factorised with Shor's algorithm, as it's the product of smallest odd prime numbers $3$ and $5$. Our implementation of the modular exponentiation needs $2l+1$ qubits scratch space with $l=\lceil\mathit{ld}(15+1)\rceil=4$. The algorithm itself needs $3l$ qubits, so a total of 21 qubits must be provided.

$ qcl -b21 -i shor.qcl 
qcl> shor(15)
: chosen random x = 4 
: measured zero in 1st register. trying again ... 
: chosen random x = 11 
: measured 128 , approximation for 0.500000 is 1 / 2 
: possible period is 2 
: 11 ^ 1 + 1 mod 15 = 12 , 11 ^ 1 - 1 mod 15 = 10 
: 15 = 5 * 3

The first try failed because 0 was measured in the first register of ${\vert\psi' \rangle}$ and $\lambda / r =0$ gives no information about the period $r$.

One might argue that this is not likely to happen, since the first register has 8 qubits and $256$ possible base-vectors, however, if a number $n$ is to be factored, one might expect a period about $\sqrt{n}$ assuming that the prime factors of $n$ are of the same order of magnitude. This would lead to a period $\frac{q}{\sqrt{n}}$ after the $\mathit{DFT}$ and the probability $p=\frac{1}{\sqrt{n}}$ to accidentally pick the basevector ${\vert \rangle}$, would be $p=25.8\%$.

In the special case of a start value $x=4$ the period of the modular exponentiation is 2 since $4^2 \,{\rm mod}\,15 = 1$, consequently the Fourier spectrum shows 2 peaks at ${\vert \rangle}$ and ${\vert 128 \rangle}$ and $p=1/2$.

The second try also had the same probability of failure since $11^2 \,{\rm mod}\,15=1$, but this time, the measurement picked the second peak in the spectrum at ${\vert 128 \rangle}$. With $128/2^8=1/2=\lambda / r$, the period $r=2$ was correctly identified and the factors $\gcd(11^{2/2} \pm 1\,,\,15)=\{3, 5\}$ to 15 have been found.

Footnotes

... desired.^3.3

If the supposed period $r'=b$ derived form the rational approximation $a/b\approx c\,2^{-2m}$ is odd or $\gcd(x^{r'/2}\pm 1,N)=1$, then one could try to expand $a/b$ by some integer factor $k$ in order to guess the actual period $r=kb$.

... number^3.4

Since both testfunctions are not part of the algorithm itself, short but inefficient implementations with $O(\sqrt{n})$ have been used