Python（ショアのアルゴリズム02）
■Qiskitでのショアのアルゴリズム。
Qiskitのチュートリアルで、ショアのアルゴリズムがあったのでそこからコードを組み立てて試した。その中では N=15（対象の数）, a=2 （Nと互いに素な数）を例として挙げていた。
ユニタリとか連分数展開とか数学的な内容が出てきて、何となくしか理解できていないけど、流れは下のような感じかと思う。
対象の数に対して互いに素になる数字 a をランダム指定。ここでは 2。0 からの連番として k。a^(2^k)の数字に対して15で割って余りを出す。Excelで確認したものが下（多分数が大きすぎて 、k = 6, 7 ではエラーが出てる）。チュートリアルのコード内では繰り返し2乗法というものを使っているそう。
a^(2^k) mod N の結果一覧（[2, 4, 1, 1, 1, 1, 1, 1]）から 2, 4に対してゲートを設定する。1のところはいらないよう。下の図で2mod15, 4mod15がそのゲート。
その後、逆フーリエ変換（IQFTの部分）を通して、測定。
この計算結果と制御ビット数（ここでは8）から少数を得る（例えば、11000000 = 192, 2^8 = 256 で192/256）。それを分数近似する（連分数展開）。公式の θ = k / r の考えから r (位相) を得る（あくまで候補として）。ここまでの計算で素因数を求める計算（a^(r/2) ± 1 と N の最大公約数）に必要な値が求められる。ここでは a = 2, r = 4 のため、a^(r/2) ± 1 = 3, 5になり、15の要素として正しい。
N=15, a=2の実行結果が下。
Register Output Phase
0 11000000(bin) = 192(dec) 192/256 = 0.75
1 10000000(bin) = 128(dec) 128/256 = 0.50
2 00000000(bin) = 0(dec) 0/256 = 0.00
3 01000000(bin) = 64(dec) 64/256 = 0.25
Phase Fraction Guess for r
0 0.75 3/4 4
1 0.50 1/2 2
2 0.00 0/1 1
3 0.25 1/4 4
{'01000000': 249, '11000000': 257, '10000000': 259, '00000000': 259}
ATTEMPT 0:
Phase: theta = 0.25
Order of 2 modulo 15 estimated as: r = 4
*** Non-trivial factor found: 3 ***
前回試した N = 26, a = 17 で実行した結果が下。
最後は24と出てるので、なんかおかしい所がありそう。結果のビット数から手計算で行うと r = 2, 6 となった。r = 6 は、前回の古典の方法で行ったのと同じ結果になっている。
ショアのアルゴリズムまで見たのでQiskitはここまで。上で出てきた繰り返し2乗法とか連分数展開とかどこかの機会で見たいかも。
最後に使ったコードをのせておく（Qiskitのチュートリアルのコードを基に多少修正したもの）。
ibmのruntimeが必要。動作環境は下のバージョン。
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from fractions import Fraction
from math import floor, gcd, log
from qiskit_aer import AerSimulator
from qiskit import transpile
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.circuit.library import QFT, UnitaryGate
from qiskit_ibm_runtime import SamplerV2 as Sampler
def M_amodN(a, N):
"""
Return a UnitaryGate implementing |x> -> |a*x mod N>
a: multiplier
N: modulus
"""
# 必要な qubit 数（N の範囲を表現できるだけ）
n = int(np.ceil(np.log2(N)))
dim = 2**n
# 置換行列を作成
U = np.zeros((dim, dim))
for x in range(N):
y = (a * x) % N
U[y, x] = 1
# N 以上の値はそのまま写す（不要なら 0 にしてもよい）
for x in range(N, dim):
U[x, x] = 1
# ユニタリゲートに変換
gate = UnitaryGate(U, label=f"{a}mod{N}")
return gate
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
print(f"a2kmodN:{a}, k:{k}, N:{N}, a~2:{a ** 2}")
"""################################################################################"""
N = 26
a = 17
# 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(a, k, N) for k in k_list]
print(b_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 != 1: # M1は恒等作用素なのでスキップ
circuit.compose(
M_amodN(b, N).control(), qubits=[qubit] + list(target), inplace=True
)
# 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)
plt.show()
# Aer シミュレーション
backend = AerSimulator()
tqc = transpile(circuit, backend)
result = backend.run(tqc, shots=1024).result()
counts = result.get_counts()
"""################################################################################"""
# 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)
# Rows to be displayed in a table
rows = []
for phase in measured_phases:
frac = Fraction(phase).limit_denominator(N)
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)
# 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
# pub = transpiled_circuit
pub = tqc
job = sampler.run([pub], shots=1024)
result = job.result()[0]
counts = result.data["out"].get_counts()
# 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)
"""################################################################################"""
FACTOR_FOUND = False
num_attempt = 0
while not FACTOR_FOUND:
print(f"\nATTEMPT {num_attempt}:")
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} ***")
