ショアのアルゴリズム

この「Qiskit in Classrooms」モジュールを実行するには、以下のパッケージがインストールされた Python 環境が必要です：

• qiskit v2.1.0 以降
• qiskit-ibm-runtime v0.40.1 以降
• qiskit-aer v0.17.0 以降
• qiskit.visualization
• numpy
• pylatexenc

上記パッケージのセットアップとインストール方法については、Qiskit のインストールガイドを参照してください。 実際の量子コンピューターでジョブを実行するには、IBM Cloud アカウントのセットアップガイドの手順に従って IBM Quantum® のアカウントを作成する必要があります。

このモジュールはテスト済みで、QPU 時間を 3 秒使用しました。これはあくまで目安であり、実際の使用量は異なる場合があります。

# Added by doQumentation — required packages for this notebook
!pip install -q numpy qiskit qiskit-ibm-runtime

# Uncomment and modify this line as needed to install dependencies
#!pip install 'qiskit>=2.1.0' 'qiskit-ibm-runtime>=0.40.1' 'qiskit-aer>=0.17.0' 'numpy' 'pylatexenc'

はじめに

1990 年代初頭、古典コンピューターには解くことが難しい問題を量子コンピューターが解けるかもしれないという期待が高まっていました。才能ある計算機科学者たちが、いくつかのニッチで特殊な問題に対して量子コンピューティングの力を示すアルゴリズムを考案していましたが、量子コンピューティングの分野を確実に革命的に変える「キラーアプリ」はまだ誰も見つけていませんでした。そこに転機が訪れたのが 1994 年、Peter Shor が大きな数を素因数分解するための現在「ショアのアルゴリズム」と呼ばれる手法を考案したときです。

当時、大きな数の素因数を求めることが古典コンピューターにとって非常に困難であることは広く知られていました。実際、インターネットのセキュリティプロトコルはこの困難さに依存していました。Shor は、将来の理論的な量子コンピューターに困難なステップの一部を委ねることで、これらの因数を指数関数的に効率よく求める方法を発見しました。

このモジュールでは、ショアのアルゴリズムを探究します。まず、アルゴリズムの背景をさらに整理し、解く問題を形式化してサイバーセキュリティとの関連を説明します。次に、モジュラー数学の入門と、それを因数分解問題にどう応用するかを解説し、因数分解が「位数探索」と呼ばれる別の問題に帰着できることを示します。前のモジュールで学んだ量子フーリエ変換と量子位相推定がどのように関わるか、そしてそれらを使って位数探索問題を解く方法を説明します。

最後に、実際の量子コンピューター上でショアのアルゴリズムを実行します！ただし、このアルゴリズムが真に役立つのは、大規模なフォールトトレラント量子コンピューターが利用できるようになってからであり、それはまだ数年先のことです。ここでは、アルゴリズムの動作を確認するために小さな数を因数分解するにとどめます。

因数分解問題

因数分解問題の目標は、数 $N$ の素因数を求めることです。一部の $N$ については、これはかなり簡単です。たとえば $N$ が偶数であれば、素因数の 1 つは 2 です。$N = p^k$（$p$ は素数）のような素数の冪であれば、$N$ の $k$ 乗根を近似し、近くの素数を探すことで $p$ を比較的容易に求められます。

しかし、古典コンピューターが苦労するのは、$N$ が奇数かつ素数の冪でない場合です。これがショアのアルゴリズムが扱うケースです。このアルゴリズムは $N = pq$ を満たす 2 つの因数 $p$ と $q$ を求めます。すべての因数が素数になるまで再帰的に適用できます。次のセクションでこの問題にどう取り組むかを見ていきます。

サイバーセキュリティとの関連

大きな数の因数分解が困難であるという事実に基づいて、多くの暗号化方式が構築されており、今日広く使われている RSA もその 1 つです。RSA では、2 つの大きな素数を掛け合わせて $N = p \cdot q$ を求めることで公開鍵を作成します。誰でもこの公開鍵を使ってデータを暗号化できますが、秘密鍵である $p$ と $q$ を持つ人だけがデータを復号できます。

$N$ を簡単に因数分解できれば、誰でも $p$ と $q$ を割り出して暗号を破ることができます。しかし実際にはそれは困難です。これは非常に難しい問題として知られています。実際、RSA1024 と呼ばれる数（2 進数で 1024 桁、10 進数で 309 桁）の素因数は、1991 年に $100,000 の懸賞が掛けられていたにもかかわらず、いまだに発見されていません。

Shor の解法

1994 年、Peter Shor は量子コンピューターが古典コンピューターよりも指数関数的に効率よく大きな数を因数分解できることに気づきました。その洞察は、この因数分解問題とモジュラー算術との関係に基づいています。まずモジュラー算術の簡単な入門を行い、次にこれを使って $N$ を因数分解する方法を見ていきます。

モジュラー算術

モジュラー算術は循環的な計数システムです。通常の 0, 1, 2, … という整数での計数から始まりますが、ある周期 $N$ の後に数え直しが始まります。例を見てみましょう。周期が 5 の場合、通常なら 5 に到達するところで 0 に戻ります：

$0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, ...$

これは「モジュロ 5」の世界では 5 が 0 と等しいからです。$5\bmod 5 \ = 0$ と表します。実際、5 のすべての倍数は $0\bmod 5$ と等しくなります。

理解度チェック

以下の問題を読み、答えを考えてから、三角形をクリックして解答を確認してください。

モジュラー算術を使って次の問題を解いてください：

あなたは午前 8 時に長距離大陸横断列車の旅に出発します。列車の旅は 60 時間です。到着時刻は何時ですか？

解答：

周期は 24 です（1 日は 24 時間のため）。したがって、この問題はモジュラー算術で次のように書けます：

$(8+60)\text{mod}(24) = 20$

つまり、目的地には 20:00（午後 8 時）に到着することになります。

$\mathbb{Z}_N$ と $\mathbb{Z}_N^*$

2 つの集合 $\mathbb{Z}_N$ と $\mathbb{Z}_N^*$ を導入すると便利です。$\mathbb{Z}_N$ は「モジュロ $N$」の世界に存在する数の集合です。たとえば、モジュロ 5 で数えていたときの集合は $\mathbb{Z}_5=\{0,1,2,3,4\}$ となります。別の例として、$\mathbb{Z}_{15} = \{0,1,2,3,4,5,6,7,8,9,10,11,12,13,14\}$ です。$\mathbb{Z}_N$ の要素に対して（モジュロ $N$ の）加算と乗算を行うと、結果も $\mathbb{Z}_N$ の要素となります。これにより $\mathbb{Z}_N$ は環と呼ばれる数学的対象になります。

ショアのアルゴリズムにとって特に重要な $\mathbb{Z}_N$ の部分集合があります。それは、各要素と $N$ の最大公約数が 1 である、つまり各要素が $N$ と「互いに素」であるような $\mathbb{Z}_N$ の部分集合です。これらの数の集合にモジュラー乗算演算を加えると、群と呼ばれる別の数学的対象が形成されます。これを $\mathbb{Z}_N^*$ と呼びます。$\mathbb{Z}_N^*$（および有限群一般）の特徴として、任意の要素 $a \in \mathbb{Z}_N^*$ を選んで $a$ を繰り返し自分自身に掛けると、必ず最終的に 1 が得られます。1 を得るために $a$ を自分自身に掛ける必要がある最小の回数を $a$ の位数と呼びます。この事実は、以下の因数分解方法の議論において非常に重要になります。

理解度チェック

以下の問題を読み、答えを考えてから、三角形をクリックして解答を確認してください。

$\mathbb{Z}_{15}^*$ とは何ですか？

解答：

$\mathbb{Z}_{15}^* = \{1,2,4,7,8,11,13,14\}$

以下の数を除外しました：

$$\begin{aligned} 3: GCD(3,15)=3 \\ 5: GCD(5,15)=5 \\ 6: GCD(6,15)=3 \\ 9: GCD(9,15)=3 \\ 10: GCD(10,15)=5 \\ 12: GCD(12,15)=3 \\ \end{aligned}$$

$\mathbb{Z}_{15}^*$ の各要素の位数はいくつですか？

解答：

位数 $r$ は、各要素 $a$ に対して $a^r\text{mod}(15)=1$ を満たす最小の数です。

$$\begin{aligned} 1^1\text{mod}(15) = 1, r=1 \\ 2^4\text{mod}(15) = 1, r=4 \\ 4^2\text{mod}(15) = 1, r=2 \\ 7^4\text{mod}(15) = 1, r=4 \\ 8^4\text{mod}(15) = 1, r=4 \\ 11^2\text{mod}(15) = 1, r=2 \\ 13^4\text{mod}(15) = 1, r=4 \\ 14^2\text{mod}(15) = 1, r=2 \\ \end{aligned}$$

なお、$\mathbb{Z}_{15}^*$ の数の位数を求めることができましたが、これはより大きな $N$ に対しては一般に簡単なことではありません。これが因数分解問題の核心であり、量子コンピューターが必要な理由です。残りのノートブックを進めていく中でその理由が明らかになります。

モジュラー算術を因数分解問題に応用する

$N = pq$ を満たす因数 $p$ と $q$ を求める鍵は、次の条件を満たす別の整数 $x$ を見つけることにあります：

$x^2 \equiv 1 \bmod N$ かつ $x \not\equiv \pm 1 \bmod N.$

$x$ を求めることがどのように $p$ と $q$ の発見につながるのでしょうか？その論理を見ていきましょう。$x^2 \equiv 1 \bmod N$ より、$x^2 - 1 \equiv 0 \bmod N$ となります。つまり、$x^2 - 1$ は $N$ の倍数です。したがって、ある整数 $l$ に対して、

$x^2 - 1 = l N$

$x^2 - 1$ を因数分解すると：

$(x+1)(x-1) = l N$

最初の仮定から $x \not\equiv \pm 1 \bmod N$ なので、$N$ は $x+1$ にも $x-1$ にも割り切れません。したがって、$N$ の 2 つの因数 $p$ と $q$ は、それぞれ $x-1$ と $x+1$ を割り切る必要があります。$p$ が $x-1$ の因数で $q$ が $x+1$ の因数か、あるいはその逆かのいずれかです。そのため、$N$ と $x-1$、$x+1$ それぞれとの最大公約数（GCD）を計算すれば、因数 $p$ と $q$ が求まります。2 つの数の GCD を計算することは古典的に容易なタスクで、たとえばユークリッドのアルゴリズムを使って実行できます。

理解度チェック

以下の問題を読み、答えを考えてから、三角形をクリックして解答を確認してください。

上記の各論理ステップを理解するのは難しいかもしれないので、具体例を使って確認してみましょう。$N=15$、$x=11$ を使います。まず $x^2 \equiv 1 \text{mod}(N)$ かつ $x \not\equiv \pm 1 \bmod N$ であることを確認してください。次に各ステップを確認し、最後に $\text{GCD}(11\pm1,15)$ を計算して、それらが $15$ の因数であることを確認してください。

解答：

$11^2 = 121$、これは $15*8 + 1$ なので、$11^2\bmod 15 = 1$。$\checkmark$

$11 - 1 = 10$ は $0\bmod 15$ と等しくありません。$\checkmark$

$11 + 1 = 12$ は $0\bmod 15$ と等しくありません。$\checkmark$

次に、ある整数 $l$ に対して $(x+1)(x-1) = l N$ が成り立つことを確認します。$x$ と $N$ を代入すると $(12)(10) = l 15$、$l = 8$ のとき成立。$\checkmark$

次に $\text{GCD}(12,15)$ と $\text{GCD}(10,15)$ を計算します。

$$\begin{aligned} \text{GCD}(12,15) = 3 \\ \text{GCD}(10,15) = 5 \end{aligned}$$

これで $15$ の因数が求まりました！

アルゴリズム

$x^2 \equiv 1\bmod N$ を満たす整数 $x$ を見つけることが $N$ の因数分解にどう役立つかを確認したところで、ショアのアルゴリズムを詳しく見ていきましょう。本質的には $x$ を求めることに帰着します：

1. ランダムな整数を選ぶ $1 < a < N$ を満たすランダムな整数 $a$ を選びます。

• $\text{GCD}(a, N)$ を古典的に計算します。
  • $\text{GCD}(a, N) > 1$ であれば、すでに因数を見つけたことになります。終了。
  • そうでなければ続けます。

2. $a$ の $N$ に対する位数 $r$ を求める $a^r \equiv 1 \pmod N$ を満たす最小の正の整数 $r$ を求めます。

3. 位数が偶数かどうかを確認する

• $r$ が奇数であれば、ステップ 1 に戻って新しい $a$ を選びます。
• $r$ が偶数であれば、ステップ 4 に進みます。

4. $x = a^{r/2} \bmod N$ を計算する

• $x \not\equiv 1 \pmod N$ かつ $x \not\equiv -1 \pmod N$ であることを確認します。
  • $x \equiv \pm 1 \pmod N$ であれば、ステップ 1 に戻って新しい $a$ を選びます。
• そうでなければ、GCD を計算して因数を取り出します：

$$p = \text{GCD}(x-1, N), \quad q = \text{GCD}(x+1, N)$$

これらが $N$ の自明でない因数となります。

5. 必要に応じて再帰的に因数分解する

• $p$ や $q$ が素数でない場合は、アルゴリズムを再帰的に適用して完全に因数分解します。
• すべての因数が素数になれば、因数分解は完了です。

この手順を見ると、量子コンピューターがなぜ必要なのかすぐにはわからないかもしれません。それが必要な理由は、ステップ 2 の「$a$ の $N$ に対する位数を求めること」が、古典的には非常に困難な問題だからです。計算量は $N$ の数に対して指数関数的に増加します。しかし量子コンピューターを使えば、量子位相推定を利用してこれを解くことができます。ステップ 4 の 2 つの整数の GCD を求めることは、古典的には比較的容易なことです。したがって、実際に量子コンピューターの力が必要なのは位数探索のステップだけです。因数分解問題は位数探索問題に「帰着」すると言います。

難しい部分：位数探索

次に、量子コンピューターを位数探索にどのように活用できるかを説明します。まず「位数」の意味を明確にしておきましょう。数学的な意味はすでに述べました：$a^r = 1 \pmod N$ を満たす最初の非ゼロの整数 $r$ です。ここではこの概念についてもう少し直感的に理解してみましょう。

$N$ が十分に小さければ、$a$ の各冪を計算し、その数の $N$ に対する剰余を取り、$a^r = 1 \text{mod}(N)$ を満たす冪 $r$ が見つかったところで止めることで、単純に位数を求めることができます。これが先ほどの例 $N=15$ でやったことです。いくつかのサンプル値 $a$ と $N$ に対するモジュラー冪のグラフを見てみましょう：

何か気づきましたか？これらは周期関数です！そして位数 $r$ は周期と同じです！つまり、位数探索は周期探索と等価です。

量子コンピューターは関数の周期を求めることに非常に適しています。そのために、量子位相推定（QPE）と呼ばれるアルゴリズムのサブルーチンを使います。QPE と量子フーリエ変換との関係については前のモジュールで説明しました。詳しい復習は、QFT モジュール、またはジョン・ワトラスの量子アルゴリズムコースの量子位相推定のレッスンを参照してください。ここではその手順の概要を説明します：

量子位相推定（QPE）では、ユニタリ $U$ とそのユニタリの固有状態 $|\psi\rangle$ から始めます。次に QPE を使って対応する固有値を近似します。演算子がユニタリなので、固有値は $e^{2\pi i \theta}$ の形をとります。したがって、固有値を求めることは周期関数における $\theta$ の値を求めることと等価です。Circuit は次のようになります：

制御 Qubit の数（上図の上部 $m$ 個の Qubit）が近似の精度を決定します。

ショアのアルゴリズムでは、QPE をユニタリ演算子 $M_a$ に対して使います：

$M_a|y\rangle \equiv |ay \mod N \rangle .$

ここで $|y\rangle$ は多 Qubit レジスタの計算基底状態を表し、Qubit の 2 進数の値が整数 $y$ に対応します。たとえば $N=15$、$y=2$ の場合、$|y\rangle$ は 4 Qubit の基底状態 $|0010\rangle$ で表されます（15 までの数を符号化するには 4 Qubit が必要なため）。（この概念が馴染みのない方は、量子状態の 2 進数符号化についての復習としてQiskit in Classrooms の入門モジュールを参照してください。）

次に、このユニタリの固有状態を特定する必要があります。状態 $|1\rangle$ から始めると、$U$ を連続して適用するたびにレジスタの状態が $a \pmod N$ 倍になり、$r$ 回適用後に再び状態 $|1\rangle$ に戻ることがわかります。たとえば $a=3$、$N=35$ の場合：

$$\begin{aligned} M_3|1\rangle &= |3\rangle & \\ M_3^2|1\rangle &= |9\rangle \\ M_3^3|1\rangle &= |27\rangle \\ & \vdots \\ M_3^{(r-1)}|1\rangle &= |12\rangle \\ M_3^r|1\rangle &= |1\rangle \end{aligned}$$

したがって、このサイクル内の状態の重ね合わせ（$|\psi_j\rangle$）で以下の形をとるものは：

$|\psi_j\rangle = \tfrac{1}{\sqrt{r}}\sum_{k=0}^{r-1}{e^{\frac{2 \pi i j k}{r}} |a^k \rangle}$

すべて $M_a$ の固有状態です。（これら以外にも固有状態はありますが、ここでは上記の形のものだけを考えます。）

理解度チェック

以下の問題を読み、答えを考えてから、三角形をクリックして解答を確認してください。

$a=2$、$N=15$ に対応するユニタリの固有状態を求めてください。

解答：

$$\begin{aligned} M_2|1\rangle &= |2\rangle & \\ M_2^2|1\rangle &= |4\rangle \\ M_2^3|1\rangle &= |8\rangle \\ M_2^4|1\rangle &= |1\rangle \\ \end{aligned}$$

よって、位数 $r=4$ です。関心のある固有状態は、上でサイクルされたすべての状態の等しい重ね合わせに様々な位相を掛けたものになります：

$$\begin{aligned} |\psi_0\rangle &= \frac{1}{2}(|1\rangle+|2\rangle+|4\rangle+|8\rangle) \\ |\psi_1\rangle &= \frac{1}{2}(e^{2 \pi i \frac{0}{4}}|1\rangle+e^{2 \pi i \frac{1}{4}}|2\rangle+e^{2 \pi i \frac{2}{4}}|4\rangle+e^{2 \pi i \frac{3}{4}}|8\rangle) \\ &= \frac{1}{2}(|1\rangle+i|2\rangle-|4\rangle-i|8\rangle) \\ |\psi_2\rangle &= \frac{1}{2}(e^{2 \pi i \frac{0}{4}}|1\rangle+e^{2 \pi i \frac{2}{4}}|2\rangle+e^{2 \pi i \frac{4}{4}}|4\rangle+e^{2 \pi i \frac{6}{4}}|8\rangle) \\ &= \frac{1}{2}(|1\rangle-|2\rangle+|4\rangle-|8\rangle) \\ |\psi_3\rangle &= \frac{1}{2}(e^{2 \pi i \frac{0}{4}}|1\rangle+e^{2 \pi i \frac{3}{4}}|2\rangle+e^{2 \pi i \frac{6}{4}}|4\rangle+e^{2 \pi i \frac{9}{4}}|8\rangle) \\ &= \frac{1}{2}(|1\rangle-i|2\rangle-|4\rangle+i|8\rangle) \\ \end{aligned}$$

これらの固有状態の 1 つに Qubit 状態を初期化できるとします（ネタバレ：実際にはできません。少なくとも簡単にはできません。その理由と代替手段についてはすぐに説明します）。そうすれば、QPE を使って対応する固有値 $\omega_j = e^{2 \pi i \theta_j}$（$\theta_j = \frac{j}{r}$）を推定できます。その後、簡単な式で位数 $r$ を決定できます：

$r = \frac{j}{\theta_j}.$

ただし、QPE は $\theta_j$ を推定するものであって、正確な値を与えるわけではありません。$r$ と $r+1$ を区別できる程度に精度の高い推定が必要です。制御 Qubit の数 $m$ が多いほど推定精度は高くなります。レッスン末尾の問題では、数 $N$ を因数分解するために必要な最小の $m$ を求めていただきます。

ここで 1 つの問題に向き合う必要があります。上記の $r$ の求め方の説明は、すべて固有状態 $|\psi_j\rangle = \tfrac{1}{\sqrt{r}}\sum_{k=0}^{r-1}{e^{\frac{2 \pi i j k}{r}} |a^k \rangle}$ の準備から始まります。しかし、$r$ を事前に知らなければこれを準備する方法がわかりません。論理が循環してしまいます。固有状態を初期化せずに固有値を推定する方法が必要です。

$M_a$ の固有状態から始める代わりに、初期状態を 2 進数の $|1\rangle$ に対応する $n$ Qubit 状態（$|000...01\rangle$）に準備することができます。この状態自体は明らかに $M_a$ の固有状態ではありませんが、すべての固有状態 $|\psi_k\rangle$ の重ね合わせになっています：

$|1\rangle = \frac{1}{\sqrt{r}} \sum\limits_{k=0}^{r-1}{|\psi_k\rangle}$

理解度チェック

以下の問いを読み、答えを考えてから三角をクリックして解答を確認してください。

$N=15$、$a=2$ の場合に前の確認問題で求めた固有状態の重ね合わせが $|1\rangle$ に等しいことを確かめてください。

解答:

4 つの固有状態は次のとおりです。

$$\begin{aligned} |\psi_0\rangle &= \frac{1}{2}(|1\rangle+|2\rangle+|4\rangle+|8\rangle) \\ |\psi_1\rangle &= \frac{1}{2}(|1\rangle+i|2\rangle-|4\rangle-i|8\rangle) \\ |\psi_2\rangle &= \frac{1}{2}(|1\rangle-|2\rangle+|4\rangle-|8\rangle) \\ |\psi_3\rangle &= \frac{1}{2}(|1\rangle-i|2\rangle-|4\rangle+i|8\rangle) \\ \end{aligned}$$

したがって、

$$\begin{aligned} \frac{1}{\sqrt{r}} \sum\limits_{k=0}^{r-1}{|\psi_k\rangle} &= \frac{1}{2}(|\psi_0\rangle + |\psi_1\rangle + |\psi_2\rangle + |\psi_3\rangle ) \\ &= \frac{1}{4}(|1\rangle+|2\rangle+|4\rangle+|8\rangle+|1\rangle+i|2\rangle-|4\rangle-i|8\rangle+|1\rangle-|2\rangle+|4\rangle-|8\rangle + |1\rangle-i|2\rangle-|4\rangle+i|8\rangle) \\ &= \frac{1}{4}(4|1\rangle) = |1\rangle \end{aligned}$$

これによってどのように位数 $r$ を求めることができるのでしょうか？出発状態が上記の形式のすべての固有状態の重ね合わせであることから、QPE アルゴリズムはこれらの固有状態に対応する各 $\theta_k$ を同時に推定します。したがって、最後に $m$ 個の制御 Qubit を測定すると、固有値 $k \in \{0,1,2,...,r-1\}$ のうちランダムに選ばれた 1 つの値に対応する $k/r$ の近似値が得られます。この Circuit を何度か繰り返して異なる $k$ の値のサンプルをいくつか取得すれば、$r$ をすばやく導き出せるようになります。

Qiskit で実装する

先述のとおり、現在のハードウェアでは RSA1024 のような大きな数を素因数分解できる段階には達していません。ここではアルゴリズムの動作を示すために小さな数を素因数分解してみましょう。このデモでは Shor のアルゴリズムのチュートリアル に掲載されているコードを簡略化したバージョンを使用します。詳細はそちらのチュートリアルを参照してください。

量子問題を解くための標準フレームワークである Qiskit パターンフレームワークを使ってアルゴリズムを実行します。このフレームワークは次の 4 つのステップで構成されています。

1. 問題を量子 Circuit にマッピングする
2. 量子ハードウェア上で実行できるよう Circuit を最適化する
3. 量子コンピューター上で Circuit を実行する
4. 測定結果を後処理する

1. マップ

$N=15$ を素因数分解します。互いに素な整数として $a=2$ を選びます。

まず、剰余乗算ユニタリ $M_a$ を実装する Circuit を構築する必要があります。これは実装全体の中で最も難しい部分であり、実装方法によっては計算コストが非常に高くなる場合があります。ここでは少し「ずる」をします。出発状態が $|1\rangle$ であることはわかっており、以前の確認問題から

$$\begin{aligned} M_2|1\rangle &= |2\rangle & \\ M_2|2\rangle &= |4\rangle \\ M_2|4\rangle &= |8\rangle \\ M_2|8\rangle &= |1\rangle \\ \end{aligned}$$

であることがわかっています。そのため、これら 4 つの状態に対して正しい操作を行い、その他の状態はそのままにするユニタリを構築します。これが「ずる」である理由は、$2 \bmod 15$ の位数に関する知識を利用してユニタリを簡略化しているからです。素因数が未知の数を実際に素因数分解しようとする場合、このような手法は使えません。

理解度チェック

以下の問いを読み、答えを考えてから三角をクリックして解答を確認してください。

$M_2$ 演算子が上記の状態をどのように変換するかの知識を活かして、2 つの Qubit の状態を入れ替える SWAP Gate の列からこの演算子を構成してください。（ヒント: 各状態 $|i\rangle$ を 2 進数で表すと助けになります。）

解答:

$M_2$ の各状態への作用を 2 進数で書き直してみましょう。

$$\begin{aligned} M_2|0001\rangle &= |0010\rangle \\ M_2|0010\rangle &= |0100\rangle \\ M_2|0100\rangle &= |1000\rangle \\ M_2|1000\rangle &= |0001\rangle \\ \end{aligned}$$

これらの変換はそれぞれ単純な SWAP で実現できます。$M_2|0001\rangle$ は Qubit $0$ と $1$ の状態を入れ替えることで得られます。$M_2|0010\rangle$ は Qubit $1$ と $2$ の状態を入れ替えることで得られます。以降も同様です。したがって、$M_2$ 行列を次の SWAP Gate の列に分解できます。

$$M_2 = SWAP(0,1)SWAP(1,2)SWAP(2,3)$$

演算子は右から左へ作用することに注意しながら、各状態に対して意図した効果が得られることを確かめましょう。

$$\begin{aligned} M_2|0001\rangle &= SWAP(0,1)SWAP(1,2)SWAP(2,3)|0001\rangle \\ &= SWAP(0,1)SWAP(1,2)|0001\rangle \\ &= SWAP(0,1)|0001\rangle \\ &=|0010\rangle \checkmark \\ M_2|0010\rangle &= SWAP(0,1)SWAP(1,2)SWAP(2,3)|0010\rangle \\ &= SWAP(0,1)SWAP(1,2)|0010\rangle \\ &= SWAP(0,1)|0100\rangle \\ &=|0100\rangle \checkmark \\ M_2|0100\rangle &= SWAP(0,1)SWAP(1,2)SWAP(2,3)|0100\rangle \\ &= SWAP(0,1)SWAP(1,2)|1000\rangle \\ &= SWAP(0,1)|1000\rangle \\ &=|1000\rangle \checkmark \\ M_2|1000\rangle &= SWAP(0,1)SWAP(1,2)SWAP(2,3)|1000\rangle \\ &= SWAP(0,1)SWAP(1,2)|0100\rangle \\ &= SWAP(0,1)|0010\rangle \\ &=|0001\rangle \checkmark \\ \end{aligned}$$

これで、この演算子と等価な Circuit を Qiskit でコーディングできます。

まず、必要なパッケージをインポートします。

# Import necessary packages

import numpy as np
from fractions import Fraction
from math import floor, gcd, log

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.circuit.library import QFTGate
from qiskit.transpiler import generate_preset_pass_manager
from qiskit.visualization import plot_histogram

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

次に、$M_2$ 演算子を作成します。

def M2mod15():
    """
    M2 (mod 15)
    """
    b = 2
    U = QuantumCircuit(4)

    U.swap(2, 3)
    U.swap(1, 2)
    U.swap(0, 1)

    U = U.to_gate()
    U.name = f"M_{b}"

    return U

# Get the M2 operator
M2 = M2mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(4)
circ.compose(M2, inplace=True)
circ.decompose(reps=2).draw(output="mpl", fold=-1)

QPE アルゴリズムは controlled-$U$ Gate を使用します。$M_2$ Circuit が完成したので、次はこれを controlled-$M_2$ Circuit にする必要があります。

def controlled_M2mod15():
    """
    Controlled M2 (mod 15)
    """
    b = 2
    U = QuantumCircuit(4)

    U.swap(2, 3)
    U.swap(1, 2)
    U.swap(0, 1)

    U = U.to_gate()
    U.name = f"M_{b}"
    c_U = U.control()

    return c_U

# Get the controlled-M2 operator
controlled_M2 = controlled_M2mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(5)
circ.compose(controlled_M2, inplace=True)
circ.decompose(reps=1).draw(output="mpl", fold=-1)

これで controlled-$U$ Gate が完成しました。しかし量子位相推定アルゴリズムを実行するには、controlled-$U^2$、controlled-$U^4$、そして controlled-$U^{2^{m-1}}$ まで必要です。ここで $m$ は位相の推定に使用する Qubit 数を表します。Qubit 数が多いほど位相推定の精度が高くなります。今回は位相推定の手続きに $m=8$ 個の制御 Qubit を使用します。したがって、次の式が必要です。

$$M_{a^{2^k}}|y\rangle \equiv |a^{2^k} y \bmod N \rangle$$

ここでインデックス $k$（$0 \le k \le m-1 = 7$）は制御 Qubit に対応します。それでは各 $k$ の値について $a^{2^k}\bmod N$ を計算しましょう。

def a2kmodN(a, k, N):
    """Compute a^{2^k} (mod N) by repeated squaring"""
    for _ in range(k):
        a = int(np.mod(a**2, N))
    return a

k_list = range(8)
b_list = [a2kmodN(2, k, 15) for k in k_list]

print(b_list)

[2, 4, 1, 1, 1, 1, 1, 1]

$k \ge 2$ では $a^{2^k} \bmod N = 1$ となるため、対応する演算子（$M_8$ 以降）はすべて恒等演算子と等価です。したがって、新たに構築する必要があるのは $M_4$ のみです。

注意: この簡略化が成り立つのは、ここでは $2 \bmod 15$ の位数が $4$ であるためです。$k=2$（つまり $2^k = 4$）になると、それ以降のすべての演算子の累乗は恒等演算子になります。一般に、より大きな数 $N$ や異なる $a$ の選択では、高次の累乗の構築を省略することはできません。これがこの例がおもちゃの例と見なされる理由の 1 つです。数が小さいため、より大きなケースでは使えない近道が取れるのです。

def M4mod15():
    """
    M4 (mod 15)
    """
    b = 4
    U = QuantumCircuit(4)

    U.swap(1, 3)
    U.swap(0, 2)

    U = U.to_gate()
    U.name = f"M_{b}"

    return U

# Get the M4 operator
M4 = M4mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(4)
circ.compose(M4, inplace=True)
circ.decompose(reps=2).draw(output="mpl", fold=-1)

先ほどと同様に、これを controlled-$M_4$ 演算子にします。

def controlled_M4mod15():
    """
    Controlled M4 (mod 15)
    """
    b = 4
    U = QuantumCircuit(4)

    U.swap(1, 3)
    U.swap(0, 2)

    U = U.to_gate()
    U.name = f"M_{b}"
    c_U = U.control()

    return c_U

# Get the controlled-M4 operator
controlled_M4 = controlled_M4mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(5)
circ.compose(controlled_M4, inplace=True)
circ.decompose(reps=1).draw(output="mpl", fold=-1)

これですべてを組み合わせて、位相推定を使った量子 Circuit で $2\bmod 15$ の位数を求めることができます。

# Order finding problem for N = 15 with a = 2
N = 15
a = 2

# Number of qubits
num_target = floor(log(N - 1, 2)) + 1  # for modular exponentiation operators
num_control = 2 * num_target  # for enough precision of estimation

# List of M_b operators in order
k_list = range(num_control)
b_list = [a2kmodN(2, k, 15) for k in k_list]

# Initialize the circuit
control = QuantumRegister(num_control, name="C")
target = QuantumRegister(num_target, name="T")
output = ClassicalRegister(num_control, name="out")
circuit = QuantumCircuit(control, target, output)

# Initialize the target register to the state |1>
circuit.x(num_control)

# Add the Hadamard gates and controlled versions of the
# multiplication gates
for k, qubit in enumerate(control):
    circuit.h(k)
    b = b_list[k]
    if b == 2:
        circuit.compose(
            M2mod15().control(), qubits=[qubit] + list(target), inplace=True
        )
    elif b == 4:
        circuit.compose(
            M4mod15().control(), qubits=[qubit] + list(target), inplace=True
        )
    else:
        continue  # M1 is the identity operator

# Apply the inverse QFT to the control register
circuit.compose(QFTGate(num_control).inverse(), qubits=control, inplace=True)

# Measure the control register
circuit.measure(control, output)

circuit.draw("mpl", fold=-1)

2. 最適化

回路のマッピングが完了したら、次のステップは特定の量子コンピューターで実行できるよう回路を最適化することです。まず、Backendを読み込む必要があります。

service = QiskitRuntimeService()

backend = service.backend("ibm_marrakesh")

アカウントの利用時間が不足している場合、または何らかの理由でシミュレーターを使いたい場合は、以下のセルを実行することで、上記で選択した量子デバイスを模倣するシミュレーターをセットアップできます。

pm = generate_preset_pass_manager(optimization_level=2, backend=backend)

transpiled_circuit = pm.run(circuit)

print(f"2q-depth: {transpiled_circuit.depth(lambda x: x.operation.num_qubits==2)}")
print(f"2q-size: {transpiled_circuit.size(lambda x: x.operation.num_qubits==2)}")
print(f"Operator counts: {transpiled_circuit.count_ops()}")
transpiled_circuit.draw(output="mpl", fold=-1, style="clifford", idle_wires=False)

2q-depth: 188
2q-size: 281
Operator counts: OrderedDict({'sx': 548, 'rz': 380, 'cz': 281, 'measure': 8, 'x': 6})

3. 実行

# Sampler primitive to obtain the probability distribution
sampler = Sampler(backend)

# Turn on dynamical decoupling with sequence XpXm
sampler.options.dynamical_decoupling.enable = True
sampler.options.dynamical_decoupling.sequence_type = "XpXm"
# Enable gate twirling
sampler.options.twirling.enable_gates = True

pub = transpiled_circuit
job = sampler.run([pub], shots=1024)

result = job.result()[0]
counts = result.data["out"].get_counts()

plot_histogram(counts, figsize=(35, 5))

00000000、01000000、10000000、11000000 に明確な4つのピークが見られます。量子コンピューターのノイズによって、他のビット列にも若干のカウントが存在します。これらは無視し、閾値を設けることで上位の4つのみを保持します。この閾値を超えるカウントのみを、ノイズを上回る真のシグナルとみなします。

# Dictionary of bitstrings and their counts to keep
counts_keep = {}
# Threshold to filter
threshold = np.max(list(counts.values())) / 2

for key, value in counts.items():
    if value > threshold:
        counts_keep[key] = value

print(counts_keep)

4. 後処理

ショアのアルゴリズムでは、アルゴリズムのかなりの部分が古典的に実行されます。そのため、量子コンピューターから測定結果を取得した後の「後処理」ステップに残りの処理をまとめます。上記の各測定値は整数に変換でき、$2^m$ で割ることで $\frac{k}{r}$ の近似値が得られます。ここで $k$ は毎回ランダムに決まります。

a = 2
N = 15

FACTOR_FOUND = False
num_attempt = 0

while not FACTOR_FOUND:
    print(f"\nATTEMPT {num_attempt}:")
    # Here, we get the bitstring by iterating over outcomes
    # of a previous hardware run with multiple shots.
    # Instead, we can also perform a single-shot measurement
    # here in the loop.
    bitstring = list(counts_keep.keys())[num_attempt]
    num_attempt += 1
    # Find the phase from measurement
    decimal = int(bitstring, 2)
    phase = decimal / (2**num_control)  # phase = k / r
    print(f"Phase: theta = {phase}")

    # Guess the order from phase
    frac = Fraction(phase).limit_denominator(N)
    r = frac.denominator  # order = r
    print(f"Order of {a} modulo {N} estimated as: r = {r}")

    if phase != 0:
        # Guesses for factors are gcd(a^{r / 2} ± 1, 15)
        if r % 2 == 0:
            x = pow(a, r // 2, N) - 1
            d = gcd(x, N)
            if d > 1:
                FACTOR_FOUND = True
                print(f"*** Non-trivial factor found: {x} ***")

ATTEMPT 0:
Phase: theta = 0.0
Order of 2 modulo 15 estimated as: r = 1

ATTEMPT 1:
Phase: theta = 0.75
Order of 2 modulo 15 estimated as: r = 4
*** Non-trivial factor found: 3 ***

まとめ

このモジュールを通じて、ピーター・ショアがこれほど巧妙なアルゴリズムを考案したことへの新たな敬意を感じていただけたかもしれません。しかし同時に、その表面的なシンプルさの裏に潜む深さへの理解も深まったのではないでしょうか。このアルゴリズムは印象的（あるいは圧倒的）な複雑さを持つように見えますが、論理の各ステップに分解してゆっくりと辿ることで、あなた自身もショアのアルゴリズムを実行できるようになります。

RSA1024 のような大きな数を因数分解するためにこのアルゴリズムを使うには程遠い状況ですが、量子コンピューターは日々進歩しており、フォールトトレランスと呼ばれる閾値に達したとき、こうしたアルゴリズムも実用化の時代が訪れるでしょう。量子コンピューティングを学ぶ今は、非常にエキサイティングな時代です！

問題

重要な概念

• 現代の暗号システムは、大きな整数の因数分解が古典的に困難であることに依存しています。
• 剰余算術（$\mathbb{Z}_N$ および $\mathbb{Z}_N^*$ の構造を含む）は、ショアのアルゴリズムの数学的基盤を提供します。
• 整数 $N$ の因数分解の問題は、$N$ を法とする数の位数を求める問題に帰着できます。
• 量子位数探索は、量子位相推定の手法を用いて関数 $a^x \mod N$ の周期を決定します。
• ショアのアルゴリズムは、底を選択し、量子位数探索を実行し、その結果から古典的に因数を計算するという古典・量子ハイブリッドのワークフローで構成されています。

正誤問題

1. 正/誤　ショアのアルゴリズムの効率性は、RSA暗号のセキュリティを脅かします。
2. 正/誤　ショアのアルゴリズムは、現代のどの量子コンピューターでも効率的に実行できます。
3. 正/誤　ショアのアルゴリズムは、量子位相推定（QPE）を主要なサブルーチンとして使用します。
4. 正/誤　ショアのアルゴリズムの古典的な部分には、最大公約数（GCD）の計算が含まれます。
5. 正/誤　ショアのアルゴリズムは偶数の因数分解にのみ機能します。
6. 正/誤　ショアのアルゴリズムの実行が成功すれば、常に正しい因数が保証されます。

短答問題

1. ショアのアルゴリズムがRSA暗号に対する潜在的な将来の脅威と考えられているのはなぜですか？
2. 剰余指数関数の周期（位数）を求めることが、ショアのアルゴリズムにおける数の因数分解にどのように役立つのですか？

発展問題

1. QPEで位数 $r$ の正しい値を求めるために必要な精度を得るには、因数分解しようとする数 $N$ に対して制御 Qubit $m$ がいくつ必要ですか？

2. ここで概説した15の因数分解の手順に従い、今度は21を因数分解してみてください。