この演習を円滑に行う為には下記のページの手順でまず開発環境を準備してください。
https://qiita.com/jg1mnv/items/810c6a31ea507db60292

開発環境が用意出来ましたら、この演習を新規のJupyterNotebookに貼り付けて演習を進めます。

Shorのアルゴリズムによる因数分解 (Qiskit 2.1対応版)

1. はじめに

このノートブックでは、Shorのアルゴリズムを使って整数を因数分解する量子アルゴリズムを学びます。このコードは、元々東京大学 ICEPP QCワークブックに掲載されていたものを、Python 3.11 および Qiskit 2.1.1 環境で動作するように更新し、解説を加えたものです。
※参考サイト
https://utokyo-icepp.github.io/qc-workbook/ja/shor.html

対象者:

●量子コンピューティングの基本を理解している方
●Shorのアルゴリズムの概念を学び、Qiskitでの実装例を見たい方

Shorのアルゴリズムとは:

Shorのアルゴリズムは、与えられた合成数 N を効率的に素因数分解する量子アルゴリズムです。古典コンピュータでは非常に時間がかかる大きな数の素因数分解を、量子コンピュータならば現実的な時間で解ける可能性を示したことで注目されています。その核心は、a^x mod N の周期を見つける問題（位数発見問題）を量子フーリエ変換 (QFT) を用いて効率的に解く点にあります。

このノートブックでは、比較的小さな数 N=15 の因数分解を例に、Shorのアルゴリズムのステップを追っていきます。

2. 必要なライブラリのインポート

まず、計算に必要なライブラリをインポートします。 qBraid環境では、これらのライブラリは通常プリインストールされています。

import matplotlib.pyplot as plt
import numpy as np
import math
from fractions import Fraction # Python標準ライブラリで連分数を扱う

# Qiskitのインポート
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile # <-- ここに追加
from qiskit.circuit.library import QFT
from qiskit_aer import AerSimulator  # Qiskit Aer のシミュレータ
from qiskit.visualization import plot_histogram

# Jupyter Notebook環境でプロットをインライン表示するためのマジックコマンド
%matplotlib inline

# バージョン情報確認のセルもそのまま実行してください
import qiskit
import qiskit_aer
import sys # Pythonのバージョン情報を取得するために追加

print(f"Python バージョン: {'.'.join(map(str, sys.version_info[:3]))}")
print(f"Qiskit バージョン: {qiskit.__version__}")
print(f"Qiskit Aer バージョン: {qiskit_aer.__version__}")

import qiskit
print(f"Qiskit バージョン: {qiskit.__version__}")
# Python 3.11, Qiskit 2.0.0 (qiskit-aer 0.14.0程度) を想定

3. Shorのアルゴリズムの構成要素

Shorのアルゴリズムは、いくつかの主要な量子回路コンポーネントと古典的な後処理から構成されます。ここでは、N=15 の因数分解を試みます。

3.1. 制御Uゲート (モジュラー指数関数)

Shorのアルゴリズムでは、a^x mod N を計算する操作が必要です。これは制御Uゲートとして実装されます。具体的には、制御量子ビットが |1⟩ の場合にのみ、ターゲット量子ビットに対して y → a*y mod N という操作（あるいは y → a^(2^j)*y mod N）を行うゲートを作成します。

元のコードでは、N=15 の場合に特化した a の値 (2, 7, 8, 11, 13) に対する U_a^k ( k は 2 のべき乗) を計算する制御ゲート c_amod15(a, k) を定義しています。この関数は a^k mod 15 を行うユニタリ操作の制御版を返します。ここで k は 2^j の形で与えられ、Shorのアルゴリズムの j 番目の制御量子ビットに対応します。

c_amod15(a, exponent) 関数の実装

この関数は、a を底とし、exponent を指数とするモジュラー指数演算 x → (a^exponent * x) mod 15 (実際には x=1 の状態から a^exponent mod 15 を生成する部分と、それに伴う他の入力への作用) を行うユニタリ演算 U を構築し、その制御版ゲート C-U を返します。元のコードでは、a の値に応じて特定のSWAPゲートシーケンスを exponent 回繰り返すことで U_{a^{exponent}} を実装しています。

Qiskit 2.0では、QuantumCircuit.to_gate() や Gate.control() メソッドが引き続き利用できます。

def c_amod15(a, exponent):
    """
    a^exponent mod 15 を計算する制御ゲートを作成します。
    この関数は、N=15 の場合に特化しており、a は {2, 7, 8, 11, 13} のいずれかである必要があります。
    操作は4量子ビットのターゲットレジスタに作用します。

    Args:
        a (int): 底 (2, 7, 8, 11, or 13)
        exponent (int): 指数 (通常は 2^j の形)

    Returns:
        qiskit.circuit.Gate: 制御されたモジュラー指数演算ゲート
    """
    if a not in [2, 7, 8, 11, 13]:
        raise ValueError("'a' must be 2, 7, 8, 11, or 13")

    # U_a (a*x mod 15) を exponent 回作用させる回路を構築
    U_circuit = QuantumCircuit(4) # ターゲット量子ビットは4つ

    for _ in range(exponent): # U_a を exponent 回繰り返す
        if a == 2:
            U_circuit.swap(0, 1)
            U_circuit.swap(1, 2)
            U_circuit.swap(2, 3)
        elif a == 7:
            U_circuit.swap(2, 3)
            U_circuit.swap(1, 2)
            U_circuit.swap(0, 1)
            for q_idx in range(4):
                U_circuit.x(q_idx)
        elif a == 8: # 8 = 2^{-1} mod 15 なので、U_8 = (U_2)^{-1}
            U_circuit.swap(2, 3)
            U_circuit.swap(1, 2)
            U_circuit.swap(0, 1)
        elif a == 11: # 11 = (2*7)^{-1} mod 15?  11*11 = 121 = 1 mod 15
                      # U_11 のSWAPシーケンス
            U_circuit.swap(0, 1)
            U_circuit.swap(2, 3)
            U_circuit.swap(1, 2)
        elif a == 13: # 13 = -2 mod 15. U_13 = X (U_2)^{-1} X (または X U_8 X)
            U_circuit.swap(1, 3)
            U_circuit.swap(0, 2)
            for q_idx in range(4):
                U_circuit.x(q_idx)

    # 回路をゲートに変換
    U_gate = U_circuit.to_gate()
    U_gate.name = f"U({a}^{exponent}mod15)"

    # 制御ゲートを作成 (制御量子ビット1つ)
    cU_gate = U_gate.control(1)
    return cU_gate

3.2. 逆量子フーリエ変換 (IQFT)

量子フーリエ変換 (QFT) は、Shorのアルゴリズムで周期を見つけるためのキーとなります。アルゴリズムの最後では、測定の前に逆量子フーリエ変換 (IQFT) を適用します。 Qiskitの qiskit.circuit.library.QFT を利用します。

def iqft_gate(n_qubits):
    """
    n量子ビットの逆量子フーリエ変換(IQFT)ゲートを作成します。
    QiskitのQFTライブラリを使用し、逆変換を指定します。
    do_swaps=True はQFTの標準的な定義（出力の量子ビット順序を反転させるSWAPを含む）に合わせます。
    """
    iqft = QFT(n_qubits, inverse=True, do_swaps=True, name="IQFT")
    return iqft

4. Shorのアルゴリズムの量子回路構築 (N=15)

Shorのアルゴリズムの量子部分は、以下のステップで構成されます。

1.量子ビットの準備:

●カウンティングレジスタ: 周期 r を推定するために使用。2n 量子ビットが目安（nはNを表すのに必要なビット数）。N=15 (< 2^4) なので n=4。カウンティングレジスタは 2n=8 量子ビットあれば十分な精度が期待できます。

●補助レジスタ (ターゲットレジスタ): a^x mod N の計算を行う。N=15 を格納するのに 4 量子ビット必要。

●古典レジスタ: カウンティングレジスタの測定結果を格納。

2.初期状態の準備:

●カウンティングレジスタの全量子ビットにアダマールゲートを適用し、重ね合わせ状態 |+⟩^⊗2n を作成。

●補助レジスタを |1⟩ の状態に初期化 (0...01)。N=15 の場合、|0001⟩。

3.制御モジュラー指数演算の適用:

●カウンティングレジスタの j 番目の量子ビットを制御ビットとし、補助レジスタに U_a^(2^j) 演算 ( a^(2^j) * y mod N ) を適用します。これを全てのカウンティング量子ビットに対して行います。

4.逆量子フーリエ変換の適用:

●カウンティングレジスタに逆量子フーリエ変換 (IQFT) を適用します。

5.測定:

●カウンティングレジスタを測定し、古典ビットに結果を格納します。

パラメータ設定

N=15 を因数分解します。a は N と互いに素な整数を選びます。ここでは a=7 を使用します。

N = 15
a = 7 # Nと互いに素な数を適当に選ぶ (例: 2, 7, 8, 11, 13 for N=15)

# カウンティング量子ビットの数
# N を表現するのに必要なビット数を n とすると (N < 2^n)
# 良い精度のためには 2n ビット程度あると望ましい
n_bits_N = math.ceil(math.log2(N))
n_counting_qubits = 2 * n_bits_N  # 例: N=15 (n=4) -> 8 カウンティング量子ビット
print(f"N = {N} (表現に必要なビット数 n = {n_bits_N})")
print(f"カウンティング量子ビット数: {n_counting_qubits}")

# 補助量子ビットの数 (Nを格納するのに必要)
n_aux_qubits = n_bits_N
print(f"補助量子ビット数: {n_aux_qubits}")

量子回路の構築

# 量子レジスタと古典レジスタの準備
counting_qreg = QuantumRegister(n_counting_qubits, name='counting')
aux_qreg = QuantumRegister(n_aux_qubits, name='aux') # a^x mod N の計算用
classical_creg = ClassicalRegister(n_counting_qubits, name='c')

# 量子回路の初期化
qc = QuantumCircuit(counting_qreg, aux_qreg, classical_creg)

# 1. 初期状態の準備
# カウンティングレジスタをHゲートで重ね合わせ状態に
for q_idx in range(n_counting_qubits):
    qc.h(counting_qreg[q_idx])

# 補助レジスタを |1> に初期化 (0...01)
# N=15 (4量子ビット) の場合、|0001> なので、aux_qreg[0] を |1> にする
qc.x(aux_qreg[0]) # aux_qreg[0] がLSB (最下位ビット) と仮定

qc.barrier() # 回路の見た目を分かりやすくするためのバリア

# 2. 制御モジュラー指数演算の適用
# counting_qreg[j] を制御ビットとして、(U_a)^(2^j) を aux_qreg に作用させる
for j in range(n_counting_qubits):
    exponent = 2**j
    cU_gate = c_amod15(a, exponent) # a^(2^j) mod 15 の制御ゲート

    # 制御ビットとターゲットビットを指定してゲートを適用
    # 制御ビット: counting_qreg[j]
    # ターゲットビット: aux_qreg全体
    control_qubit = counting_qreg[j]
    target_qubits = [aux_qreg[i] for i in range(n_aux_qubits)]
    
    qc.append(cU_gate, [control_qubit] + target_qubits)

qc.barrier()

# 3. 逆量子フーリエ変換 (IQFT) をカウンティングレジスタに適用
# IQFTゲートの適用 (do_swaps=Trueに注意)
# QiskitのQFTはデフォルトでdo_swaps=Trueなので、出力の量子ビット順序が反転される。
# IQFTを適用する量子ビットはカウンティングレジスタ
qc.append(iqft_gate(n_counting_qubits), counting_qreg)

qc.barrier()

# 4. カウンティングレジスタの測定
qc.measure(counting_qreg, classical_creg)

# 回路を描画
print("Shorのアルゴリズム用量子回路:")
qc.draw(output='mpl', fold=-1) # 'text', 'mpl', 'latex' など

補足: qc.draw('mpl') は Matplotlib を使用して描画します。環境によっては表示に時間がかかったり、フォント設定が必要な場合があります。'text' はアスキーアートで表示します。

5. シミュレーションの実行¶

構築した量子回路をシミュレータで実行し、測定結果を得ます。

# 5. シミュレーションの実行

# AerSimulator の準備
simulator = AerSimulator()

# ショット数 (測定回数)
shots = 2048

# --- 修正箇所ここから ---

# 回路 (qc) に含まれるカスタムゲートを明示的に分解します。
# decompose() メソッドは、回路内の "定義を持つ" カスタムゲートを
# その定義に基づいて1レベル展開します。
# これにより、cU_gate のようなカスタムゲートが、より基本的なゲートの組み合わせに置き換わります。
print("元の回路の深さ:", qc.depth()) # 参考情報
decomposed_qc = qc.decompose()
print("分解後の回路の深さ:", decomposed_qc.depth()) # 参考情報

# 次に、分解された回路をシミュレータ用にトランスパイルします。
# backend=simulator を指定することで、AerSimulator がサポートする
# 基本ゲートセットに最適化・変換されます。
print("シミュレータ用にトランスパイル中...")
transpiled_qc = transpile(decomposed_qc, backend=simulator)
print("トランスパイル後の回路の深さ:", transpiled_qc.depth()) # 参考情報

# オプション: トランスパイル後の回路の内容を確認（デバッグ用）
# print("\nトランスパイル後の回路:")
# print(transpiled_qc.draw(output='text', fold=-1))
# --- 修正箇所ここまで ---


print(f"\nシミュレータ ({simulator.name}) で {shots} ショット実行します...")
# トランスパイル済みの回路 (transpiled_qc) を実行します
job = simulator.run(transpiled_qc, shots=shots)
result = job.result()

# カウント結果を取得する際も、実行した回路オブジェクト (transpiled_qc) を指定するのが確実です。
# あるいは、result.get_counts() と引数なしで呼び出すと、最初の実験結果のカウントを返します。
counts = result.get_counts(transpiled_qc)
# または
# counts = result.get_counts() # これでも通常はOK

print("\n測定結果 (カウント):")
print(counts)

# 結果をヒストグラムで表示
plot_histogram(counts)
plt.show()

6. 結果の解析 (古典的な後処理)

量子計算で得られた測定結果 (カウンティングレジスタの値 m) から、周期 r を推定します。 m / 2^(n_counting_qubits) が s / r の良い近似になっていることを利用します (s はある整数)。この近似分数を連分数展開することで、周期 r の候補を見つけます。

6.1. 連分数展開による周期の推定

測定結果として得られたビット列 (例えば '10000000') は10進数に変換して m とします。 m / Q (ここで Q = 2^(n_counting_qubits)) を連分数展開し、その近似分数の分母が周期 r の候補となります。

以下の continued_fraction_expansion 関数と get_period_from_measurement 関数は、元のコードの ContinuedFractions クラスのロジックを参考に、測定結果から周期を推定する処理を実装します。

def continued_fraction_expansion(numerator, denominator):
    """
    分数 numerator/denominator の連分数展開の係数をリストとして返します。
    例: 7/15 = [0; 2, 7] -> [0, 2, 7]
    """
    result = []
    while denominator:
        integer_part = numerator // denominator
        result.append(integer_part)
        remainder = numerator % denominator
        numerator = denominator
        denominator = remainder
    return result

def get_convergents(cf_coeffs):
    """
    連分数展開の係数リストから、近似分数を生成します。
    戻り値は (分子, 分母) のリスト。
    """
    convergents_list = []
    # p_{-2}=0, p_{-1}=1, q_{-2}=1, q_{-1}=0
    p_prev2, p_prev1 = 0, 1
    q_prev2, q_prev1 = 1, 0
    for i, coeff_a in enumerate(cf_coeffs):
        p_curr = coeff_a * p_prev1 + p_prev2
        q_curr = coeff_a * q_prev1 + q_prev2
        
        # 0除算を避ける (q_currが0になることは稀だが念のため)
        if q_curr == 0:
            # 通常、連分数の係数が非負整数であればq_currは増加するはず
            # 例外的なケースや入力エラーを考慮
            print(f"Warning: Denominator q_curr is zero at step {i}. Coefficients: {cf_coeffs}")
            break

        convergents_list.append((p_curr, q_curr))
        
        p_prev2, p_prev1 = p_prev1, p_curr
        q_prev2, q_prev1 = q_prev1, q_curr
        
    return convergents_list

def get_period_candidates_from_measurement(measurement_decimal, q_val, max_period_N):
    """
    測定値 (10進数) から周期 r の候補を求めます。
    Args:
        measurement_decimal (int): 測定された値 (10進数)。
        q_val (int): カウンティングレジスタの取りうる状態の総数 (2^n_counting_qubits)。
        max_period_N (int): 探す周期の上限 (通常はN)。
    Returns:
        list: 周期rの候補リスト。
    """
    if measurement_decimal == 0: # m=0 の場合は情報なし (または周期が非常に大きい)
        return []

    # m/Q を分数として表現
    phase = Fraction(measurement_decimal, q_val)
    
    # 連分数展開 (PythonのFractionはlimit_denominatorで近似分数を得られる)
    # limit_denominator(max_denominator) は、分母が max_denominator 以下の最良近似分数を返す
    # ここでは周期rがN未満であると仮定して max_denominator = N とする
    approx_frac = phase.limit_denominator(max_period_N)
    
    r_candidate = approx_frac.denominator
    
    # 周期rは1より大きい必要がある
    if r_candidate > 1:
        return [r_candidate]
    else: # limit_denominatorが適切にrを見つけられない場合もある (例: r=1 になる)
          # より詳細な連分数展開からの候補探索も可能だが、ここではシンプルにする
          # 例えば、get_convergentsを使って全ての近似分母を調べるなど
        cf_coeffs = continued_fraction_expansion(measurement_decimal, q_val)
        convergents = get_convergents(cf_coeffs)
        potential_r = []
        for p, q in convergents:
            if q > 1 and q < max_period_N: # qが周期の候補
                # 検証: a^q mod N == 1 かどうか
                # ここでは候補としてqを返すに留める
                potential_r.append(q)
        return list(set(potential_r)) # 重複除去

6.2. 周期の検証と因数の発見

得られた周期 r の候補を使って、実際に a^r mod N == 1 であるかを確認します。もし r が偶数で、かつ a^(r/2) mod N != -1 mod N (つまり a^(r/2) mod N != N-1) ならば、 gcd(a^(r/2) - 1, N) と gcd(a^(r/2) + 1, N) が N の非自明な因数を与える可能性があります。

print(f"\n--- 周期の推定と因数分解 (N={N}, a={a}) ---")

# Q = 2^(カウンティング量子ビット数)
Q_val = 2**n_counting_qubits

# 最も頻繁に測定された結果から周期を推定
# (実際には複数の測定結果を試すのが良い)
sorted_counts = sorted(counts.items(), key=lambda item: item[1], reverse=True)

found_factors = False
tested_measurements_count = 0

for measured_str, count in sorted_counts:
    if found_factors:
        break
    
    # 測定結果が非常に少ない場合はスキップ (ノイズの可能性)
    if count < shots / 100 and tested_measurements_count > 5 : # 例: 全体の1%未満は無視（上位5件は試す）
        continue

    tested_measurements_count += 1
    measured_decimal = int(measured_str, 2) # 2進数文字列を10進数に
    
    print(f"\n測定結果: '{measured_str}' (10進数: {measured_decimal}, 回数: {count})")
    
    if measured_decimal == 0:
        print("  m=0 のため、周期推定スキップ。")
        continue

    # 周期rの候補を取得
    # N より大きい周期は意味がないので、max_period_N = N
    r_candidates = get_period_candidates_from_measurement(measured_decimal, Q_val, N)
    print(f"  m/Q = {measured_decimal}/{Q_val} ~= {Fraction(measured_decimal, Q_val)}")
    
    if not r_candidates:
        print("  この測定値から有効な周期候補が見つかりませんでした。")
        continue

    print(f"  周期rの候補: {r_candidates}")

    for r_cand in r_candidates:
        if r_cand == 0: continue # 周期0は無意味

        # 周期の検証: a^r_cand mod N == 1 か？
        if pow(a, r_cand, N) != 1:
            print(f"  候補 r={r_cand} は検証失敗 (a^r mod N != 1)。")
            continue
        
        print(f"  候補 r={r_cand} は検証成功 (a^r mod N == 1)。")

        # r が奇数の場合、このrからは因数を見つけられない
        if r_cand % 2 != 0:
            print(f"  r={r_cand} は奇数です。この周期からは直接因数を見つけられません。")
            continue

        # r が偶数の場合
        # x = a^(r/2) mod N
        # x == 1 または x == N-1 (-1 mod N) の場合、自明な因数しか得られない
        term = pow(a, r_cand // 2, N)
        if term == 1 or term == (N - 1):
            print(f"  r={r_cand} の場合、a^(r/2) mod N = {term} となり、自明な因数しか得られません。")
            continue
            
        # 非自明な因数が見つかる可能性
        factor1 = math.gcd(term - 1, N)
        factor2 = math.gcd(term + 1, N)

        print(f"  発見された因数候補: gcd({term-1}, {N}) = {factor1}, gcd({term+1}, {N}) = {factor2}")

        if factor1 != 1 and factor1 != N:
            print(f"  !!! 非自明な因数が見つかりました: {factor1} !!!")
            if N // factor1 == factor2 :
                 print(f"  N = {N} = {factor1} * {N//factor1}")
                 found_factors = True
                 break # 因数が見つかったらループを抜ける
        
        if factor2 != 1 and factor2 != N:
            print(f"  !!! 非自明な因数が見つかりました: {factor2} !!!")
            if N // factor2 == factor1:
                print(f"  N = {N} = {factor2} * {N//factor2}")
                found_factors = True
                break # 因数が見つかったらループを抜ける
    
if not found_factors:
    print("\n今回の試行では、Nの非自明な因数を見つけることができませんでした。")
    print("考えられる理由:")
    print("- `a` の選び方が悪かった (周期が奇数になるなど)。")
    print("- シミュレーションのショット数が不足し、正しい周期に対応する測定結果が得られなかった。")
    print("- 乱択アルゴリズムであるため、確率的に失敗することがある。")
    print("異なる `a` で試すか、ショット数を増やして再試行してください。")

7. まとめと実践

このノートブックでは、ShorのアルゴリズムのQiskit 2.0による実装例 (N=15) を見てきました。主なステップは以下の通りです:

1.適切な数の量子ビット（カウンティング用、補助用）を用意。
2.カウンティングレジスタをアダマールゲートで初期化。補助レジスタを |1⟩ に初期化。
3.制御モジュラー指数演算 (C-U_a^{2^j}) を適用。
4.カウンティングレジスタに逆量子フーリエ変換 (IQFT) を適用。
5.カウンティングレジスタを測定。
6.測定結果から連分数展開を用いて周期 r を推定。
7.周期 r から N の因数を計算。

注意点:

●Shorのアルゴリズムは確率的アルゴリズムであり、一度の実行で必ずしも正しい周期や因数が見つかるとは限りません。複数回の試行や、異なる a の選択が必要になることがあります。
●N=15 のような小さな数では、古典的なアルゴリズムの方がはるかに高速です。Shorのアルゴリズムの真価は、非常に大きな数を対象とした場合に発揮されます。
●実際の量子コンピュータでは、ノイズの影響により、シミュレータ通りの結果が得られないことがあります。誤り訂正符号などの技術が重要になります。

演習課題 (オプション):

●a の値を他のもの (例: 2, 8, 11, 13) に変更して実行し、結果がどう変わるか観察してみましょう。
●n_counting_qubits の数を変更して、周期推定の精度にどのような影響があるか調べてみましょう。
●N=21 の因数分解に挑戦してみましょう (ヒント: c_amodN 関数を一般化する必要があります)。