Shor's Algorithm for Quantum Factorization

Motivation

In contrast to finding and multiplying of large prime numbers, no efficient classical algorithm for the factorization 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 factorization 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 [15] and for a more detailed dicussion of the QCL implementation, look at [25].

Modular Exponentiation

Let $N=n_1 n_2$ with the greatest common divisor $\gcd (n_1,n_2)=1$ be the number to be factorized, $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$$ (4.19)

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$$ (4.20)

$$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}\... ...f}\; a {\rm mod} b \neq 0 \\ \end{array}\right. \quad\mathrm{with}\quad a>b$$ (4.21)

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 factorize pure prime powers (and of course to recognize a factor of 2), an efficient probabilistic algorithm for factorization can be found if the period $r$ of the modular exponentiation can be determined in polynomial time.

Period of a Function

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}^{N-1} c_i {\vert i,0 \rangle... ...t c_i\vert={1 \over \sqrt{N}} \quad\mathrm{and}\quad N=2^{2n}$$ (4.22)

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... ...gle} = {1 \over 2^n} \sum_{i=0}^{N-1} {\vert i,f(i) \rangle}$$ (4.23)

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 N/r \right\rc... ...uad c_j' = {\left\lceil N \over r \right\rceil}^{-\frac{1}{2}}$$ (4.24)

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$. This problem can be solved by performing a discrete Fourier transform (see 4.2.3)

$$\mathit{DFT}:{\vert x \rangle} \to \frac{1}{\sqrt{N}} \s... ...\mathrm{e}^{\frac{2\pi \mathrm{i}}{N} xy} {\vert y \rangle}$$ (4.25)

on the register, as the probability spectrum of the transformed state is invariant to the offset (i.e. only the phases but not the absolute value of the complex amplitudes are effected).

$${\vert\tilde{\psi}' \rangle} = \mathit{DFT} {\vert\psi' \rangle} = \sum_{i=0}^{N-1} \tilde{c}_i' {\vert i,k \rangle}$$ (4.26)

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

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

If $N=2^{2n}$ is a multiple of $r$ then $\tilde{c}_i'={\rm e}^{\phi_i} / \sqrt{r}$ if $i$ is a multiple of $N/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 $N/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.$$ (4.28)

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 N/r$ with $\lambda \in {\bf Z}_r$. This can be written as $c / N = 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 [16] and is implemented in the denominator function (appendix B.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.^4.1

Quantum Fourier Transform

For a $N$ dimensional vector ${\vert\psi \rangle}$, the discrete Fourier transform is defined as

$$\mathit{DFT}:{\vert x \rangle} \to \frac{1}{\sqrt{N}} \s... ...\mathrm{e}^{\frac{2\pi \mathrm{i}}{N} xy} {\vert y \rangle}$$ (4.29)

Since ${\vert\psi \rangle}$ is a combined state of $n$ qubits, $N$ is always a power of 2. The classical fast Fourier Transform ($\mathit{FFT}$) uses a binary decomposition of the exponent to perform the transformation in $O(n2^n)$ steps.

As suggested by Coppersmith [7], the same principle could adapted be to quantum computers by using a combination of Hadamard transformations $H$ (see 3.4.4.3) and conditional phase gates $V$ (see 3.4.4.4). The indices below indicate the qubits operated on:

$$\mathit{DFT'}= \prod_{i=1}^{n-1}\left( H_{n-i-1}(\frac{\pi... ...1} V_{n-i-1,n-j-1}(\frac{2\pi}{2^{i-j+1}}) \right)\;H_{n-1}$$ (4.30)

$\mathit{DFT'}$ iterates the qubits form the MSB to the LSB, ``splits'' the qubits with the Hadamard transformation and then conditionally applies phases according to their relative binary position ( $\mathrm{e}^\frac{2\pi\mathrm{i}}{2^{i-j+1}}$) to each already split qubit.

The base-vectors of the transformed state ${\vert\tilde{\psi}' \rangle}=\mathit{DFT'} {\vert\psi \rangle}$ are given in reverse bit order, so to get the actual $\mathit{DFT}$, the bits have to be flipped.

operator dft(qureg q) { // main operator
  const n=#q;           // set n to length of input
  int i; int j;         // declare loop counters
  for i=0 to n-1 {
    for j=0 to i-1 {    // apply conditional phase gates
      CPhase(2*pi/2^(i-j+1),q[n-i-1] & q[n-j-1]);
    }
    Mix(q[n-i-1]);      // qubit rotation
  }    
  flip(q);              // swap bit order of the output
}

Modular Arithmetic

The most difficult part in implementing Shor's algorithm is the construction of an efficient quantum function for modular exponentiation.

$$\mathtt{expn}_{a,n}(\mathbf{b},\mathbf{e}): {\vert b \rang... ...e}_{\mathbf{b}}{\vert a^b {\rm mod} n \rangle}_{\mathbf{e}}$$ (4.31)

Assuming we already have an implementation for modular addition, we could use it to construct modular multiplication and finally exponentiation since

$$ab {\rm mod} n = \sum_{i=0}^{\lceil\mathrm{ld} b\rceil} b_i\left(2^i a {\rm mod} n\right) \quad\mathrm{with}\quad b_i\in{\bf B}$$ (4.32)

$$a^b {\rm mod} n = \prod_{i=0}^{\lceil\mathrm{ld} b\rceil} \left(a^{2^i b_i} {\rm mod} n\right) \quad\mathrm{with}\quad b_i\in{\bf B}$$ (4.33)

Modular Addition

The addition modulo $n$ of a classic integer $a$ and a quantum register $\mathbf{b}$ can result in either $a+b$ or $(a-n)+b)$, depending on the particular base-vector ${\vert b \rangle}$.

While for $b<n$ the operation is revertible, this is not the case for $b\geq n$, so, if $n$ doesn't happen to be a power of $2$, besides the target resister $\mathbf{y}_s$ for the sum, we need an additional flag-qubit $\mathbf{y}_y$ to allow for a quantum function addn which is both, unitary and invariant to $\mathbf{b}$:

$$\mathtt{addn}_{a,n} : {\vert b \rangle}_\mathbf{b}{\vert ... ...}_{flag}} & \;\mbox{if} a+b \geq n \\ \end{array}\right.$$ (4.34)

The actual implementation of addn can be found in appendix B.5.

Since $\mathtt{addn}_{n-a,n}$ is a quantum function for modular subtraction and thus implements the inverse function $f^{-1}_{a,n}(b)=b-a {\rm mod} n$ to $f_{a,n}=a+b {\rm mod} n$, we can construct an overwriting version oaddn of modular addition, by using the method introduced in 3.5.2.3:

$$F': {\vert i,0 \rangle} \stackrel{U_f}{\longrightarrow} ... ...U}^\dagger _{f^{-1}}}{\longrightarrow} {\vert f(i),0 \rangle}$$ (4.35)

$\mathtt{addn}_{n-a,n}$ doesn't invert the overflow flag $\mathbf{y}_f$, so we have to switch it manually:

$${U}^\dagger _{f^{-1}}=\mathtt{addn}_{n-a,n}(\mathbf{b},\mathbf{y}_s, \mathbf{y}_f)$$ (4.36)

The original target registers $\mathbf{y}_s$ and $\mathbf{y}_f$ can now be allocated as unmanaged local scratch.

qufunct oaddn(int a,int n,qureg sum,quconst e) {
  qureg j[#sum];
  qureg f[1];
  
  addn(a,n,sum,f,j,e);             // junk -> a+b mod n
  Swap(sum,j);                     // swap junk and sum
  CNot(f,e);                       // toggle flag
  !addn(n-a,n,sum,f,j,e);          // uncompute b to zero
}

The register $\mathbf{e}$ is an enable register (see 2.2.2.6), so addn and oaddn are in fact conditional operators which only have an effect on eigenvectors of the form ${\vert x \rangle}{\vert 111\ldots \rangle}_{\mathbf{e}}$.

Modular Multiplication

Modular multiplication is merely a composition of conditional additions for each qubit of $\mathbf{b}$. The first summand can be slightly optimized, since the accumulator (prod) is still empty.

qufunct muln(int a,int n,quconst b,qureg prod,quconst e) {
  int i;
  
  for i=0 to #prod-1 {
    if bit(a,i) { CNot(prod[i],b[0] & e); }
  }
  for i=1 to #b-1 {
    oaddn(2^i*a mod n,n,prod,b[i] & e);
  }
}

As above, we can construct an overwriting version, if an implementation of the inverse function exists. This is the case if $\gcd(a,n)=1$ so $a$ and $n$ are relatively prime, because then the modular inverse $a^{-1}$ with $a^{-1}a {\rm mod} n=1$ exists. Since we intend to use the operator for the Shor algorithm which demands that $\gcd(a^k,n)=1$, this is good enough for us.

By using two conditional XOR gates defined as

$$\mathtt{cxor} : {\vert a \rangle}_\mathbf{a}{\vert b \ran... ...le}_\mathbf{e} & \;\mbox{otherwise} \\ \end{array}\right.$$ (4.37)

for swapping the registers^4.2we can implement a conditional overwriting version of muln defined as

$$\mathtt{omuln}_{[[\mathbf{e}]],a,n}{\vert b \rangle}\to{\vert ab {\rm mod} n \rangle}$$ (4.38)

qufunct omuln(int a,int n,qureg b,quconst e) {
  qureg j[#b];

  muln(a,n,b,j,e);
  !muln(invmod(a,n),n,j,b,e);
  cxor(j,b,e);
  cxor(b,j,e);
}

Modular Exponentiation

As with muln, we can construct modular exponentiation by conditionally applying omuln with the qubits of the exponents as enable string. Before we can start the iteration, the accumulator (ex) has to be initialized by 1.

qufunct expn(int a,int n,quconst b,quvoid ex) {
  int i;
  
  Not(ex[0]);                            // start with 1
  for i=0 to #b-1 {
    omuln(powmod(a,2^i,n),n,ex,b[i]);    // ex -> ex*a^2^i mod n
  }
}

Implementation

Auxiliary Functions

The implementation of the Shor algorithm uses the following functions:

• boolean testprime(int n)
  Tests whether $n$ is a prime number
• boolean testprimepower(int n)
  Tests whether $n$ is a prime power^4.3
• 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 B.2.

The Procedure shor

The procedure shor checks whether the integer number is suitable for quantum factorization, 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 1st 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 factorized with Shor's algorithm, as it's the product of the smallest odd prime numbers $3$ and $5$. Our implementation of the modular exponentiation needs $2l+1$ qubits scratch space with $l=\lceil\mathrm{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.^4.1

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$.

... registers^4.2

normally, 3 XOR operations are necessary to swap a register, but since one register is empty, 2 XORs suffice.

... power^4.3

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