はじめに

みなさん、こんにちは。IBM Bob + IBM Q の第二弾です。いやー、Bobのおかげでやってみたかったことが、次々とチャレンジできています。
今回は、量子コンピュータの用途で代表格といえる素因数分解をやってみます。

ちなみに、Bobにプログラムにかかわるところの記事まで書いてもらいました。

本記事では、小さな例として N = 15 を IBM Q 上で実際に素因数分解します。

実行環境: Python 3.9 以上、Qiskit 1.x、qiskit-aer、qiskit-ibm-runtime

ソースコード

GitHubにアップロードしています。ご自由にダウンロードしてください。

shor-factor15-ibmq/
├── shor_factor15.py      # メインプログラム
├── requirements.txt      # 依存ライブラリ
└── README.md             # セットアップ手順

1. 素因数分解のポイント - Shorのアルゴリズム

ショアのアルゴリズムは、解くのに何億年もかかる「巨大な数の素因数分解」を、「波の周期（繰り返しのリズム）を見つける問題」に変換して量子コンピュータで一瞬で解く手法です。

15を「3 × 5」に分解する流れを、4つのステップで追ってみます。

①. 適当な数を選ぶ

まず、15より小さくて1以外の数をランダムに1つ選びます。ここでは 7 を選んだとします。

※. ショアのアルゴリズムでは、15と「共通の約数（1以外の割り切れる数）」を持たない数なら、どれを選んでも成り立ちます。

②. 7を何乗かして「15で割った余り」を並べる

7を累乗していくと、あまりの数字に決まったパターンが現れます。
$7^1 = 7 \quad \rightarrow$ 15で割った余りは 7
$7^2 = 49 \quad \rightarrow$ 15で割った余りは 4
$7^3 = 343 \quad \rightarrow$ 15で割った余りは 13
$7^4 = 2401 \quad \rightarrow$ 15で割った余りは 1
$7^5 = 16807 \quad \rightarrow$ 15で割った余りは 7 （ここで最初に戻る）

$r$	$7^r \bmod 15$
0	1
1	7
2	4
3	13
4	1 ← 周期 $r=4$

余りは「7, 4, 13, 1」のパターンを永遠に繰り返します。この「4回ごとに元に戻る」という周期 $r = 4$ を見つけることが最大の鍵です。

③. 量子コンピュータで周期を一瞬で見つける

桁数が巨大になると、普通のコンピュータは順番に計算するため途方もない時間がかかります。これが RSA 暗号の安全性の根拠です。
ショアのアルゴリズムでは、量子重ね合わせを使って「すべての計算を同時に行い」、波の打ち消し合いと強め合い（量子フーリエ変換）を利用して周期 $r = 4$ だけを浮き立たせて一瞬で検出します。

④. 周期から素数（答え）を導く

周期 $r = 4$ さえ分かれば、あとは普通の算数です。
● 選んだ数（7）の「周期の半分乗」を計算： $7^{(4/2)} = 7^2 = 49$
● 49に「1を足した数（50）」と「1を引いた数（48）」を作る ※1
● 50と15の最大公約数 $\rightarrow$ 5
● 48と15の最大公約数 $\rightarrow$ 3

これで素因数である 3 と 5 が見つかりました。計算が困難な「素因数分解」を、量子計算が得意とする「規則的な周期探し」にすることが、このアルゴリズムのポイントです。

※1の補足

$x^2 - 1 = (x - 1)(x + 1)$ という中学校で習うこの公式の形に強引に当てはめるために「周期の半分乗」と「$\pm 1$」が必要になります。

A. 「あまりが 1」を式で表す
周期 $r$ が見つかったということは、「$a^r$ を 15 で割るとあまりが 1 になる」ということです。これを引き算に直すと、$a^r - 1$ は必ず 15 の倍数になるという性質が成り立ちます。

※. 数学では「何乗もしていく計算」のスタート地点を0乗と考えます。
0乗（スタート）：どんな数も0乗は $1$ です（あまりも 1）。
あまりが「1」に戻ること自体を「1周した（周期）」と定義しているからです。

B. 因数分解の公式を使う
ここで、$a^r$ を「周期の半分乗の2乗」つまり $(a^{r/2})^2$ と言い換えます。すると、あの公式がそのまま使えます。$$a^r - 1 = (a^{r/2})^2 - 1^2 = (a^{r/2} - 1)(a^{r/2} + 1)$$

2. 量子回路の構造

この項目は、Bobにお願いしましたが、うーん。。。。。勉強しましょう！！

Shor のアルゴリズムの量子部分は 量子位相推定 (QPE) です。

2-1. 使用するゲート一覧

ゲート	記号	行列表現	役割
アダマール	H	$\frac{1}{\sqrt{2}}\begin{pmatrix}1&1\1&-1\end{pmatrix}$	重ね合わせ状態を生成
パウリ X	X	$\begin{pmatrix}0&1\1&0\end{pmatrix}$	ビット反転（NOT）
CNOT	CX	$\begin{pmatrix}1&0&0&0\0&1&0&0\0&0&0&1\0&0&1&0\end{pmatrix}$	制御 NOT
CSWAP	Fredkin	制御 SWAP	乗算ユニタリの実装
制御位相	CP($\theta$)	$\begin{pmatrix}1&0&0&0\0&1&0&0\0&0&1&0\0&0&0&e^{i\theta}\end{pmatrix}$	QFT の中核
量子フーリエ逆変換	IQFT	$F_N^\dagger$	位相を測定値に変換

2-2. count レジスタ（位相推定レジスタ）

$n$ ビットの count レジスタは、まず全ビットにアダマールゲートを適用して

$$|0\rangle^{\otimes n} \xrightarrow{H^{\otimes n}} \frac{1}{\sqrt{2^n}} \sum_{x=0}^{2^n-1} |x\rangle$$

の重ね合わせ状態を作ります。

はい！最終的に周期のヒントとなる「数字」を読み取るためのメモ帳（兼・分度器）です。

2-3. target レジスタ（作業レジスタ）

4 量子ビットで $|0\rangle$ から $|15\rangle$ を表現します。
$f(x) = 7^x \bmod 15$ の固有ベクトルを近似するため、$|1\rangle = |0001\rangle$ に初期化します。

qc.x(qr_target[0])  # |0000⟩ → |0001⟩

はい！「7を掛けて15で割る」という算数を実際に実行する電卓です。

2-4. 制御ユニタリ

$$U|y\rangle = |ay \bmod N\rangle$$

count レジスタの第 $k$ ビットを制御ビットとして、$U^{2^k}$ を target レジスタに適用します。

$U^{2^k}$ は $a^{2^k} \bmod N$ の乗算を実装します：

$k$	$7^{2^k} \bmod 15$	操作
0	7	CSWAP + CX の組合せ
1	4	CSWAP ×2
2	1	恒等（何もしない）
3	1	恒等
⋮	1	恒等

はい！、「Countレジスタのコインが表だったら、Targetレジスタで7を掛ける」という連動型の仕掛けです。この計算を行うと、 Countレジスタ側に「掛け算のリズム（波の角度のズレ）」が巻き戻るように記録されます。

2-5. 量子フーリエ逆変換 (IQFT)

IQFT は QFT のユニタリ共役 $F_N^\dagger$ です。
count レジスタに蓄積された位相情報を計算基底の振幅に変換し、測定可能にします。

$$\text{IQFT}|j\rangle = \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1} e^{-2\pi i jk/N}|k\rangle$$

n ビットの IQFT は制御位相ゲートとアダマールで分解されます：

q[n-1] ─── H ─── CP(-π/2) ─── CP(-π/4) ─── ⋯ ─── SWAP
q[n-2] ─────────────●──────────────────────
q[n-3] ────────────────────────●───────────
  ⋮

はい！　回路の中で最も重要なパーツです。あちこちを向いて混ざり合った「波の角度」の中から、強め合う波（リズム）だけを残し、「周期を表すひとつの数字」へと一瞬でまとめ上げます。

3. セットアップ

3-1. ライブラリのインストール

pip install qiskit qiskit-aer qiskit-ibm-runtime matplotlib pylatexenc numpy

3-2. IBM Quantum のアカウント設定

IBM Quantum でアカウントを作成
マイページで API トークンをコピー
環境変数または引数として渡す

# 環境変数で設定する場合
export IBM_QUANTUM_TOKEN="your_token_here"

4. 実装コードの解説

4-1. 制御ユニタリゲートの実装

def c_amod15(a: int, power: int) -> QuantumCircuit:
    """
    制御ユニタリ: |y⟩ → |(a^power) * y mod 15⟩
    q[0] = 制御ビット, q[1..4] = 4 ビット target
    """
    U_val = pow(a, power, 15)   # a^power mod 15 を事前計算
    qc = QuantumCircuit(5)
    ctrl, t = 0, [1, 2, 3, 4]

    if U_val == 4:
        # 4 * y mod 15: 1→4→1, 2→8→2 のサイクル
        qc.cswap(ctrl, t[0], t[2])   # ビット 0 ↔ ビット 2
        qc.cswap(ctrl, t[1], t[3])   # ビット 1 ↔ ビット 3
    elif U_val == 7:
        # 7 * y mod 15: 1→7→4→13→1 のサイクル
        qc.cswap(ctrl, t[1], t[3])
        qc.cswap(ctrl, t[0], t[2])
        qc.cswap(ctrl, t[0], t[3])
        qc.cx(ctrl, t[0]); qc.cx(ctrl, t[1])
        qc.cx(ctrl, t[2]); qc.cx(ctrl, t[3])
    # ... 他の U_val も同様
    return qc

ポイント: $U^{2^k}$ は $a^{2^k} \bmod N$ を予め計算してしまえば、有限個の CSWAP/CX の組合せで表現できます。

4-2. QPE 回路の構築

def build_qpe_circuit(n_count: int = 8, a: int = 7) -> QuantumCircuit:
    qr_count  = QuantumRegister(n_count, name="count")   # 位相推定レジスタ
    qr_target = QuantumRegister(4,       name="target")  # 作業レジスタ
    cr        = ClassicalRegister(n_count, name="meas")
    qc = QuantumCircuit(qr_count, qr_target, cr)

    # ① count レジスタをアダマールで重ね合わせ
    qc.h(qr_count)
    # ② target レジスタを |1⟩ に初期化
    qc.x(qr_target[0])
    qc.barrier()

    # ③ 制御 U^(2^k) を適用
    for k in range(n_count):
        U_gate = c_amod15(a, 2**k).to_gate(label=f"U^{2**k}")
        qc.append(U_gate, [qr_count[k]] + list(qr_target))
    qc.barrier()

    # ④ IQFT を適用
    iqft = QFT(n_count, inverse=True, do_swaps=True).decompose()
    qc.append(iqft, qr_count)
    qc.barrier()

    # ⑤ 測定
    qc.measure(qr_count, cr)
    return qc

4-3. 周期の復元（連分数展開）

測定値 $m$ は位相 $\varphi \approx s/r$ ($s$ は整数) を近似します。

$$\varphi = \frac{m}{2^n}$$

Fraction クラスの limit_denominator(N) メソッドで連分数展開し、分母から $r$ を推定します：

from fractions import Fraction

phase = measured_int / (2 ** n_count)
frac  = Fraction(phase).limit_denominator(N)
r     = frac.denominator   # 推定された周期

5. 実行方法

5-1. シミュレーターで実行

# デフォルト: n_count=8, shots=2048, a=7
python shor_factor15.py --mode sim

# shots を増やして精度向上
python shor_factor15.py --mode sim --shots 8192

# 回路図も保存する
python shor_factor15.py --mode sim --diagram

実行例:

============================================================
  [シミュレーター] N=15, a=7, n_count=8
============================================================

回路サイズ: 12 量子ビット, 深さ=142, ゲート数=316

測定回数: 2048 shots

--- 素因数導出 ---
  測定値=  64  位相=0.2500  r=4  gcd(7^2-1,15)=3  gcd(7^2+1,15)=5  → 15 = 3 × 5  (確率 23.8%)
  測定値= 192  位相=0.7500  r=4  gcd(7^2-1,15)=3  gcd(7^2+1,15)=5  → 15 = 3 × 5  (確率 22.9%)
  測定値= 128  位相=0.5000  r=2  gcd(7^1-1,15)=3  gcd(7^1+1,15)=1  → ... (r=2 は非自明解なし)
  測定値=   0  位相=0.0000  r=1  → スキップ

✓ 素因数分解成功: 15 = 3 × 5

ヒストグラム保存: sim_histogram.png

測定値が 64 や 192（= 256/4 × 1、256/4 × 3）に集中しているのは、位相 $\varphi = 1/4, 3/4$ に対応し、分母 4 が周期 $r = 4$ そのものだからです。

5-2. 実機で実行

# IBM Quantum API トークンを指定
python shor_factor15.py --mode real --token "YOUR_IBM_QUANTUM_TOKEN"

# バックエンドを指定する場合
python shor_factor15.py --mode real --token "TOKEN" --backend ibm_kyoto

# シミュレーターと実機を両方実行
python shor_factor15.py --mode both --token "TOKEN"

注意: 実機実行では n_count=4 に自動制限されます。量子ビット数とデコヒーレンスの制約のためです。待機時間によっては数分〜数十分かかります。

5-3. シミュレーターと実機の結果比較

項目	シミュレーター	実機
n_count	8	4
shots	2048	1024
理想的な測定値	64, 192 (n=8)	4, 12 (n=4)
ノイズ	なし	あり（NISQ）
成功確率	≈ 100%	50〜80%
実行時間	数秒	数分〜

実機では量子ビットのデコヒーレンスやゲートエラーにより、ノイズが乗った測定値が得られます。それでも最頻値付近に正しい位相に対応する測定値が現れるため、素因数分解に成功することが多いです。

実機での実行結果


使用バックエンド: ibm_marrakesh
量子ビット数  : 156
待機ジョブ数  : 0

トランスパイル後: 深さ=251, ゲート数=413

ジョブ送信完了: job_id = da7p7trsq5js73bk7u20
結果を待機中... (数分かかる場合があります)

測定回数: 1024 shots

--- 素因数導出 ---
  測定値=  12  位相=0.7500  r=4  gcd(7^2-1,15)=3  gcd(7^2+1,15)=5  -> 15 = 3 x 5  (確率 18.07%)
  測定値=   4  位相=0.2500  r=4  gcd(7^2-1,15)=3  gcd(7^2+1,15)=5  -> 15 = 3 x 5  (確率 17.58%)

[OK] 素因数分解成功: 15 = 3 x 5

Histogram saved: real_histogram.png

IBM Quantum Platformのジョブ実行結果

測定値が 4, 8, 12 という4の倍数の位置に集中していることが周期 r=4 を示しています。0100(=4) と 1100(=12) から連分数展開すると分母が4になり、これが周期 r=4 の証拠です。

その他の小さなバー（0001, 0010 など）は量子ノイズ（実機のゲートエラー・デコヒーレンス）によるものです。

6. 結果の見方

測定ヒストグラムの解釈

sim_histogram.png を見ると、いくつかの整数値に確率が集中しています。

$n = 8$ ビットの場合、理論上の測定値は：

$$m = \frac{s}{r} \times 2^n = \frac{s}{4} \times 256 \in {0, 64, 128, 192}$$

$m = 0$: $s = 0$ → 位相 0 → 周期が求まらない
$m = 64$: $s = 1$ → 位相 $1/4$ → $r = 4$ ✓
$m = 128$: $s = 2$ → 位相 $1/2$ → $r = 2$（有効な周期だが因数が自明になることがある）
$m = 192$: $s = 3$ → 位相 $3/4$ → $r = 4$ ✓

7. Shor アルゴリズムの計算量

ステップ	量子ビット数	ゲート数
QPE (count レジスタ)	$O(\log N)$	$O((\log N)^2)$
制御 $U^{2^k}$ (target)	$O(\log N)$	$O((\log N)^3)$ per gate
IQFT	$O(\log N)$	$O((\log N)^2)$
合計	$O(\log N)$	$O((\log N)^3)$

$N = 15$ の場合、理論上は $\log_2 15 \approx 4$ 量子ビットで足りますが、十分な精度を得るために count レジスタを 8 ビット以上にすることが一般的です。

8. まとめ

本記事では前回からかなりレベルアップしまして、IBM Bobに、IBM Q を使った $15 = 3 \times 5$ を素因数分解する Shor のアルゴリズムを実装してもらいました。

量子位相推定 (QPE) で $f(x) = 7^x \bmod 15$ の周期 $r = 4$ を確率的に発見
連分数展開で測定値から周期を正確に復元
GCD 計算で非自明な因数を導出

Shorのアルゴリズムや量子コンピュータの勉強は簡単にはしていましたが、実装となるとほぼ無理でした。
とういか、勉強のほうも完全に理解しているわけではないので、量子回路をどうやったら記述できるのかもわかっていませんでした。

Bobにお願いすると簡単に実装して、実機を使った検証が行え、素因数分解が行えましたね。

現在は、N=15 のような小さな数しか実用的に素因数分解できませんが、誤り訂正量子コンピュータが実現すれば大規模な RSA 暗号を解読できる可能性がありますね。

今回は、以上となります。

参考文献

P. W. Shor, "Algorithms for quantum computation: discrete logarithms and factoring," Proceedings 35th Annual Symposium on Foundations of Computer Science, 1994.
M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information, Cambridge University Press, 2010.
Qiskit Textbook — Shor's Algorithm
IBM Quantum Documentation

参考

SPSS Modeler ノードリファレンス目次

SPSS Modeler 逆引きストリーム集

SPSS funさん記事集

IBM 斎藤さんの記事集

SPSS連載ブログバックナンバー

SPSSヒモトクブログなどは以下のTechXchangeのコミュニティに統合されました。
ご興味がある方は、ぜひiBM IDを登録して参加してみてください！！！お待ちしています。

IBM TechXchange Data Science Japan

IBM Bob紹介ページ

IBM Quantum Platform