[qiskit] Shore's Algorithm

Shore’s algorithm

I will be using qiskit to study quantum algorithms in my own way. Since this is a record of my personal study, I may have left out a lot of explanations.

I am following the qiskit website.

https://qiskit.org/textbook/ja/ch-algorithms/quantum-phase-estimation.html

I would like to proceed to the Shore algorithm, which has the potential to break the RSA cipher, where quantum computers are expected to have the biggest impact. The RSA currently used in web systems is easy to calculate the product of two primes, but the inverse factorization takes $O(N)$ of computation, and nowadays, at 1000 bits (about 300 digits), it takes as much time as the history of the universe, even with a supercomputer, making it practically impossible This means that it is practically impossible.

github

The file in jupyter notebook format is here

google colaboratory

If you want to run it on google colaboratory here

Author’s environment

!sw_vers

ProductName: Mac OS X
ProductVersion: 10.14.6
BuildVersion: 18G103

Python -V

Python 3.8.5

Import the basic libraries and check their versions.

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

import qiskit
import json

import matplotlib.pyplot as plt
import numpy as np
from qiskit import QuantumCircuit, Aer, execute
from qiskit.visualization import plot_histogram
from math import gcd
from numpy.random import randint
from tabulate import tabulate
from fractions import Fraction

from qiskit import IBMQ, Aer, transpile, assemble
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister

from qiskit.visualization import plot_histogram
from qiskit_textbook.tools import array_to_latex

dict(qiskit.__qiskit_version__)

{'qiskit-terra': '0.17.4',
 'qiskit-aer': '0.8.2',
 'qiskit-ignis': '0.6.0',
 'qiskit-ibmq-provider': '0.13.1',
 'qiskit-aqua': '0.9.1',
 'qiskit': '0.26.2',
 'qiskit-nature': None,
 'qiskit-finance': None,
 'qiskit-optimization': None,
 'qiskit-machine-learning': None}

Review of RSA cryptography.

Preparation

First, I would like to briefly review the theoretical part of RSA cryptography, which is the reason why quantum computers are attracting attention. When I was in college, I was a science major and studied mathematics such as group theory, but now I am just an IT engineer, so I cannot guarantee the accuracy. I will skip the explanation of the related theorems, etc.

Let $p,q$ be prime numbers, $n=pq$, $\varphi(n)=(p-1)(q-1)$ and Euler’s $\varphi$ function.

Let $a$ be an integer prime to $n$ and each other, then from Fermat’s Little Theorem

$$ a^{\varphi(n)}\equiv 1\left(\text{mod} n\right) $$

is established.

We also calculate $d$ so that $e$ and $de\equiv 1\left(\text{mod} \varphi(n)\right)$$ are satisfied. $d$ can be calculated by using some integer $k$.

$$ d e=k \cdot \varphi(n)+1 $$

and multiply by.

The receiver of the message publishes the public key $(n, e)$ and keeps the private key $d$ private.

Sending a message

Encrypt a message $m$ with the message sender as $c \equiv m^e(\bmod n)$ and send $c$.

Receive message.

The recipient of $c$ computes $m^{ed}$ as

$$ m^{ e d}=m^{k \varphi(n)+1} \equiv m \cdot (m^{\varphi(n)})^{k} \equiv m \quad \text{mod} n $$

and $m$ is obtained. We use Fermat’s Little Theorem along the way.

How to Decrypt

To decrypt the RSA cipher, we need to

find $p$ and $q$ from $n$.
find $m$ from $c \equiv m^e(\bmod n)$ ($c, e, n$ are known)

Shore’s algorithm will be an attempt to break the RSA cipher through factorization of 1.

Prime factorization with rank

Prepare $n$, the target of prime factorization, and $a$, which are prime to each other.

$$ a^r \equiv 1 \operatorname{mod} n $$

Find an even number $r>1$ such that $r=2r’$$. Using that $r=2r’$, we get

$$ a^{2r’} - 1 = \left(a^{r’} + 1\right)\left(a^{r’} - 1\right) \equiv 0 \operatorname{mod} n $$

So, $\left(a^{r’} + 1\right)\left(a^{r’} - 1\right) $ contains all the prime factors of $n$. Therefore, by finding the greatest common divisor of each value and the original number $n$, there is a possibility that we can do a prime factorization. However, if you are unlucky and one of them may contain $n$, you will have to choose another $a$ and recalculate $r$ again.

$$ p=\operatorname{gcd}\left(a^{r’}+1, n\right), q=\operatorname{gcd}\left(a^{r’}-1, n\right) $$

Implemented in python

Let’s see if we can implement the above algorithm in python using simple $N$ and $a$.

We will use $N=52$ and $X=17$ to find the number of positions.

n = 52
x = 17

ret_list = []
max_iter = 50
for a in range(max_iter):
  ret_list.append(x**a % n)

plt.grid()
plt.plot(range(max_iter), ret_list)
plt.show()

ret_list.index(1, 10)

r = [i for i, v in enumerate(ret_list) if v == 1]
print('rank r = {}'.format(r[1] - r[0]))

Rank r = 6

Since the rank happens to be 6, which is an even number, we find the greatest common divisor of 3 and 52.

import math

r = 6

math.gcd(int(x**(r/2) + 1), n)

math.gcd(int(x**(r/2) - 1), n)

We now have 26 and 4, which are 52 prime factors. If we recursively factorize 26, then 4, and so on, we can finally factorize 52.

Find the rank r

Define u_s

The main part of the above algorithm is to find the rank $r$. It seems that no algorithm has been found to find this problem in polynomial time on a classical computer. By using quantum phase estimation in quantum computers, we can find this $r$.

The core of Shore’s algorithm is to use unitary operators like the following.

$$ U_{a, N}|y\rangle = |a y \bmod N\rangle \cdots (1) $$

This unitary operator is the same as

$$ U_{a, N}|1\rangle = |a \bmod N\rangle \\ U_{a, N}^2|1\rangle = | a^2\bmod N\rangle \\ U_{a, N}^3|1\rangle = | a^3 \bmod N\rangle \\ \vdots \\ U_{a, N}^r|1\rangle = | a^r \bmod N\rangle = | 1 \bmod N\rangle \\ $$

and so on in a circular fashion. Using this circularity, by summing, we get a new vector $|u_0\rangle$.

$$ \left|u_{0}\right\rangle=\frac{1}{\sqrt{r}}\sum_{k=0}^{r-1}\left|a^{k}\bmod N\right\rangle $$

is the eigenvector of the unitary operator $U$.

$$ U\left|u_{0}\right\rangle=\frac{1}{\sqrt{r}}\sum_{k=0}^{r-1}U\left|a^{k}\bmod N\right\rangle = |u_0\rangle $$

Next, consider the following vector that mimics this $|u_0\rangle$.

$$ \left|u_{1}\right\rangle=\frac{1}{\sqrt{r}}\sum_{k=0}^{r-1} e^{-\frac{2 \pi i k}{r}}\left|a^{k}\bmod N\right\rangle $$

Acting on this vector with $U$, we get

$$ \begin{aligned} &U\left|u_{1}\right\rangle=\frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} e^{-\frac{2 \pi i}{r} k} U\left|a^{k} \operatorname{mod} N\right\rangle\\ &=\frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} e^{-\frac{2 \pi \cdot i}{r} \cdot k}\left|a^{k+1} \bmod N\right\rangle\\ &=\frac{1}{\sqrt{r}} e^{\frac{2 \pi i}{r}} \sum_{k=0}^{r-1} e^{-\frac{2 \pi i}{r}(k+1)}\left|a^{k+1} \bmod N\right\rangle \\ &=e^{\frac{2 \pi i}{r} } |u_1\rangle \end{aligned} $$

This is also an eigenvector of $U$ and the eigenvalue is $\displaystyle e^{\frac{2 \pi i}{r}}$.

Generalizing, we get

$$ \left|u_{s}\right\rangle=\frac{1}{\sqrt{r}}\sum_{k=0}^{r-1} e^{-\frac{2 x i k}{r}s}\left|a^{k}\bmod N\right\rangle $$

and $|u_s\rangle$ can be defined, and acting on this with the unitary operator

$$ \begin{aligned} &U\left|u_{s}\right\rangle=\frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} e^{-\frac{2 \pi i}{r} ks} U\left|a^{k} \operatorname{mod} N\right\rangle\\ &=\frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} e^{-\frac{2 \pi \cdot i }{r} ks}\left|a^{k+1} \bmod N\right\rangle\\ &=\frac{1}{\sqrt{r}} e^{\frac{2 \pi i}{r} s} \sum_{k=0}^{r-1} e^{-\frac{2 \pi i}{r}(k+1)s}\left|a^{k+1} \bmod N\right\rangle \\ &=e^{\frac{2 \pi i}{r}s } |u_s\rangle \cdots (2) \end{aligned} $$

and we get the eigenvalue $e^{\frac{2 \pi i}{r}s } $.

This $u_s$ is an eigenvector of $U$, and we can use it with an eye to quantum phase estimation.

Another representation of $u_s

This $u_s$ is defined from $\left|a^{k}\bmod N\right\rangle$, but we can express $\left|a^{k}\bmod N\right\rangle$ using $u_s$ in the following way.

$$ \begin{aligned} &|u_s\rangle=\frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} e^{-\frac{2 \pi i k}{r} s} \mid a^{k} \bmod N \rangle\\ \end{aligned} $$

$\frac{1}{\sqrt{r}} e^{\frac{2 \pi i k^{\prime}}{r} s}\left|u_{s}\right\rangle$ over both sides.

$$ \begin{aligned} &\frac{1}{\sqrt{r}} e^{\frac{2 \pi i k^{\prime}}{r} s}\left|u_{s}\right\rangle=\frac{1}{r} \sum_{k=0}^{r-1} e^{-\frac{2 \pi i^{\prime} S}{r }\left(k-k^{\prime}\right)}\left|a^{k}\bmod N\right\rangle \end{aligned} $$

Sum about $S$

$$ \begin{aligned} &\frac{1}{\sqrt{r}} \sum_{s=0}^{r-1} e^{\frac{2 \pi i k^{\prime}}{r} s}\left|u_{s}\right\rangle=\frac{1}{r} \sum_{s=0}^{r-1} \sum_{k=0}^{r -1} e^{-\frac{2 \pi i s}{r}\left(k-k^{\prime}\right)}\left|a^{k}\bmod N\right\rangle \end{aligned} $$

The right side can be simplified from the completion of $k$ and $k’$.

When $k=k’$.

$$ \sum_{s=0}^{r-1} e^{-\frac{2 \pi_{n} i s}{r}\left(k-k^{\prime}\right)}=r $$

When $k\neq k’$. $$ \sum_{s=0}^{r-1} e^{-\frac{2 \pi i s}{t}\left(k-k’\right)}=0 $$

$$ \frac{1}{\sqrt{r}}\sum_{s=0}^{r-1} e^{\frac{2 \pi i k}{r}s}\left|u_{s}\right\rangle=\left|a^{k}\bmod N\right\rangle \cdots (3) $$

Important Properties of $u_s

There is one more important property about $u_s$. In (3), if $k’=0$, then

$$ \frac{1}{\sqrt{r}}\sum_{s=0}^{r-1}\left|u_{s}\right\rangle=|1\rangle \cdots (4) $$

The result is Summing over all of the eigenvectors, we get $|1\rangle$. Since $r$ is not known when solving the actual problem, we can use the superposition of the eigenvectors, $|1\rangle$, and quantum phase estimation to determine the phase of the eigenvalue $\displaystyle \boldsymbol{e}^{\frac{2 \pi i}{r} s}$ of $u_s$. frac{s}{r}$.

Quantum circuit

Now that we are ready, we will proceed with the calculation according to the following quantum circuit. Basically, it is the same as the quantum phase estimation.

We have two quantum registers, the first one is for storing the phase and the second one is for calculating the power-surplus.

The first register initializes all qubits to $|0\rangle$ and the second register initializes them to $|0\cdots 1\rangle$.

$$ |\varphi_1 \rangle = |0\cdots 0\rangle | 0\cdots 1\rangle $$

Apply adamantine gates.

Apply an adamantine gate to all qubits in the first register.

$$ \begin{aligned} \left|\varphi_{2}\right\rangle &=\frac{1}{\sqrt{2^{n}}}(|0\rangle+|1\rangle)^{\otimes n} \otimes|0 \cdots 1\rangle \\ &=\frac{1}{\sqrt{2^{n}}}\sum_{k=0}^{2^{n}-1}|k\rangle \otimes |0 \cdots 1\rangle \end{aligned} $$

Power Multiplication (Controlled Unitary Arithmetic)

Perform the control unitary operation. Here, we use (3) and denote it by $|u_s\rangle$.

$$ \begin{aligned} &\left|\varphi_{3}\right\rangle=\frac{1}{\sqrt{2^{n}}} \sum_{k=0}^{2^{n}-1}|k\rangle\otimes\left|a^{k} \bmod n\right\rangle\\ &=\frac{1}{\sqrt{2^{n}}} \sum_{k=0}^{2^{n}-1} |k\rangle \otimes\frac{1}{\sqrt{r}} \sum_{s=0}^{r-1} e^{\frac{2 \pi i k}{r}} s\left|u_{s}\right\rangle\\ &=\frac{1}{\sqrt{r}} \sum_{s=0}^{r-1}\frac{1}{\sqrt{2^{n}}} \sum_{k=0}^{2^{n}-1} e^{\frac{2 \pi i k}{r} s}|k\rangle\otimes\left|u_{s}\right\rangle \end{aligned} $$

Inverse Quantum Fourier Transform

Next, we perform an inverse quantum Fourier transform on the first register to obtain

$$ \begin{aligned} &\frac{1}{\sqrt{2^{n}}}\sum_{k=0}^{2^{n}-1} e^{\frac{2 \pi i k}{r} s}|k\rangle \\ \rightarrow &\frac{1}{{2^{n}}} \sum_{k=0}^{2^{n}-1} \sum_{t=0}^{r-1} e^{\frac{2 \pi i k}{r} s} e^{-2 \pi i k t}|t\rangle \\ =&\frac{1}{{2^{n}}} sum_{k=0}^{2^{n}-1} sum_{t=0}^{r-1} e^{2\pi i k\left(\frac{s}{r}-t\right)}|t\rangle \end{aligned} $$

Here, measuring the first register gives the highest probability that the phase of $t=\frac{r}{s}$ will be measured, and the phase of $\frac{s}{k}$ will be stored in each qubit of the first register.

$$ \left|\varphi_{4}\right\rangle=\frac{1}{\sqrt{r}} \sum_{s=0}^{r-1}\left|\frac{s}{r}\right\rangle\otimes\left|u_{s}\right\rangle $$

However, since this state is observed in superposition for $s$, it is random (sampling from a uniform distribution) which $s$ will be observed, so the result must be checked. It will fail if $s=0$ or if $s,r$ is not a prime number.

Implemented with qiskit

This is mostly a copy, but I’ll go ahead and implement it while understanding the code on the qiskit website. We will find the period of $a=7, N=15$, but first we will run it in python to find the answer.

n = 15
a = 7

ret_list = [].
max_iter = 50
for k in range(max_iter):
  ret_list.append(a**k % n)

ret_list.index(1, 10)

r = [i for i, v in enumerate(ret_list) if v == 1].
print('rank r = {}'.format(r[1] - r[0]))

Rank r = 4

Function to find the power of $a$.

def c_amod15(a, power):
    if a not in [2,7,8,11,13]:
        raise ValueError("'a' must be 2,7,8,11 or 13")
    U = QuantumCircuit(4)
    for iteration in range(power):
        if a in [2,13]:
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
        if a in [7,8]:
            U.swap(2,3)
            U.swap(1,2)
            U.swap(0,1)
        if a == 11:
            U.swap(1,3)
            U.swap(0,2)
        if a in [7,11,13]:
            for q in range(4):
                U.x(q)
    U = U.to_gate()
    U.name = "%i^%i mod 15" % (a, power)
    c_U = U.control()
    return c_U

This function may look confusing at first glance, but the circuit to multiply $a$ in a quantum circuit looks like this We have prepared a register for the calculation and let it calculate there.

# Specify variables
n_count = 8 # number of counting qubits
a = 7

# Compute the inverse quantum Fourier transform　
def qft_dagger(n):
    """Apply the inverse QFT of n qubits to the first n qubits of the circuit."""
    qc = QuantumCircuit(n)
    # Don't forget the Swaps!
    for qubit in range(n//2):
        qc.swap(qubit, n-qubit-1)
    for j in range(n):
        for m in range(j):
            qc.cu1(-np.pi/float(2**(j-m)), m, j)
        qc.h(j)
    qc.name = "QFT†"
    return qc

This completes the preparation, and the above functions are executed sequentially.

# With n_counts of qubits for measurement and 4 qubits to manipulate U
# Create a quantum circuit
qc = QuantumCircuit(n_count + 4, n_count)

# Initialize the measurement qubits to the
# Initialize the measurement qubit to the # |+> state
for q in range(n_count):
    qc.h(q)

# Rightmost of the second register.
# Put the ancillary registers into the |1> state
qc.x(3+n_count)

# Manipulate control U
for q in range(n_count):
    qc.append(c_amod15(a, 2**q),
             [q] + [i+n_count for i in range(4)])

# Operate the inverse QFT
qc.append(qft_dagger(n_count), range(n_count))

# Measure the circuit
qc.measure(range(n_count), range(n_count))
qc.draw('mpl')

<qiskit.circuit.gate.Gate object at 0x13264f520>
<qiskit.circuit.gate.Gate object at 0x132316040>
<qiskit.circuit.gate.Gate object at 0x132612a00>
<qiskit.circuit.gate.Gate object at 0x132bba400>
<qiskit.circuit.gate.Gate object at 0x1336bab20>
<qiskit.circuit.gate.Gate object at 0x132626fd0>
<qiskit.circuit.gate.Gate object at 0x13369cfa0>
<qiskit.circuit.gate.Gate object at 0x133a5e160>

backend = Aer.get_backend('qasm_simulator')
results = execute(qc, backend, shots=2048).result()
counts = results.get_counts()
plot_histogram(counts)

The measurement yields four qubits. We will try to calculate which phase these correspond to.

rows, measured_phases = [], []
for output in counts:
    decimal = int(output, 2) # convert binary numbers to decimal
    phase = decimal/(2**n_count) # find eigenvalues
    measured_phases.append(phase)
    # Append these values to the table rows: rows.append(["phase"])
    rows.append(["%s(bin) = %i(dec)" % (output, decimal),
                 "%i/%i = %.2f" % (decimal, 2**n_count, phase)])
# Use tabulate to print the rows as an ASCII table: # print(tabulate(dec))
print(tabulate(rows,
               headers=["Register Output", "Phase"],
               colalign=("left", "right")))

Register Output Phase
------------------------ --------------
00000000(bin) = 0(dec) 0/256 = 0.00
11000000(bin) = 192(dec) 192/256 = 0.75
10000000(bin) = 128(dec) 128/256 = 0.50
01000000(bin) = 64(dec) 64/256 = 0.25

We will use the phase obtained here to find the form of the fraction using the continued fraction algorithm, and finally find $r$. A continued fraction is a method of approximating a real number by using integers, as shown below.

$$ a_{0}+\frac{1}{a_{1}+\frac{1}{a_{2}+\frac{1}{\cdots+\frac{1}{a_{m}}}}} $$

python implements continuous fraction expansion, e.g.

Fraction(0.25)

Fraction(1, 4)

Fraction(0.875)

Fraction(7, 8)

Fraction(7, 8). However, for values with infinite decimals, the

Fraction(0.3333333)

Fraction(6004798902680711, 18014398509481984)

which are hard to understand, we can put an upper bound on the denominator value.

Fraction(0.3333333).limit_denominator(100000000)

Fraction(3333333, 10000000)

Fraction(0.333333).limit_denominator(1000000000000)

Fraction(243429975650, 730289999979)

A maximum of 15 is sufficient here.

Fraction(0.33333333).limit_denominator(15)

Fraction(1, 3)

Now we can actually use the obtained Phase to find $r$ using the continuous fraction algorithm.

rows = [].
for phase in measured_phases:
    frac = Fraction(phase).limit_denominator(15)
    rows.append([phase, "%i/%i" % (frac.numerator, frac.denominator), frac.denominator])
# Display the ASCII table
print(tabulate(rows,
               headers=["Phase", "Fraction", "Guess for r"],
               colalign=('right','right','right')))

  Phase Fraction Guess for r
------- ---------- -------------
      0 0/1 1
   0.75 3/4 4
    0.5 1/2 2
   0.25 1/4 4

As you can see, we have obtained the correct results for the two Phases. As mentioned above, it is possible that the exact $r$ may not be obtained depending on the $s$ being randomly sampled. In that case, repeat the experiment again.

Find p and q

Once the even number $r$ is obtained, the prime numbers $p and q$ can be obtained by factorizing the following and finding the greatest common divisor.