素因数分解アルゴリズムを学習する

この実習ではショアのアルゴリズムを学習します。名前を聞いたことがある人もいるかもしれませんが、ショアのアルゴリズム[Sho99,NC00]は最も有名な量子アルゴリズムと言っても良いでしょう。ショアのアルゴリズムの元になっている量子位相推定と呼ばれる手法を学んだ後、ショアのアルゴリズムの各ステップを実例とともに紹介します。最後に、Qiskitを使用してショアのアルゴリズムを実装し、実際に素因数分解を行ってみます。

\(\newcommand{\ket}[1]{|#1\rangle}\) \(\newcommand{\modequiv}[3]{#1 \equiv #2 \pmod{#3}}\) \(\newcommand{\modnequiv}[3]{#1 \not\equiv #2 \pmod{#3}}\)

はじめに

古典計算をはるかに上回る能力を持つ量子計算の一つの例として、最も広く知られているのがショアの量子計算アルゴリズムでしょう。このアルゴリズムが考える問題は、大きな正の数を二つの素数の積に分解するというもので、問題自体はシンプルです。しかし、古典計算では素因数分解の有効なアルゴリズムが知られておらず、数が大きくなると指数関数的に計算量が増えると予想されています。 ショアのアルゴリズムを用いれば、同じ問題を多項式時間で解くことができると考えられています(一般的に、問題のサイズに対して多項式で増える計算時間で解くことができれば、それは有効なアルゴリズムだと見なされます)。

古典計算での素因数分解の難しさは、現在広く使われている鍵暗号技術の元になっています。なので、指数関数的に高速なショアのアルゴリズムが量子コンピュータで実現されれば、秘密の情報が暴かれる可能性があります。ショアのアルゴリズムが大きく注目される理由はそこにあります。

量子位相推定

まず、ショアのアルゴリズムの元になっている「量子位相推定」Quantum Phase Estimation, QPE)と呼ばれる手法について学びましょう。ショアのアルゴリズムを理解していけば、ショアのアルゴリズムの核心部分は、実はほぼQPEそのものであることが見えてくると思います。QPEの理解には「量子フーリエ変換」(Quantum Fourier Transform, QFT)の理解が欠かせませんが、QFTについては、この実習の問題7、もしくは参考文献[1]を参照してください。

QPEはとても大事な計算手法で、ショアのアルゴリズムだけでなく、いろいろな量子アルゴリズムのサブルーチンとしても使われています。

QPEが考える問題は、「あるユニタリー演算\(U\)に対して\(U\ket{\psi}=e^{2\pi i\theta}\ket{\psi}\)となる固有ベクトル\(\ket{\psi}\)が与えられるとして、その固有値\(e^{2\pi i\theta}\)の位相\(\theta\)を求めることができるか?」という問題です。

1量子ビットの位相推定

まず、下図にあるような量子回路を考えてみましょう。上側の量子ビットは\(\ket{0}\)、下側の量子ビットには\(U\)の固有ベクトル\(\ket{\psi}\)が初期状態として与えられているとします。

qpe_1qubit

この場合、量子回路の1-3の各ステップでの量子状態は、以下のように表現できます。

  • ステップ 1 : \(\frac{1}{\sqrt{2}}(\ket{0}\ket{\psi}+\ket{1}\ket{\psi})\)

  • ステップ 2 : \(\frac{1}{\sqrt{2}}(\ket{0}\ket{\psi}+\ket{1} e^{2\pi i\theta}\ket{\psi})\)

  • ステップ 3 : \(\frac{1}{2}\left[(1+e^{2\pi i\theta})\ket{0}+(1-e^{2\pi i\theta})\ket{1}\right]\ket{\psi}\)

この状態で上側の量子ビットを測定すると、\(|(1+e^{2\pi i\theta})/2|^2\)の確率で0、\(|(1-e^{2\pi i\theta})/2|^2\)の確率で1を測定するでしょう。つまり、この確率の値から位相\(\theta\)を求めることができるわけです。 しかし、\(\theta\)の値が小さい場合(\(\theta\ll1\))、ほぼ100%の確率で0を測定し、ほぼ0%の確率で1を測定することになるため、100%あるいは0%からのずれを精度良く決めるためには測定を何度も繰り返す必要が出てきます。これではあまり優位性を感じませんね。。

Hint

この制御\(U\)ゲートの部分がショアのアルゴリズムの「オラクル」に対応しています。その後ろには逆量子フーリエ変換が来ますが、1量子ビットの場合は\(H\)ゲート(\(H=H^\dagger\))です。つまり、1量子ビットの量子フーリエ変換は\(H\)ゲートそのものということですね。

本題に戻って、では少ない測定回数でも精度良く位相を決める方法は、何か考えられるでしょうか。。

\(n\)量子ビットの位相推定

そこで、上側のレジスタを\(n\)量子ビットに拡張した量子回路(下図)を考えてみましょう。

qpe_wo_iqft

それに応じて、下側のレジスタにも\(U\)を繰り返し適用することになりますが、鍵になるのは\(U\)を2の羃乗回適用するところです。それがどういう意味を持つのかを理解するために、準備として\(U^{2^x}\ket{\psi}\)が以下のように書けることを見ておきます(まあ当然と言えば当然ですが)。

\[\begin{split} \begin{aligned} U^{2^x}\ket{\psi}&=U^{2^x-1}U\ket{\psi}\\ &=U^{2^x-1}e^{2\pi i\theta}\ket{\psi}\\ &=U^{2^x-2}e^{2\pi i\theta2}\ket{\psi}\\ &=\cdots\\ &=e^{2\pi i\theta2^x}\ket{\psi} \end{aligned} \end{split}\]

この量子回路に対して、同様に1, 2, … \(n+1\)の各ステップでの量子状態を追跡していくと、以下のようになることが分かります。

  • ステップ 1 : \(\frac{1}{\sqrt{2^n}}(\ket{0}+\ket{1})^{\otimes n}\ket{\psi}\)

  • ステップ 2 : \(\frac{1}{\sqrt{2^n}}(\ket{0}+e^{2\pi i\theta2^{n-1}}\ket{1})(\ket{0}+\ket{1})^{\otimes n-1}\ket{\psi}\)

  • \(\cdots\)

  • ステップ \(n+1\) : \(\frac{1}{\sqrt{2^n}}(\ket{0}+e^{2\pi i\theta2^{n-1}}\ket{1})(\ket{0}+e^{2\pi i\theta2^{n-2}}\ket{1})\cdots(\ket{0}+e^{2\pi i\theta2^0}\ket{1})\ket{\psi}\)

ステップ \(n+1\)後の\(n\)ビットレジスタの状態をよく見ると、この状態はQFTで\(j\)\(2^n\theta\)としたものと同等であることが分かります。つまり、この\(n\)ビットレジスタに逆フーリエ変換\(\rm{QFT}^\dagger\)を適用すれば、状態\(\ket{2^n\theta}\)が得られるわけです!この状態を測定すれば、\(2^n\theta\)つまり固有値の位相\(\theta\)(を\(2^n\)倍したもの)を求めることができるというのがQPEです(下図)。

qpe

ただし、一般に\(2^n \theta\)が整数であるとは限りません。非整数値に対する逆フーリエ変換については、補足ページを参照してください。

ショアのアルゴリズム

では、そろそろ本題であるショアのアルゴリズムに入っていきましょう。ショアのアルゴリズムが考えるのは「ある正の合成数\(N\)を、自明ではない素数の積\(N=qp\)に分解する」という問題です。

まず、整数の剰余についての表記法をおさらいしておきます。以下のような整数\(x\)の並びを考えたとき、例えば3で割った余りを\(y\)とすると

x

0

1

2

3

4

5

6

7

8

9

y

0

1

2

0

1

2

0

1

2

0

ですね。この時、\(\modequiv{x}{y}{3}\)と書くものとします(\(k\)を0以上の整数とすれば、\(x=3k+y\)と書くこともできます)。

ショアのアルゴリズムの流れを書くと、以下のようなフローチャートになります。黒字の部分は古典計算で実行し、青地の部分を量子コンピュータで実行することになります。「アルゴリズムの一部でしか量子計算を使わないのはなぜ?」と思うかもしれませんが、この青地の部分を古典計算で実行するのが難しいというのがその理由です。つまり、古典計算で十分なところは古典計算にまかして、古典計算では難しい部分を量子計算で実行するというのがその発想の大元にあります。なぜ青地の部分が古典計算では難しいのかは追々明らかになります。

shor_flow

素因数分解の例

簡単な例として、\(N=15\)の素因数分解をこのアルゴリズムに沿って考えてみましょう。

例えば、15に素な数として仮に\(a=7\)を選んだとします。そこで\(7^x\)を15で割った余りを\(y\)とすると

x

0

1

2

3

4

5

6

\(\cdots\)

y

1

7

4

13

1

7

4

\(\cdots\)

のようになります。つまり、\(\modequiv{7^r}{1}{15}\)を満たす最小の非自明な\(r\)は4になることが分かります。 \(r=4\)は偶数なので、\(\modequiv{x}{7^{4/2}}{15}\)が定義でき、\(x=4\)です。\(x+1 = \modnequiv{5}{0}{15}\)なので、

\[ \{p,q\}=\{\gcd(5,15), \gcd(3,15)\}=\{5,3\} \]

となって、\(15=5\times3\)が得られました!

量子回路

では次に、\(N=15\)の素因数分解を実装する量子回路を考えていきましょう。いきなり答えを書いてしまうようですが、回路自体は以下のような構成をしています。

shor

上段にある4個の量子ビットが測定用のレジスタ、下段の4個の量子ビットが作業用のレジスタに対応します。それぞれのレジスタが4つずつなのは、15が4ビット(\(n=4\))で表現できるからです(15の2進数表記 = \(1111_2\))。状態は全て\(\ket{0}\)に初期化されているものとして、測定用ビットの状態を\(\ket{x}\)、作業用ビットの状態を\(\ket{w}\)とします。 \(U_f\)は以下のようなオラクル

shor_oracle2

で、作業用ビットの出力状態が\(\ket{w\oplus f(x)}\)になるものと理解しておきます(詳細は後で説明します)。関数\(f(x)\)\(f(x) = a^x \bmod N\)とします。

では、同様に回路のステップ 1-5ごとに量子状態を見ていきましょう。まずステップ 1で測定用量子ビットの等しい重ね合わせ状態を生成します。各計算基底は0から15までの整数で書いておくことにします。

  • ステップ 1 :\(\frac{1}{\sqrt{2^4}}\left[\sum_{j=0}^{2^4-1}\ket{j}\right]\ket{0}^{\otimes 4} = \frac{1}{4}\left[\ket{0}+\ket{1}+\cdots+\ket{15}\right]\ket{0}^{\otimes 4}\)

オラクル\(U_f\)を適用した後の状態は、オラクルの定義から以下のようになります。

  • ステップ 2 :

\[\begin{split} \begin{aligned} &\frac{1}{4}\left[\ket{0}\ket{0 \oplus (7^0 \bmod 15)}+\ket{1}\ket{0 \oplus (7^1 \bmod 15)}+\cdots+\ket{15}\ket{0 \oplus (7^{15} \bmod 15)}\right]\\ =&\frac{1}{4}\left[\ket{0}\ket{1}+\ket{1}\ket{7}+\ket{2}\ket{4}+\ket{3}\ket{13}+\ket{4}\ket{1}+\cdots+\ket{15}\ket{13}\right] \end{aligned} \end{split}\]

ステップ 2の後に作業用ビットを測定します。\(\ket{w}\)\(\ket{7^x \bmod 15}\)、つまり\(\ket{1}\), \(\ket{7}\), \(\ket{4}\), \(\ket{13}\)のどれかなので、例えば測定の結果13が得られたとします。その場合、測定用ビットの状態は

  • ステップ 3 :\(\frac{1}{2}\left[\ket{3}+\ket{7}+\ket{11}+\ket{15}\right]\)

となります。次に、測定用ビットに逆量子フーリエ変換\(\rm{QFT}^\dagger\)を適用します。逆量子フーリエ変換はある状態\(\ket{j}\)\(\ket{j} \to \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}e^{\frac{-2\pi ijk}{N}}\ket{k}\)に変換するので、

  • ステップ 4 :

\[\begin{split} \begin{aligned} &\frac{1}{2}QFT^\dagger\left[\ket{3}+\ket{7}+\ket{11}+\ket{15}\right]\\ =&\frac{1}{2}\frac1{\sqrt{2^4}}\sum_{k=0}^{2^4-1}\left[e^{\frac{-2\pi i\cdot3k}{2^4}}+e^{\frac{-2\pi i\cdot7k}{2^4}}+e^{\frac{-2\pi i\cdot11k}{2^4}}+e^{\frac{-2\pi i\cdot15k}{2^4}}\right]\ket{k}\\ =&\frac{1}{8}\left[4\ket{0}+4i\ket{4}-4\ket{8}-4i\ket{12}\right] \end{aligned} \end{split}\]

となります。ここで、状態として\(\ket{0}\), \(\ket{4}\), \(\ket{8}\), \(\ket{12}\)しか出てこないところが鍵で、量子状態の干渉を使って不要な答えの振幅を小さくしているわけです。

  • ステップ 5 :最後に測定用ビットを測定すると、0, 4, 8, 12がそれぞれ1/4の確率で得られます。

ステップ 2で\(7^x \bmod 15\)を計算しているので想像がつきますが、すでに繰り返しの兆候が現れていますね。

測定結果の解析

この測定結果の意味を考察してみましょう。ショアのアルゴリズムの回路\(n\)量子ビット位相推定の回路の類似性から、ここではこの2つが同一の働きをするものと仮定してみます(以下で補足説明します)。その場合、測定用レジスタは固有値\(e^{2\pi i\theta}\)の位相\(\theta\)\(2^4=16\)倍したものになっているはずです。つまり、例えば測定用レジスタを測定した結果が仮に4の場合、位相\(\theta\)\(\theta=4/16=0.25\)です。この位相の値は何を意味しているのでしょうか?

ショアのアルゴリズムの量子回路として、これまで\(\ket{w}=\ket{0}^{\otimes n}\)を初期状態として、\(U_f\ket{x}\ket{w}=\ket{x}\ket{w\oplus f(x)}\) \((f(x) = a^x \bmod N)\) となるオラクル\(U_f\)を考えてきました。この\(U_f\)を実装するために、以下のようなユニタリー演算子\(U\)を考えてみます。

\[\begin{split} U\ket{m} = \begin{cases} \ket{am \bmod N)} & 0 \leq m \leq N - 1 \\ \ket{m} & N \leq m \leq 2^n-1 \end{cases} \end{split}\]

このユニタリーは、

\[ U^{x}\ket{1} = U^{x-1} \ket{a \bmod N} = U^{x-2} \ket{a^2 \bmod N} = \cdots = \ket{a^x \bmod N} \]

を満たすので、\(w=0\)とした\(U_f\ket{x}\ket{0}\)\(U\)を使って実装することができます。

\[\begin{split} \begin{aligned} U_f\ket{x}\ket{0}&=\ket{x}\ket{0 \oplus (a^x \bmod N)}\\ &=\ket{x}\ket{a^x \bmod N}\\ &=U^x\ket{x}\ket{1} \end{aligned} \end{split}\]

ここで、天下り的ではありますが

\[ \ket{\psi_s} \equiv \frac{1}{\sqrt{r}}\sum_{k=0}^{r-1}e^{-2\pi isk/r}\ket{a^k \bmod N} \]

\(s\)\(0<s<r-1\)の整数)となるベクトル\(\ket{\psi_s}\)を定義すると、

\[ \frac{1}{\sqrt{r}}\sum_{s=0}^{r-1}\ket{\psi_s}=\ket{1} \]

が導けると同時に、\(\ket{\psi_s}\)\(U\)の固有ベクトルであり、固有値\(e^{2\pi is/r}\)を持つことが分かります。

\[ U\ket{\psi_s}=e^{2\pi is/r}\ket{\psi_s} \]

つまり、ショアのアルゴリズムのオラクル\(U_f\)による操作は、固有値\(e^{2\pi is/r}\)を持つ固有ベクトル\(\ket{\psi_s}\)の重ね合わせ状態\(\ket{1}\)にユニタリー\(U\)\(x\)回適用することと同等なわけです。量子位相推定の回路と比較すれば、やっていることはほぼ同じですね。その後に逆QFTを掛けるということまで考えれば、まさにQPEそのものの操作を行っていることに対応しています。

QPEで得られる位相は何だったかを思い出すと、それは\(U\ket{\psi}=e^{2\pi i\theta}\ket{\psi}\)なるユニタリー演算\(U\)と固有ベクトル\(\ket{\psi}\)に対する固有値\(e^{2\pi i\theta}\)に含まれる位相\(\theta\)でした。以上のことから、ショアのアルゴリズムから得られる位相\(\theta\)は、\(s/r\)(の整数倍)の意味を持つことも分かるでしょう。

連分数展開

以上の考察から、測定の結果得られる位相は\(\theta \approx s/r\)であることが分かりました。この結果から位数\(r\)を求めるために連分数展開という手法を使いますが、詳細は他の文献に委ねます([NC00]にもp.230に記述があります)。この手法を使うことで、\(\theta\)に最も近い分数として\(s/r\)を求めることができます。

例えば\(\theta=0.25\)の場合、\(r=4\)が得られます(頻度は小さいですが\(r=8\)が得られる可能性もあります)。ここまでできれば、あとは古典計算のみで求める素因数に分解することができますね(ここを参照)。

剰余指数化

オラクル\(U_f\)による操作\(U_f\ket{x}\ket{w}=\ket{x}\ket{w\oplus f(x)}\)とはどういうものか、もう少し考えてみましょう。\(f(x) = a^x \bmod N\)は、\(x\)の2進数表記

\[ x=(x_1x_2\cdots x_n)_2 = 2^{n-1}x_1+2^{n-2}x_2+\cdots+2^0x_n \]

を使って

\[\begin{split} \begin{aligned} f(x) & = a^x \bmod N \\ & = a^{2^{n-1}x_1+2^{n-2}x_2+\cdots+2^0x_n} \bmod N \\ & = a^{2^{n-1}x_1}a^{2^{n-2}x_2}\cdots a^{2^0x_n} \bmod N \end{aligned} \end{split}\]

と書くことができます。つまり、この関数は以下のようなユニタリー演算を考えれば実装することができます。

shor_oracle

\(n\)量子ビットQPEの回路と比較すれば、このユニタリーはQPEの\(U^{2^x}\)演算を実装しているものだと分かるでしょう。このように、第2レジスタ(上図では一番下のワイヤに繋がるレジスタ)の内容に、第1レジスタの各ビットで制御された\(a^x \bmod N\)を適用してQPEの\(U^{2^x}\)演算を実現する手法を、剰余指数化と呼びます。

アルゴリズムの実装

ここから、ショアのアルゴリズムを実装していきます。

# Tested with python 3.7.9, qiskit 0.23.5, numpy 1.20.1

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 fractions import Fraction

位数の発見

まず最初に、繰り返しの位数(周期)を発見するアルゴリズムを見てみます。

\(N\)を正の整数として、関数\(f(x) = a^x \bmod N\)の振る舞いを考えます。ショアのアルゴリズムに立ち返ってみると、 ここで\(a\)\(N\)と互いに素な\(N\)未満の正の整数で、位数\(r\)\(\modequiv{a^r}{1}{N}\)を満たす非ゼロの最小の整数でした。 以下のグラフにこの関数の例を示します。 ポイント間の線は周期性を確認するためのものです。

N = 35
a = 3

# プロットするデータを計算する
xvals = np.arange(35)
yvals = [np.mod(a**x, N) for x in xvals]

# matplotlibを使って描画
fig, ax = plt.subplots()
ax.plot(xvals, yvals, linewidth=1, linestyle='dotted', marker='x')
ax.set(xlabel='$x$', ylabel='$%i^x$ mod $%i$' % (a, N),
       title="Example of Periodic Function in Shor's Algorithm")
try: # グラフ上にrをプロット
    r = yvals[1:].index(1) + 1
    plt.annotate(text='', xy=(0,1), xytext=(r,1), arrowprops=dict(arrowstyle='<->'))
    plt.annotate(text='$r=%i$' % r, xy=(r/3,1.5))
except:
    print('Could not find period, check a < N and have no common factors.')
_images/shor_12_0.png

オラクルの実装

以下では、\(N=15\)を素因数に分解してみます。上で説明したように、\(U\ket{m}=\ket{am \bmod N}\)となるユニタリー\(U\)\(x\)回繰り返すことで、オラクル\(U_f\)を実装します。

練習問題として、\(U\ket{m}=\ket{am \bmod 15}\)を実行する関数U_amod15を以下に実装してください(U_amod15は制御ゲートですが、標的ビットのユニタリー演算に対応する部分を書いてみてください)。

def c_amod15(a, power):
    """mod 15による制御ゲート"""
    if a not in [2,4,7,8,11,13,14]:
        raise ValueError("'a' must be 2,4,7,8,11,13 or 14")

    U = QuantumCircuit(4)

    ##################
    ### EDIT BELOW ###
    ##################
    
    #U.?

    ##################
    ### EDIT ABOVE ###
    ##################

    # 以下で制御ゲートに変換
    U = U.to_gate()
    U.name = "%i^%i mod 15" % (a, power)
    c_U = U.control()
    return c_U

powerは繰り返しの回数を表します。

回路全体の実装

測定用ビットとして、8量子ビットを使います。

# 測定用ビットの数
n_count = 8

a = 7

次に、逆QFTの回路を考えます。この実習の問題7を参考、QFTの逆回路qft_dagger(n)を書いてみてください。引数の\(n\)は測定用ビットの数n_countが入ることに注意します。

def qft_dagger(n):
    """n量子ビットの逆量子フーリエ変換回路を書く"""
    qc = QuantumCircuit(n)

    ##################
    ### EDIT BELOW ###
    ##################
    
    #qc.?

    ##################
    ### EDIT ABOVE ###
    ##################

    qc.name = "QFT^dagger"
    return qc

ショアのアルゴリズムを実装する量子回路は、以上の要素を使って構築することができます。

# n_count個の測定用量子ビットと、Uを操作するための4つの作業用量子ビットで量子回路を作る
qc = QuantumCircuit(n_count+4, n_count)

# 測定用量子ビットにHゲートをかけて初期化
qc.h(list(range(n_count)))

# 最後の作業用量子ビットを|1>の状態にする
qc.x(3+n_count)

# 制御Uゲートを適用
for q in range(n_count):
    qc.append(c_amod15(a, 2**q), [q]+[i+n_count for i in range(4)])

# 逆QFTを適用
qc.append(qft_dagger(n_count), list(range(n_count)))

# 回路を測定
qc.measure(list(range(n_count)), list(range(n_count)))
qc.draw('mpl')
_images/shor_20_0.png

シミュレータで実行して、結果を確認してみます。

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

def show_distribution(answer):
    n = len(answer)
    x = [int(key,2) for key in list(answer.keys())]
    y = list(answer.values())

    fig, ax = plt.subplots()
    rect = ax.bar(x,y)

    def autolabel(rects):
        for rect in rects:
            height = rect.get_height()
            ax.annotate('{:.3f}'.format(height/sum(y)),
                        xy=(rect.get_x()+rect.get_width()/2, height),xytext=(0,0),
                        textcoords="offset points",ha='center', va='bottom')
    autolabel(rect)
    plt.ylabel('Probabilities')
    plt.show()

show_distribution(answer)
_images/shor_22_0.png

計算結果の解析

出力された結果から、位相を求めてみます。

rows, measured_phases = [], []
for output in answer:
    decimal = int(output, 2)  # 10進数に変換
    phase = decimal/(2**n_count)
    measured_phases.append(phase)
    # これらの値をテーブルの行に追加:
    rows.append(["%i" % (decimal),
                 "%i/%i = %.2f" % (decimal, 2**n_count, phase)])
# 結果を表示
print('Register Output              Phase')
print('----------------------------------')

# 回路を実装できたら、以下のコードをコメントアウトして結果を確認
#for i in range(len(rows)):
#    print('%15s %18s' % (rows[i][0],rows[i][1]))
Register Output              Phase
----------------------------------

得られた位相の情報から、連分数アルゴリズムを使用して\(s\)\(r\)を見つけることができます。Pythonの組み込みのfractions(分数)モジュールを使用して、小数をFractionオブジェクトに変換できます。

rows = []
for phase in measured_phases:
    frac = Fraction(phase).limit_denominator(15)
    rows.append([phase, "%i/%i" % (frac.numerator, frac.denominator), frac.denominator])

# 結果を表示
print('     Phase   Fraction     Guess for r')
print('-------------------------------------')

# 回路を実装できたら、以下のコードをコメントアウトして結果を確認
#for i in range(len(rows)):
#    print('%10s %10s %15s' % (rows[i][0],rows[i][1],rows[i][2]))
     Phase   Fraction     Guess for r
-------------------------------------

limit_denominatorメソッドを使って、分母が特定の値(ここでは15)を下回る分数で、最も位相の値に近いものを得ています。

測定された結果のうち、2つ(64と192)が正しい答えである\(r=4\)を与えたことが分かります。

参考文献

NC00(1,2)

Michael A Nielsen and Isaac L Chuang. Quantum Computation and Quantum Information, chapter 5.3, pages 226–234. Cambridge University Press, 2000.

Sho99

Peter W. Shor. Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer. SIAM Review, 41(2):303–332, 1999. URL: https://doi.org/10.1137/S0036144598347011, arXiv:https://doi.org/10.1137/S0036144598347011, doi:10.1137/S0036144598347011.