Pythonでのショアの素因数分解アルゴリズムの実装

この記事では、量子コンピュータを用いて素因数分解を行うショアのアルゴリズム（Shor's algorithm）をPythonで実装する方法を解説します。Qiskitライブラリを使用し、実際のコード例とともにアルゴリズムの動作を説明します。

ショアのアルゴリズムとは
- 古典アルゴリズムとの比較
アルゴリズムの仕組み
- 量子フーリエ変換 (QFT)
具体例: 15の素因数分解
Python (Qiskit) での実装
まとめ

ショアのアルゴリズムとは

ショアのアルゴリズムは、1994年にPeter Shorによって考案された、量子コンピュータ上で効率的に素因数分解を行うためのアルゴリズムです。従来のコンピュータ（古典コンピュータ）では困難な大きな数の素因数分解を、量子コンピュータを用いることで高速に解ける可能性があります。これは、RSA暗号などの公開鍵暗号システムの安全性を脅かす可能性があるなど、非常に重要なアルゴリズムです。

古典アルゴリズムとの比較

試行除算: 素因数分解の最も単純な方法は、小さい数から順に割っていく方法です。しかし、非常に効率が悪く、大きな数に対しては現実的な時間で解けません。
より高度なアルゴリズム: 数体ふるい法など、より効率的な古典アルゴリズムも存在しますが、それでも計算量は指数時間であり、非常に大きな数の素因数分解は困難です。

ショアのアルゴリズムの優位性

ショアのアルゴリズムは、量子コンピュータ上で多項式時間で素因数分解を解くことができます。具体的には、数のビット数を n としたとき、O(n³)程度の計算量で解けるとされています。これは、古典コンピュータにおける最良のアルゴリズムよりもはるかに高速です。

アルゴリズムの仕組み

ショアのアルゴリズムは、以下の2つの主要部分から構成されます。

古典コンピュータによる前処理
- 素因数分解したい数 N に対し、位数 r を求める問題を解くことに帰着させます。
  - 位数とは: ある数 a に対して、a^r ≡ 1 (mod N)を満たす最小の正の整数 r のこと。
- この位数 r から、N の非自明な約数（1と N 以外の約数）を効率的に見つけられます。
量子コンピュータによる位数発見
- 量子フーリエ変換を利用して、周期関数から位数 r を効率的に求めます。

量子フーリエ変換 (QFT)

量子フーリエ変換は、ショアのアルゴリズムにおいて中心的な役割を果たす、量子コンピュータ特有の演算ですので少し詳しく解説します。これは、古典コンピュータにおける高速フーリエ変換（FFT）に相当するもので、特定の周期性を持つ状態の周期を効率的に見つけ出すことができます。

量子フーリエ変換 (QFT) の働き

QFTは、例えるなら、音の波形を周波数成分に分解するような働きをします。音波に様々な高さの音（周波数）が含まれているように、量子状態にも様々な「周期」が重ね合わされた状態で含まれていることがあります。QFTは、この重ね合わされた状態を、「周期」ごとに分解する役割を果たします。

具体的には、ある周期 r を持つ量子状態に対してQFTを適用すると、その状態は、1/rの倍数の位置にピークを持つ状態に変換されます。つまり、QFTの出力状態を測定することで、元の状態の周期 r に対応する情報を得ることができるのです。

逆QFT

QFTの逆変換（逆QFT）は、QFTで周波数成分に分解された状態を、元の状態に戻す操作です。ショアのアルゴリズムでは、位数発見の最後に逆QFTを適用することで、「周期 r に対応するピーク」の位置を明確にし、測定によってその値を読み出します。

具体例: 15の素因数分解

ここでは、N=15 を例に、アルゴリズムの手順を説明します。（はてなの数式がうまく書けませんでした、すみません。）

前処理:
- N=15 とする。
- a を N と互いに素な数（1 < a < N）としてランダムに選ぶ。ここでは、a=7 を選ぶ。
- a^r ≡ 1 (mod N)を満たす最小の正の整数 r（位数）を求める。この r を求める部分が、量子コンピュータの得意とするところです。
- r が偶数であれば、x ≡ a^(r/2) (mod N)を計算する。
- x+1 または x-1 と N の最大公約数(GCD)を計算すると、N の非自明な約数が求まる。
量子コンピュータによる位数発見:
- 量子フーリエ変換を用いて、a^x mod N の周期 r を見つけます。
後処理:
- r が求まったら、gcd(a^(r/2) - 1, N)とgcd(a^(r/2) + 1, N)を計算します。
- この例では、r=4 が得られるので、gcd(7^(4/2) - 1, 15) = gcd(48, 15) = 3 および gcd(7^(4/2) + 1, 15) = gcd(50, 15) = 5 となり、15の素因数3と5が得られます。

Python (Qiskit) での実装

量子コンピュータはまだ一般的ではありませんが、Qiskitのような量子プログラミングフレームワークを用いることで、古典コンピュータ上に量子コンピュータのシミュレータを構築し、量子アルゴリズムの動作を検証できます。これにより、手元のPCでも量子アルゴリズムの学習や実験が可能です。Qiskitを用いて、ショアのアルゴリズムを実装していきます。

1. ライブラリのインポート

ここでは、必要なライブラリをインポートします。QuantumCircuit, Aer, execute, plot_histogramはQiskitの基本的な機能を利用するためのモジュールです。mathは、最大公約数(GCD)の計算など、数学関数を使用するためのモジュールです。

import math
from qiskit import QuantumCircuit, Aer, execute
from qiskit.visualization import plot_histogram
import pandas as pd
from fractions import Fraction

2. 古典的前処理

ここでは、素因数分解する数 N と、N と互いに素な数 a を設定します。古典的アルゴリズムで位数 r を求め、求めた r を用いて、N の約数を計算できるか試します。

N = 15
a = 7  # Nと互いに素な数

# 古典的な方法で位数rを求める（確認用）
def classical_order_finding(a, N):
    r = 1
    while (a**r) % N != 1:
        r += 1
    return r

r = classical_order_finding(a, N)
print(f"Classical order finding: r = {r}")

# 約数の計算
if r % 2 == 0:
    x = a**(r // 2) % N
    factor1 = math.gcd(x - 1, N)
    factor2 = math.gcd(x + 1, N)
    print(f"Factors: {factor1}, {factor2}")
else:
    print("Order is odd. Cannot find factors.")

classical_order_finding(a,N): 古典的アルゴリズムで位数 r を求めます。
求めた r を用いて、N の約数を計算します。

3. 量子回路の準備

量子コンピュータで位数 r を求める部分を作成します。ここでは、[tex:a^x \mod N]を計算する量子回路（モジュラーべき乗演算）が必要です。 Qiskitでは、これを直接実装する関数は提供されていないため、自分で作成する必要があります。ここでは、モジュラーべき乗演算を行う回路を制御ゲート化する関数c_amod15を定義します。

def c_amod15(a, power):
    """Controlled multiplication by a mod 15"""
    if a not in [2,7,8,11,13]:
        raise ValueError("'a' must be 2,7,8,11 or 13")
    U = QuantumCircuit(4)
    for iteration in range(power):
        if a in [2,13]:
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
        if a in [7,8]:
            U.swap(2,3)
            U.swap(1,2)
            U.swap(0,1)
        if a == 11:
            U.swap(1,3)
            U.swap(0,2)
        if a in [7,11,13]:
            for q in range(4):
                U.x(q)
    U = U.to_gate()
    U.name = "%i^%i mod 15" % (a, power)
    c_U = U.control()
    return c_U

aとpowerを受け取り、対応する制御ゲートを返します。
Uは4量子ビットの量子回路として定義されます。
aの値に応じて、swapゲートやxゲートでビットの入れ替えや反転を行います。

4. 量子フーリエ変換 (QFT)

位数発見に必要となる量子フーリエ変換(QFT)と、逆量子フーリエ変換を定義します。ここでは、n 量子ビットの逆量子フーリエ変換を行う回路を作成する関数qft_daggerを定義します。

def qft_dagger(n):
    """n-qubit QFTdagger the first n qubits in circ"""
    qc = QuantumCircuit(n)
    # Don't forget the Swaps!
    for qubit in range(n//2):
        qc.swap(qubit, n-qubit-1)
    for j in range(n):
        for m in range(j):
            qc.cp(-math.pi/float(2**(j-m)), m, j)
        qc.h(j)
    qc.name = "QFT†"
    return qc

qc.swap(qubit, n-qubit-1): 量子ビットの順序を反転させます。これはQFTの定義に現れるビット順の反転に対応します。
qc.cp(-math.pi/float(2**(j-m)), m, j): 制御位相ゲートを適用します。
qc.h(j): アダマールゲートを適用します。

5. 量子回路の構築

ここでは、位数発見のための量子回路を構築します。初期状態を設定し、重ね合わせ状態を作成し、モジュラーべき乗演算、逆量子フーリエ変換を適用し測定します。

# 量子ビット数
n_count = 8  # 位相推定に使う量子ビット数
# 量子回路
qc = QuantumCircuit(n_count + 4, n_count)

# 初期状態 |1>
qc.x(n_count)

# 重ね合わせ状態の作成
for q in range(n_count):
    qc.h(q)

# モジュラーべき乗演算
for q in range(n_count):
    qc.append(c_amod15(a, 2**q), [q] + [i+n_count for i in range(4)])

# 逆量子フーリエ変換
qc.append(qft_dagger(n_count), range(n_count))

# 測定
qc.measure(range(n_count), range(n_count))
qc.draw(output='mpl')

n_count: 位相推定に使用する量子ビット数を設定します。
量子回路qcを作成します。最初のn_count個の量子ビットにアダマールゲートを適用し、重ね合わせ状態を作成します。
モジュラーべき乗演算、逆量子フーリエ変換を適用し測定を行います。

6. 回路の実行と結果の解析

ここでは、作成した量子回路をシミュレータで実行し、測定結果を取得します。測定結果をヒストグラムで表示し、測定された位相から位数 r の候補を求めます。

backend = Aer.get_backend('qasm_simulator')
results = execute(qc, backend, shots=2048).result()
counts = results.get_counts()
plot_histogram(counts)

qasm_simulatorで回路を実行し、結果を取得します。
測定結果をヒストグラムで表示します。

rows, measured_phases = [], []
for output in counts:
    decimal = int(output, 2)  # binary to decimal
    phase = decimal/(2**n_count)  # phase
    measured_phases.append(phase)
    rows.append([f"{output}(bin) = {decimal:>3}(dec)",
                 f"{decimal}/{2**n_count} = {phase:.2f}"])
# Print the rows in a table
headers=["Register Output", "Phase"]
df = pd.DataFrame(rows, columns=headers)
print(df)

測定された位相をmeasured_phasesに格納します。
位相から、位数 r の候補を求めます。

# 測定された位相から位数rを求める
frac = Fraction(0.25).limit_denominator(15) #例
s = frac.numerator
r = frac.denominator
print(r)

測定された位相を分数で近似します。
分母が位数 r の候補となります。

7. 素因数の取得

ここでは、得られた位数 r が偶数であるかを確認し、素因数を計算します。

if r % 2 == 0:
    x = a**(r // 2) % N
    factor1 = math.gcd(x - 1, N)
    factor2 = math.gcd(x + 1, N)
    print(f"Factors: {factor1}, {factor2}")
else:
    print("Order is odd. Cannot find factors.")

得られた位数 r が偶数であれば、N の約数を計算します。
gcd(x-1, N)とgcd(x+1, N)を計算し、素因数分解の結果を得ます。

まとめ

この記事では、PythonとQiskitを用いてショアの素因数分解アルゴリズムを実装する方法を解説しました。ショアのアルゴリズムは、量子コンピュータの能力を示す重要なアルゴリズムであり、RSA暗号の解読など、暗号技術に大きな影響を与える可能性があります。ただし、現在の量子コンピュータはまだ発展途上であり、大規模な数の素因数分解を実用的に行えるようになるには、さらなる技術革新が必要です。他の量子アルゴリズム「グローバーのアルゴリズム」については下記の記事で解説しています。

pydocument.hatenablog.com

問題解決のための「アルゴリズム×数学」が基礎からしっかり身につく本 [ 米田優峻 ]