DAL Tech BlogPublicationへの投稿
➗

量子コンピュータでShorのアルゴリズムを実行する

Gen Sekita
に公開
量子コンピュータ
Qiskit
quantum
shorのアルゴリズム
tech

はじめに

データソリューション事業部の関田です。
今回は、量子アルゴリズムの1つであるShorのアルゴリズムにより、21を素因数分解します。
量子計算用ライブラリのQiskitを用いて実装し、IBMの量子コンピュータ実機で計算してみました。

本記事はDAL Tech Blog Advent Calendar 2025として投稿しました。全ての記事は以下からご確認いただけます。

https://adventar.org/calendars/12288

Shorのアルゴリズム概要

Shorのアルゴリズムは、素因数分解を行う量子アルゴリズムとして知られています。
素因数分解を行う具体的な手段として、以下の関数 f(x) の周期 r を求めます（ここでは詳細を略しますが、都合の良い r が求められたら素因数分解はほぼ完了といえます）。

f(x) = a^x \ ( \rm{mod} \ N )

ここで、N を素因数分解したい値、a を N と互いに素な整数として選びます。
例えば、この後設定する N=21, \ a=4 では、以下のようになります。

\begin{aligned} f(0) &= 1, \\ f(1) &= 4, \\ f(2) &= 16, \\ f(3) &= 1, \\ f(4) &= 4. \end{aligned}

このように 1, 4, 16 という3つの値を繰り返すため、周期は 3 です。
この周期を発見するところで量子計算を行うのですが、量子計算部分の詳しい解説はここではスキップさせていただきます。
キーワードは次章の実装コード内に記載しているので、詳細が気になる方はそちらを参考にShorのアルゴリズムについて調べていただけると幸いです。

ここで求められた周期 r と a の値によっては素因数分解できない場合がありますが、今回は都合の良い (N, a) を選択したため素因数分解することができます。

21を素因数分解する

実際に N=21, \ a=4 でShorのアルゴリズムを実行してみます。
実装ではQiskitというライブラリを用いています。
Qiskit関連ライブラリのバージョンは次の通りです。
※Qiskitはバージョンによってプログラムが大きく変わるので要注意です。

  • qiskit 2.2.3
  • qiskit-aer 0.17.2
  • qiskit-ibm-runtime 0.43.1

まず、コードは以下です。

Shorのアルゴリズム
import time
from fractions import Fraction
from math import gcd

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.circuit.library import QFT
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

# シミュレータ用
from qiskit_aer import AerSimulator

# 実機用 (IBM Quantum)
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler

# ==========================================
# 1. パラメータ設定
# ==========================================
N = 21
a = 4  # aとNは互いに素である必要があります。a=4は周期r=3を持ちます(4^1=4, 4^2=16, 4^3=64=1 mod 21)
n_target = 5  # N=21を表現するためのビット数
n_count = 6   # 精度を決めるビット数

# 実行モード切り替え (Trueにすると実機へ送信されます。要API設定)
USE_REAL_HARDWARE = False
API_TOKEN = YOUR_API_TOKEN # 実機を使う場合は設定してください

# ==========================================
# 2. 回路の構築関数
# ==========================================
def c_amod21(power):
    """
    a = 4 の場合の制御ユニタリゲート (Controlled-U^power)。
    4^x mod 21 の動作：
    1 -> 4 -> 16 -> 1 ... というサイクルを作ります。
    |00001> -> |00100> -> |10000> -> |00001>
    """
    U = QuantumCircuit(n_target)
    for _ in range(power):
        # 簡易的な置換ゲートとして実装 (SWAP論理でサイクルを作る)
        # Note: ここではa=4に特化した実装を行っています
        # 1(00001) -> 4(00100)
        # 4(00100) -> 16(10000)
        # 16(10000) -> 1(00001)
        # これを実現する簡易SWAPロジック
        U.swap(0, 2)
        U.swap(2, 4)

    U_gate = U.to_gate()
    U_gate.name = f"{a}^{power} mod 21"
    return U_gate.control()

def create_shor_circuit(n_count, n_target):
    # レジスタの作成
    c = QuantumRegister(n_count, 'count')
    t = QuantumRegister(n_target, 'target')
    cr = ClassicalRegister(n_count, 'outcome')
    qc = QuantumCircuit(c, t, cr)

    # 1. カウントレジスタを重ね合わせ状態に初期化
    qc.h(c)

    # 2. ターゲットレジスタを |1> に初期化
    qc.x(t[0])

    # 3. 制御モジュラー指数計算 (Controlled-U)
    for i in range(n_count):
        # 2^i 回の繰り返し適用
        qc.append(c_amod21(2**i), [c[i]] + list(t))

    # 4. 逆量子フーリエ変換 (IQFT)
    qc.append(QFT(n_count, inverse=True).to_gate(), c)

    # 5. 測定
    qc.measure(c, cr)
    return qc

# ==========================================
# 3. メイン実行処理
# ==========================================
print(f"Shorのアルゴリズム: N={N}, a={a}")

# 回路作成
qc = create_shor_circuit(n_count, n_target)

# --- 情報出力: 回路図 ---
qc.draw(output='mpl', filename=OUTPUT_PATH, fold=-1)

# --- 実行準備 ---
if USE_REAL_HARDWARE:
    service = QiskitRuntimeService(channel="ibm_quantum_platform", token=API_TOKEN)
    backend = service.least_busy(operational=True, simulator=False)
    print(f"\n実機での実行: {backend.name}")
else:
    # 実機を模したシミュレータ、または完全な理想シミュレータを使用
    backend = AerSimulator()
    print(f"\nシミュレータでの実行: {backend.name}")

# --- 計算実行と時間計測 ---
start_time = time.time()

# Sampler V2 による実行
sampler = Sampler(mode=backend)
pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
isa_circuit = pm.run(qc)

job = sampler.run([isa_circuit], shots=2048)
result = job.result()

end_time = time.time()
elapsed_time = end_time - start_time

# --- 結果処理 ---
pub_result = result[0]
counts = pub_result.data.outcome.get_counts()
sorted_counts = sorted(counts.items(), key=lambda x: x[1], reverse=True)

print(f"\n実行時間: {elapsed_time:.4f} 秒")
print(f"測定回数: {sorted_counts}")

# ==========================================
# 4. 後処理（位相から周期を推定）
# ==========================================
print("\n位相から周期を推定")
print("測定回数が多い3つの結果を表示")

for output_bitstring, count in sorted_counts[:3]: # 上位3つの結果を表示
    # 10進数に変換
    decimal_measured = int(output_bitstring, 2)
    phase = decimal_measured / (2**n_count)

    # 連分数展開で近似分数を見つける
    frac = Fraction(phase).limit_denominator(N)
    r_guess = frac.denominator

    print("==========================================")
    print(f"測定された状態: {output_bitstring} (10進数: {decimal_measured})")
    print(f" -> 位相: {phase:.4f}")
    print(f" -> 連分数展開で得られる分数: {frac.numerator}/{frac.denominator}")
    print(f" -> 周期 r = {r_guess}")

    if r_guess == 3:
        guess_factor1 = gcd(int(2**(r_guess) - 1), N)
        guess_factor2 = gcd(int(2**(r_guess) + 1), N)
        print(f"素因数分解成功！ {N} = {guess_factor1} × {guess_factor2}")
    else:
        print("周期が間違っています")
print("==========================================")

まずは、このコードをそのまま実行してみます。すると、以下の出力が得られます。

出力1
Shorのアルゴリズム: N=21, a=4

シミュレータでの実行: aer_simulator

実行時間: 0.1071 秒
測定回数: [('000000', 633), ('101011', 502), ('010101', 470), ('010110', 141), ('101010', 117), ('101100', 33), ('010100', 32), ('010111', 22), ('101001', 19), ('010011', 16), ('101101', 9), ('101000', 6), ('100111', 5), ('100110', 3), ('011000', 3), ('011100', 3), ('010010', 3), ('100000', 3), ('101110', 2), ('011011', 2), ('010001', 2), ('001110', 2), ('100101', 2), ('110000', 2), ('101111', 2), ('111011', 2), ('110101', 1), ('010000', 1), ('001111', 1), ('011010', 1), ('110110', 1), ('011101', 1), ('100010', 1), ('001100', 1), ('110001', 1), ('011001', 1), ('111101', 1), ('011110', 1)]

位相から周期を推定
測定回数が多い3つの結果を表示
==========================================
測定された状態: 000000 (10進数: 0)
 -> 位相: 0.0000
 -> 連分数展開で得られる分数: 0/1
 -> 周期 r = 1
周期が間違っています
==========================================
測定された状態: 101011 (10進数: 43)
 -> 位相: 0.6719
 -> 連分数展開で得られる分数: 2/3
 -> 周期 r = 3
素因数分解成功！ 21 = 7 × 3
==========================================
測定された状態: 010101 (10進数: 21)
 -> 位相: 0.3281
 -> 連分数展開で得られる分数: 1/3
 -> 周期 r = 3
素因数分解成功！ 21 = 7 × 3
==========================================

==で囲まれた部分に着目してください。
3パターンの結果を出力していますが、そのうち2つで周期を正しく求め、素因数分解することができています。
※複数パターン出力している理由は、量子の性質によりどうしても不確実性が生じてしまったり、位相0というあまり意味を持たない解が出てくることにより、1パターンのみだと失敗してしまう可能性があるためです。
また、上記はシミュレータ（※量子コンピュータは用いていない）での実行結果ですが、実行時間は約0.1秒と瞬時に計算できていることもわかります。

パラメータ変更

コードの中で、以下の2つのパラメータを変更して挙動を見てみます。

  • n_count: 結果を測定するための量子ビットの数で、精度に関わります。
  • USE_REAL_HARDWARE: 量子コンピュータの実機を用いるか否かを決めます。

測定用量子ビット数

測定用量子ビット数を変更すると、結果がどのように変化するかを見てみます。
出力1では、測定用量子ビット数 (n_count) を6と設定していました。
これを3に設定してみます。測定用量子ビット数を減らすと、測定の精度が下がり、所望する結果を得ることができない可能性があります。
結果は以下の通りです。

出力2
Shorのアルゴリズム: N=21, a=4

シミュレータでの実行: aer_simulator

実行時間: 0.1000 秒
測定回数: [('000', 737), ('011', 493), ('101', 459), ('110', 112), ('010', 106), ('100', 67), ('001', 40), ('111', 34)]

位相から周期を推定
測定回数が多い3つの結果を表示
==========================================
測定された状態: 000 (10進数: 0)
 -> 位相: 0.0000
 -> 連分数展開で得られる分数: 0/1
 -> 周期 r = 1
周期が間違っています
==========================================
測定された状態: 011 (10進数: 3)
 -> 位相: 0.3750
 -> 連分数展開で得られる分数: 3/8
 -> 周期 r = 8
周期が間違っています
==========================================
測定された状態: 101 (10進数: 5)
 -> 位相: 0.6250
 -> 連分数展開で得られる分数: 5/8
 -> 周期 r = 8
周期が間違っています
==========================================

測定用量子ビット数を3に減らすと、正しい周期を得ることができませんでした。このように、量子ビット数が不足していると正しい結果を得ることが難しくなります。

量子コンピュータ実機での実行

次に、量子コンピュータの実機を用いて計算してみます。
出力1を得たコードで USE_REAL_HARDWARE = True にして実行します（測定用量子ビット数は6）。
結果は以下の通りです。

出力3
Shorのアルゴリズム: N=21, a=4

実機での実行: ibm_torino

実行時間: 8.8505 秒
測定回数: [('001011', 107), ('101111', 95), ('001010', 93), ('001111', 91), ('101011', 81), ('101101', 77), ('101110', 75), ('001101', 74), ('001110', 72), ('001100', 68), ('101010', 64), ('001001', 60), ('101001', 56), ('101100', 54), ('011000', 52), ('011001', 51), ('001000', 43), ('111001', 42), ('000110', 40), ('111011', 36), ('011011', 36), ('101000', 33), ('011010', 31), ('111000', 31), ('000100', 27), ('100110', 27), ('000011', 26), ('000111', 25), ('000101', 24), ('000010', 24), ('100100', 22), ('100111', 22), ('100101', 21), ('011111', 21), ('010001', 19), ('100010', 18), ('111010', 18), ('100000', 16), ('110011', 16), ('100011', 15), ('011110', 14), ('011100', 14), ('110000', 14), ('011101', 14), ('111111', 14), ('000000', 14), ('111101', 14), ('110001', 13), ('100001', 13), ('010000', 12), ('010010', 11), ('111110', 10), ('000001', 10), ('110100', 9), ('110101', 9), ('010011', 8), ('010100', 8), ('110111', 8), ('110110', 8), ('111100', 7), ('010110', 6), ('110010', 6), ('010111', 5), ('010101', 4)]

位相から周期を推定
測定回数が多い3つの結果を表示
==========================================
測定された状態: 001011 (10進数: 11)
 -> 位相: 0.1719
 -> 連分数展開で得られる分数: 3/17
 -> 周期 r = 17
周期が間違っています
==========================================
測定された状態: 101111 (10進数: 47)
 -> 位相: 0.7344
 -> 連分数展開で得られる分数: 11/15
 -> 周期 r = 15
周期が間違っています
==========================================
測定された状態: 001010 (10進数: 10)
 -> 位相: 0.1562
 -> 連分数展開で得られる分数: 3/19
 -> 周期 r = 19
周期が間違っています
==========================================

まず実行時間を見ると、実機で計算すると8秒以上かかっており、シミュレータよりも遅くなっています。
実行時間の内訳は次のようになっています。

  • 私が量子コンピュータに投入したjobが実行されるまでに生じる待ち時間: 1秒
  • 量子コンピュータでの計算: 3秒
  • その他通信などのオーバーヘッド: 4秒

実行時間などの詳細は、IBM Quantum Platformの各ユーザーアカウント内で閲覧することができます（要登録）。

続いて、得られた周期を見ると、正しい周期を求めることができていません。正しい結果が得られなかった主な原因は、量子コンピュータの性能にあると考えられます。
量子コンピュータは計算の過程で様々な原因でノイズが発生し、デタラメなふるまいをしてしまいます。これを制御するためには多くの量子ビットが必要とされていますが、現状の量子コンピュータにはエラーを訂正できるほどの量子ビットがありません。
ノイズの影響をなるべく軽減する他の手段として、量子回路を浅くする、ということが可能です。つまり、エラーの発生箇所である「演算」の数を少なくするという発想です。これはQiskit上で簡単に試すことができ、generate_preset_pass_manager(backend=backend, optimization_level=1)の optimization_level を3などに上げることで、自動で演算を少なくするように工夫してくれます。
generate_preset_pass_manager(backend=backend, optimization_level=3)として実機で実行すると、なんとか正しい周期が3番目に出てきました。
ただし正しい状態が高確率で得られているわけではないので、うまく計算できているとはいえません。

出力4
Shorのアルゴリズム: N=21, a=4

実機での実行: ibm_fez

実行時間: 11.9768 秒
測定回数: [('110000', 52), ('100000', 51), ('101011', 50), ('100001', 49), ('110010', 48), ('010000', 43), ('111010', 42), ('100011', 42), ('101001', 41), ('101110', 40), ('111000', 40), ('101010', 39), ('011010', 38), ('001001', 38), ('110011', 38), ('011001', 38), ('100010', 37), ('011000', 37), ('100100', 37), ('101100', 36), ('011011', 36), ('001110', 35), ('110001', 35), ('101111', 35), ('101000', 35), ('010100', 34), ('001000', 34), ('000000', 33), ('000110', 32), ('001010', 31), ('111110', 31), ('100101', 31), ('100110', 31), ('001100', 31), ('010101', 31), ('010010', 30), ('000001', 30), ('010001', 30), ('110110', 29), ('001011', 29), ('000011', 29), ('110100', 29), ('111100', 28), ('110111', 28), ('011111', 27), ('001111', 26), ('011100', 25), ('111001', 25), ('000101', 25), ('000100', 25), ('010111', 25), ('100111', 25), ('110101', 24), ('011110', 24), ('011101', 23), ('000010', 23), ('010011', 23), ('101101', 23), ('010110', 22), ('001101', 21), ('111011', 20), ('111101', 18), ('111111', 17), ('000111', 14)]

位相から周期を推定
測定回数が多い3つの結果を表示
==========================================
測定された状態: 110000 (10進数: 48)
 -> 位相: 0.7500
 -> 連分数展開で得られる分数: 3/4
 -> 周期 r = 4
周期が間違っています
==========================================
測定された状態: 100000 (10進数: 32)
 -> 位相: 0.5000
 -> 連分数展開で得られる分数: 1/2
 -> 周期 r = 2
周期が間違っています
==========================================
測定された状態: 101011 (10進数: 43)
 -> 位相: 0.6719
 -> 連分数展開で得られる分数: 2/3
 -> 周期 r = 3
素因数分解成功！ 21 = 7 × 3
==========================================

おわりに

今回は21を素因数分解するShorのアルゴリズムを実装し、量子コンピュータ実機とシミュレータで計算してみました。
シミュレータで正しく動作することを確認した後、同様に量子コンピュータでも動かしてみましたが、ハードウェアの性能の問題で正確な結果を得ることができませんでした。
個人的には予想通りの結果でしたが、実際に手を動かしてみるとアルゴリズムやライブラリの使い方の理解が深まり、より一層面白く感じました！
Shorのアルゴリズムや量子コンピューティングに興味がある方は、本記事などを参考に是非手を動かしてみてください！

DAL Tech Blog
DAL Tech BlogPublication

データアナリティクスラボ株式会社（DAL）のテックブログです！　会社情報、採用情報、JOURNAL（専門的な技術記事）、オウンドメディアは以下リンクからご覧ください。

Discussion