ショアのアルゴリズム

推定使用時間：Eagle r3 プロセッサーで3秒（注：これはあくまでも目安です。ランタイムは異なるかもしれない)。

学習成果

このチュートリアルを終えた後、ユーザーは以下の点を理解できるようになります：

整数因数分解のためのショアのアルゴリズムの数学的背景
このアルゴリズムのサンプルインスタンスをハードウェア上で実行する方法

前提条件

このチュートリアルを進める前に、以下の内容についてあらかじめ理解しておいていただくことをお勧めします：

量子アルゴリズムの基礎。
位相推定と因数分解。このチュートリアルでは、その内容の一部を取り上げます。

背景

1994年にピーター・ショアによって開発されたショアのアルゴリズムは、整数を多項式時間で因数分解する画期的な量子アルゴリズムである。その重要性は、既知のどの古典的アルゴリズムよりも指数関数的に高速に大きな整数を因数分解できる点にあり、大きな数の因数分解の困難さに依存するRSAのような広く利用されている暗号システムの安全性を脅かしている。十分に高性能な量子コンピュータを用いてこの問題を効率的に解くことができれば、ショアのアルゴリズムは暗号学、サイバーセキュリティ、計算数学といった分野に革命をもたらし、量子計算が持つ変革的な力を浮き彫りにすることになるだろう。

このチュートリアルでは、量子コンピューター上で15の因数分解を行うことにより、ショールのアルゴリズムを実証することに焦点を当てる。

まず、秩序発見問題を定義し、量子位相推定プロトコルから対応する回路を構成する。次に、トランスパイル可能な最短深度回路を使って、実際のハードウェア上で順序探索回路を実行する。最後のセクションでは、順序探索問題を整数因数分解に結びつけることで、ショールのアルゴリズムを完成させる。

このチュートリアルの最後には、実際のハードウェア上でのShorのアルゴリズムの他のデモについて、一般的な実装と、15や21のような特定の整数の因数分解に合わせた実装の両方に焦点を当てて議論する。

注：このチュートリアルは、ショールのアルゴリズムに関する回路の実装とデモンストレーションに重点を置いています。この教材に関する詳細な教育リソースについては、ジョン・ワトラス博士による量子アルゴリズムの基礎コース、および参考文献セクションの論文を参照してください。

要件

このチュートリアルを始める前に、以下のものがインストールされていることを確認してください：

Qiskit SDK v2.0 またはそれ以降、可視化サポート付き
Qiskit Runtime v0.40 またはそれ以降 (pip install qiskit-ibm-runtime)

セットアップ

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

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.circuit.library import QFT, UnitaryGate
from qiskit.transpiler import CouplingMap, generate_preset_pass_manager
from qiskit.visualization import plot_histogram

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

ステップ1：古典的な入力を量子問題にマッピングする

整数因数分解のためのショールのアルゴリズムは、 次数発見問題として知られる中間問題を利用する。このセクションでは、 量子位相推定を用いた次数探索問題の解き方を示す。

位相推定問題

位相推定問題では、 $n$ 量子ビットの量子状態 $\ket{\psi}$ 、 $n$ 量子ビットに作用するユニタリー量子回路が与えられる。私たちは、 $\ket{\psi}$ が回路の作用を記述するユニタリー行列 $U$ の固有ベクトルであることを約束している。私たちの目標は、 $\ket{\psi}$ が対応する固有値 $\lambda = e^{2 \pi i \theta}$ を計算または近似することである。言い換えれば、この回路は $\theta \in [0, 1)$ を満たす近似値を出力しなければならない。 $U \ket{\psi}= e^{2 \pi i \theta} \ket{\psi}.$ 位相推定回路の目標は、 $\theta$ を $m$ ビットで近似することである。数学的に言えば、 $\theta \approx y / 2^m$ 、ここで $y \in {0, 1, 2, \dots, 2^{m-1}}$ となるような $y$ を見つけたい。次の図は、 $m$ の量子ビットを測定することで、 $m$ の $y$ を推定する量子回路を示しています。

上記の回路では、トップ $m$ の量子ビットが $\ket{0^m}$ の状態で開始され、ボトム $n$ の量子ビットが $\ket{\psi}$ の状態で開始されます。これは $U$ の固有ベクトルであることが約束されています。位相推定回路の最初の構成要素は、対応する制御量子ビットへの位相キックバックを実行する制御されたユニタリ演算です。これらの制御されたユニタリーは、最下位ビットから最上位ビットまで、制御量子ビットの位置に応じて指数化される。 $\ket{\psi}$ は $U$ の固有ベクトルであるため、ボトム $n$ の量子ビットの状態はこの操作の影響を受けないが、固有値の位相情報はトップ $m$ の量子ビットに伝搬する。

制御されたユニタリーによる位相キックバック操作の後、ユニタリー $U$ の各固有ベクトル $\ket{\psi}$ に対して、トップ $m$ の量子ビットのすべての可能な状態が互いに直交正規であることがわかった。したがって、これらの状態は完全に区別可能であり、それらが形成する基底を回転させて計算基底に戻し、測定を行うことができる。数学的解析によれば、この回転行列は $2^m$ -次元ヒルベルト空間における逆量子フーリエ変換（QFT）に対応する。この背景には、モジュラー指数化演算子の周期構造が量子状態に符号化されており、QFTはこの周期性を周波数領域で測定可能なピークに変換するという直観がある。

ショールのアルゴリズムにQFT回路が採用されている理由をより深く理解するためには、「量子アルゴリズムの基礎」コースを参照されたい。

これで、次数探索に位相推定回路を使う準備が整った。

注文検索の問題

次数発見問題を定義するために、まず数論の概念から始める。まず、任意の正の整数 $N$ に対して、集合 $\mathbb{Z}_N$ を次のように定義する。 $\mathbb{Z}_N = \{0, 1, 2, \dots, N-1\}.$ $\mathbb{Z}_N$ のすべての算術演算は $N$ をモジュロとして行われる。特に、 $N$ とコプライムであるすべての要素 $a \in \mathbb{Z}_n$ は特別であり、 $\mathbb{Z}^*_N$ を構成する。 $\mathbb{Z}^*_N = \{ a \in \mathbb{Z}_N : \mathrm{gcd}(a, N)=1 \}.$ ある要素 $a \in \mathbb{Z}^*_N$ に対して、 $a^r \equiv 1 \; (\mathrm{mod} \; N)$ が $a$ の次数として定義されるような最小の正の整数 $r$ は、 $N$ をモジュロしたものである。後で見るように、 $a \in \mathbb{Z}^*_N$ の次数を求めれば、 $N$ を因数分解することができる。

位相推定回路から次数検出回路を構成するには、2つの考慮が必要である。まず、次数 $r$ を求めるためのユニタリー $U$ を定義する必要がある。次に、位相推定回路の初期状態を準備するために、 $U$ の固有ベクトル $\ket{\psi}$ を定義する必要がある。

次数発見問題を位相推定につなげるために、古典状態が $\mathbb{Z}_N$ に対応する系で定義される演算を考える。ここで、固定要素 $a \in \mathbb{Z}^*_N$ を乗算する。特に、この乗算演算子 $M_a$ を、各 $x \in \mathbb{Z}_N$ に対して $M_a \ket{x} = \ket{ax \; (\mathrm{mod} \; N)}$ となるように定義する。式の右辺のケットの内部で、 $N$ のモジュロ積をとっていることが暗黙の了解となっていることに注意。数学的分析によれば、 $M_a$ はユニタリー演算子である。さらに、 $M_a$ には固有ベクトルと固有値のペアがあり、 $a$ の次数 $r$ を位相推定問題に結びつけることができる。具体的には、 $j \in \{0, \dots, r-1\}$ の任意の選択に対して、 $\ket{\psi_j} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \omega^{-jk}_{r} \ket{a^k}$ は $M_a$ の固有ベクトルで、対応する固有値は $\omega^{j}_{r}$ である。 $\omega^{j}_{r} = e^{2 \pi i \frac{j}{r}}.$

観察によると、便利な固有ベクトル／固有値のペアは、 $\omega^{1}_{r} = e^{2 \pi i \frac{1}{r}}$ を持つ状態 $\ket{\psi_1}$ であることがわかる。したがって、もし固有ベクトル $\ket{\psi_1}$ を見つけることができれば、量子回路で位相 $\theta=1/r$ を推定することができ、したがって次数 $r$ の推定値を得ることができる。しかし、それは簡単ではないので、別の方法を考える必要がある。

初期状態として計算状態 $\ket{1}$ を用意した場合、回路はどのような結果になるかを考えてみよう。これは $M_a$ の固有状態ではないが、上で説明した固有状態の一様な重ね合わせである。つまり、以下の関係が成り立つ。 $\ket{1} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \ket{\psi_k}$ 上式の意味するところは、初期状態を $\ket{1}$ に設定すれば、 $k \in \{ 0, \dots, r-1\}$ を一様にランダムに選び、 $\ket{\psi_k}$ を位相推定回路の固有ベクトルとして使った場合と全く同じ測定結果が得られるということである。言い換えれば、トップ $m$ 量子ビットの測定は、値 $k / r$ （ $k \in \{ 0, \dots, r-1\}$ は一様にランダムに選ばれる）に対する近似値 $y / 2^m$ をもたらす。これによって、私たちの目標であった、何度かの独立した実行の後、高い信頼性をもって $r$ 。

モジュラー累乗演算子

これまでのところ、量子回路の $U = M_a$ と $\ket{\psi} = \ket{1}$ を定義することで、位相推定問題を次数探索問題と結びつけてきた。したがって、最後に残る要素は、 $M_a$ のモジュラー指数 $M_a^k$ for $k = 1, 2, 4, \dots, 2^{m-1}$ を定義する効率的な方法を見つけることである。この計算を行うには、 $k$ の任意のべき乗について、 $M_a$ の回路に $k$ を反復するのではなく、 $b = a^k \; \mathrm{mod} \; N$ を計算し、 $M_b$ の回路を使うことで、 $M_a^k$ の回路を作ることができる。必要なのは2のべき乗だけなので、反復2乗を使えば古典的に効率よく計算できる。

ステップ2：量子ハードウェア実行に向けた問題の最適化

$N = 15$ と $a=2$ の具体的な例

ここで一旦立ち止まって、具体的な例について議論し、 $N=15$ の順序探索回路を構築することができる。 $N=15$ に対して自明でない $a \in \mathbb{Z}_N^*$ は $a \in \{2, 4, 7, 8, 11, 13, 14 \}$ である。この例では $a=2$ を選ぶ。 $M_2$ 演算子とモジュラー指数演算子 $M_2^k$ を構成する。

計算基礎状態に対する $M_2$ の作用は以下の通り。 $M_2 \ket{0} = \ket{0} \quad M_2 \ket{5} = \ket{10} \quad M_2 \ket{10} = \ket{5}$ $M_2 \ket{1} = \ket{2} \quad M_2 \ket{6} = \ket{12} \quad M_2 \ket{11} = \ket{7}$ $M_2 \ket{2} = \ket{4} \quad M_2 \ket{7} = \ket{14} \quad M_2 \ket{12} = \ket{9}$ $M_2 \ket{3} = \ket{6} \quad M_2 \ket{8} = \ket{1} \quad M_2 \ket{13} = \ket{11}$ $M_2 \ket{4} = \ket{8} \quad M_2 \ket{9} = \ket{3} \quad M_2 \ket{14} = \ket{13}$ 観察によって、基底状態がシャッフルされていることがわかる。スワップゲートを使えば、4量子ビットでこの演算を構成できる。以下では、 $M_2$ と controlled- $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)

Output:

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)

Output:

つ以上の量子ビットに作用するゲートは、さらに2量子ビットのゲートに分解される。

circ.decompose(reps=2).draw(output="mpl", fold=-1)

Output:

ここでモジュラー指数演算子を構成する必要がある。位相推定に十分な精度を得るため、推定測定には8量子ビットを使用する。従って、各 $k = 0, 1, \dots, 7$ に対して $b = a^{2^k} \; (\mathrm{mod} \; N)$ で $M_b$ を構成する必要がある。

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)

Output:

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

$b$ の値のリストからわかるように、先に構築した $M_2$ に加えて、 $M_4$ と $M_1$ も構築する必要がある。 $M_1$ は計算基底状態に対して三重に作用するので、単純に同一演算子であることに注意。

$M_4$ は、次のように計算基底状態に作用する。 $M_4 \ket{0} = \ket{0} \quad M_4 \ket{5} = \ket{5} \quad M_4 \ket{10} = \ket{10}$ $M_4 \ket{1} = \ket{4} \quad M_4 \ket{6} = \ket{9} \quad M_4 \ket{11} = \ket{14}$ $M_4 \ket{2} = \ket{8} \quad M_4 \ket{7} = \ket{13} \quad M_4 \ket{12} = \ket{3}$ $M_4 \ket{3} = \ket{12} \quad M_4 \ket{8} = \ket{2} \quad M_4 \ket{13} = \ket{7}$ $M_4 \ket{4} = \ket{1} \quad M_4 \ket{9} = \ket{6} \quad M_4 \ket{14} = \ket{11}$

したがって、この順列は以下のスワップ操作で構成できる。

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)

Output:

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)

Output:

つ以上の量子ビットに作用するゲートは、さらに2量子ビットのゲートに分解される。

circ.decompose(reps=2).draw(output="mpl", fold=-1)

Output:

与えられた $b \in \mathbb{Z}^*_N$ に対する $M_b$ 演算子は並べ替え演算であることを見た。 $N=15$ 、必要な量子ビットはわずか4個であるため、ここで扱う並べ替え問題のサイズは比較的小さく、検査によって SWAP ゲートでこれらの演算を直接合成することができた。一般的に、これはスケーラブルなアプローチではないかもしれない。その代わりに、順列行列を明示的に構築し、Qiskitの UnitaryGate クラスと転置メソッドを使用して、この順列行列を合成する必要があるかもしれません。しかし、これでは回路がかなり深くなってしまう。以下に例を示します。

def mod_mult_gate(b, N):
    """
    Modular multiplication gate from permutation matrix.
    """
    if gcd(b, N) > 1:
        print(f"Error: gcd({b},{N}) > 1")
    else:
        n = floor(log(N - 1, 2)) + 1
        U = np.full((2**n, 2**n), 0)
        for x in range(N):
            U[b * x % N][x] = 1
        for x in range(N, 2**n):
            U[x][x] = 1
        G = UnitaryGate(U)
        G.name = f"M_{b}"
        return G

# Let's build M2 using the permutation matrix definition
M2_other = mod_mult_gate(2, 15)

# Add it to a circuit
circ = QuantumCircuit(4)
circ.compose(M2_other, inplace=True)
circ = circ.decompose()

# Transpile the circuit and get the depth
coupling_map = CouplingMap.from_line(4)
pm = generate_preset_pass_manager(coupling_map=coupling_map)
transpiled_circ = pm.run(circ)

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

Output:

qubits: 4
2q-depth: 94
2q-size: 96
Operator counts: OrderedDict({'cx': 45, 'swap': 32, 'u': 24, 'u1': 7, 'u3': 4, 'unitary': 3, 'circuit-335': 1, 'circuit-338': 1, 'circuit-341': 1, 'circuit-344': 1, 'circuit-347': 1, 'circuit-350': 1, 'circuit-353': 1, 'circuit-356': 1, 'circuit-359': 1, 'circuit-362': 1, 'circuit-365': 1, 'circuit-368': 1, 'circuit-371': 1, 'circuit-374': 1, 'circuit-377': 1, 'circuit-380': 1})

これらのカウントを、 $M_2$ ゲートの手動実装のコンパイル済み回路深度と比較してみよう。

# Get the M2 operator from our manual construction
M2 = M2mod15()

# Add it to a circuit
circ = QuantumCircuit(4)
circ.compose(M2, inplace=True)
circ = circ.decompose(reps=3)

# Transpile the circuit and get the depth
coupling_map = CouplingMap.from_line(4)
pm = generate_preset_pass_manager(coupling_map=coupling_map)
transpiled_circ = pm.run(circ)

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

Output:

qubits: 4
2q-depth: 9
2q-size: 9
Operator counts: OrderedDict({'cx': 9})

見てわかるように、順列行列のアプローチでは、 $M_2$ ゲート1つでも、手作業で実装したものに比べてかなり深い回路になった。従って、 $M_b$ 、これまでのオペレーションを継続する。

これで、先に定義した制御モジュラー指数演算子を使って、完全な次数発見回路を構成する準備が整った。以下のコードでは、Qiskit Circuitライブラリから QFT回路もインポートしている。この回路では、各クビットにハダマードゲート、一連の controlled-U1 （位相によってはZ）ゲート、スワップゲートのレイヤーを使用している。

# 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(QFT(num_control, inverse=True), qubits=control, inplace=True)

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

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

Output:

なお、 $M_1$ は恒等演算子であるため、残りの制御量子ビットからは制御されたモジュラー指数演算を省略した。

このチュートリアルの後半では、この回路を ibm_marrakesh バックエンドで実行することになる。そのために、この特定のバックエンドに従って回路をトランスパイルし、回路の深さとゲート数を報告する。

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
)

Output:

2q-depth: 187
2q-size: 260
Operator counts: OrderedDict({'sx': 521, 'rz': 354, 'cz': 260, 'measure': 8, 'x': 4})

ステップ3: `Qiskit primitives`を使用して実行する

まず、この回路を理想的なシミュレーターで実行した場合、理論的にどのような結果が得られるかを説明する。以下に、1024ショットを使用した上記回路のシミュレーション結果を示す。見てわかるように、制御量子ビット上の4つのビット列にほぼ一様な分布が得られる。

# Obtained from the simulator
counts = {"00000000": 264, "01000000": 268, "10000000": 249, "11000000": 243}

plot_histogram(counts)

Output:

制御量子ビットを測定することで、 $M_a$ 演算子の8ビット位相推定が得られる。この2進数表現を10進数に変換して、測定された位相を求めることができる。上のヒストグラムからわかるように、4つの異なるビット列が測定され、それぞれが以下のような位相値に対応している。

# Rows to be displayed in table
rows = []
# Corresponding phase of each bitstring
measured_phases = []

for output in counts:
    decimal = int(output, 2)  # Convert bitstring to decimal
    phase = decimal / (2**num_control)  # Find corresponding eigenvalue
    measured_phases.append(phase)
    # Add these values to the rows in our table:
    rows.append(
        [
            f"{output}(bin) = {decimal:>3}(dec)",
            f"{decimal}/{2 ** num_control} = {phase:.2f}",
        ]
    )

# Print the rows in a table
headers = ["Register Output", "Phase"]
df = pd.DataFrame(rows, columns=headers)
print(df)

Output:

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

$k$ が $\{0, 1, \dots, r-1 \}$ から一様にランダムにサンプリングされた $\theta = k / r$ に対応することを思い出してください。したがって、 $k$ と次数 $r$ を見つけようとするために、連続分数アルゴリズムを使うことができます。 Python にはこの機能が組み込まれています。例えば、 fractions モジュールを使って、浮動小数点数を Fraction オブジェクトに変換することができる：

Fraction(0.666)

Output:

Fraction(5998794703657501, 9007199254740992)

これは、正確に結果を返す分数（この場合は、 0.6660000... ）を与えるため、上記のような厄介な結果を与える可能性があります。 .limit_denominator() 、分母がある値以下の浮動小数点に最も近い分数を求めることができる：

# Get fraction that most closely resembles 0.666
# with denominator < 15
Fraction(0.666).limit_denominator(15)

Output:

Fraction(2, 3)

この方がずっといい。次数(r)はNより小さくなければならないので、分母の最大値を 15 とする：

# Rows to be displayed in a table
rows = []

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

# Print the rows in a table
headers = ["Phase", "Fraction", "Guess for r"]
df = pd.DataFrame(rows, columns=headers)
print(df)

Output:

   Phase Fraction  Guess for r
0   0.00      0/1            1
1   0.25      1/4            4
2   0.50      1/2            2
3   0.75      3/4            4

測定された固有値のうち2つが正しい結果を示していることがわかる： $r=4$ そして、Shorの次数探索アルゴリズムが失敗する可能性があることがわかる。これらの悪い結果は、 $k = 0$ 、あるいは $k$ と $r$ が共素数でなく、 $r$ の代わりに $r$ の因数が与えられているためである。これに対する最も簡単な解決策は、 $r$ について満足のいく結果が得られるまで実験を繰り返すことである。

これまで、位相推定回路を用いた $N=15$ （ $a=2$ ）の次数探索問題をシミュレータ上で実装した。ショールのアルゴリズムの最後のステップは、次数発見問題を整数因数分解問題に関連づけることである。このアルゴリズムの最後の部分は純粋に古典的なものであり、量子コンピューターから位相測定値を得た後、古典コンピューターで解くことができる。従って、アルゴリズムの最後の部分は、実際のハードウェア上でどのように順序探索回路を実行できるかを示した後まで延期する。

ハードウェアの実行

これで、以前 ibm_marrakesh についてトランスパイルした順序探索回路を実行することができる。ここでは、エラー抑制を目的としたダイナミック・デカップリング（DD）と、エラー緩和を目的としたゲート・ツワーリングに注目する。 DDは、正確にタイミングを合わせた制御パルスを量子デバイスに連続的に印加することで、不要な環境との相互作用やデコヒーレンスを効果的に平均化する。一方、ゲート・ツワーリングは、特定の量子ゲートをランダム化し、コヒーレント誤差をパウリ誤差に変換する。この2つの技術は、量子計算のコヒーレンスと忠実度を高めるためにしばしば組み合わされる。

# 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

# Assign tags before executing
sampler.options.environment.job_tags = ["TUT_SA"]

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))

Output:

見ての通り、同じビット列が最も高いカウント数で得られた。量子ハードウェアにはノイズがあるため、他のビット列へのリークもあるが、それは統計的にフィルタリングできる。

# 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)

Output:

{'00000000': 58, '01000000': 41, '11000000': 42, '10000000': 40}

ステップ4：後処理を行い、結果を希望の古典形式で返す

整数の因数分解

ここまでは、位相推定回路を使って次数探索問題を実装する方法について説明した。さて、次数探索問題を整数分解に結びつけると、ショールのアルゴリズムが完成する。アルゴリズムのこの部分は古典的であることに注意。

$N = 15$ と $a = 2$ を例に、このことを実証する。 $k / r$ ここで $a^r \; (\textrm{mod} \; N) = 1$ と $k$ は $0$ と $r - 1$ の間のランダムな整数である。この式から、 $(a^r - 1) \; (\textrm{mod} \; N) = 0,$ があり、これは $N$ が $a^r-1$ を割らなければならないことを意味する。 $r$ も偶数であれば、 $a^r -1 = (a^{r/2}-1)(a^{r/2}+1).$ と書くことができる。 $r$ が偶数でなければ、それ以上進めないので、 $a$ の値を変えて再挑戦しなければならない。そうでなければ、 $N$ と $a^{r/2}-1$ 、または $a^{r/2}+1$ の最大公約数が、 $N$ のべき乗である確率が高い。

統計的に失敗するアルゴリズムもあるので、 $N$ のファクターが少なくとも1つ見つかるまでこのアルゴリズムを繰り返す。

以下のセルは、 $N=15$ の因子が少なくとも1つ見つかるまでアルゴリズムを繰り返す。上記のハードウェアの実行結果を使用して、各反復における位相と対応する係数を推測する。

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} ***")

Output:


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

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

ディスカッション

アルゴリズムのスケーリング

RSA暗号化には通常、2048ビットから4096ビットのオーダーの鍵サイズが必要である。ショールのアルゴリズムで2048ビットの因数分解を行おうとすると、エラー訂正のオーバーヘッドを含めて数百万量子ビットの量子回路が必要となり、回路の深さは10億のオーダーになる。そのため、ショールのアルゴリズムが現代の暗号システムの解読に実用的であるためには、最適化された回路構成法か、ロバストな量子エラー訂正が必要となる。 Shorのアルゴリズムのリソース推定に関するより詳細な議論については、 [9] を参照されたい。

課題

チュートリアル終了おめでとう！今こそ理解度を試す絶好の機会だ。 21を因数分解する回路を作ってみてください $a$ 、お好きなものをお選びいただけます。アルゴリズムのビット精度を決めて量子ビット数を選択し、モジュラー指数演算子 $M_a$ を設計する必要があります。 [6] のFig.9と [7] のFig.2に示されている方法論について読んでみてください。

def M_a_mod21():
    """
    M_a (mod 21)
    """

    # Your code here
    pass

参照

Shor, Peter W. "量子コンピュータにおける素因数分解と離散対数の多項式時間アルゴリズム " SIAM review 41.2 (1999)：303-332.
IBM Quantum ジョン・ワトラウス博士による「量子アルゴリズムの基礎」講座。
Vandersypen, Lieven MK, et al. "Experimental realization of Shor's quantum factoring algorithm using nuclear magnetic resonance " Nature 414.6866 (2001)：883-887.
Martin-Lopez, Enrique, et al. "Experimental realization of Shor's quantum factoring algorithm using qubit recycling." Nature photonics 6.11 (2012)：773-776.
Monz, Thomas, et al. "Realization of scalable Shor algorithm " Science 351.6277 (2016)：1068-1070.
Amico, Mirko, Zain H. Saleem, and Muir Kumph. " IBM Q Experienceを用いたShorのファクタリングアルゴリズムの実験的研究 " Physical Review A 100.1 (2019)： 012305.
スコサナ、ウナティ、マーク・タメ。 "Demonstration of Shor's factoring algorithm for N=21 on IBM quantum processors." Scientific reports 11.1 (2021)： 16599.
キタエフ、A. Yu. "量子測定とアベル安定化問題 " arXiv preprint quant-ph/9511026 (1995).
ギドニー、クレイグ、マーティン・エケロ。 "2000万個のノイズの多い量子ビットを使って、8時間で2048ビットのRSA整数を因数分解する方法 " クォンタム5（2021年）： 433.

このページは役に立ちましたか？

バグや誤字の報告、またはコンテンツの要求はGitHubで行ってください。

ショアのアルゴリズム

学習成果

前提条件

背景

要件

セットアップ

ステップ1：古典的な入力を量子問題にマッピングする

位相推定問題

注文検索の問題

モジュラー累乗演算子

ステップ2：量子ハードウェア実行に向けた問題の最適化

N=15N = 15N=15 と a=2a=2a=2 の具体的な例

ステップ3: Qiskit primitivesを使用して実行する

ハードウェアの実行

ステップ4：後処理を行い、結果を希望の古典形式で返す

整数の因数分解

ディスカッション

関連研究

アルゴリズムのスケーリング

課題

参照

$N = 15$ と $a=2$ の具体的な例

ステップ3: `Qiskit primitives`を使用して実行する