New!! 栗まんじゅうでわかるショアのアルゴリズム（①15個溜まったら捨てる）

import numpy as np
import matplotlib.pyplot as plt
from qiskit_aer import AerSimulator
from qiskit import QuantumCircuit, transpile
from qiskit.circuit.library import UnitaryGate

#コントロール8ビット + ターゲット4ビット
qc = QuantumCircuit(12)

#ターゲットビットの最下位を |1> に初期化
qc.x(8)

#コントロールビットすべてにアダマールを適用
for i in range(8):
    qc.h(i)

# --- モジュロべき乗演算 ---
#行列を用いて x -> (x * a_power) mod N の変換ゲートを作る関数
def get_mod_matrix(a_power):
    U = np.zeros((16, 16)) #ターゲットが4ビットなので 2^4
    for x in range(16):
        if x < 15:
            y = (x * a_power) % 15
            U[y, x] = 1
        else:
            U[x, x] = 1 #N以上は恒等変換
    return U

#コントロールビットすべてにゲートを適用
for i in range(8):
    #aを決め、a^(2^i) mod N を計算
    a_power = pow(2, 2**i, 15)
    U_matrix = get_mod_matrix(a_power)
    
    #行列を制御ゲートに変換
    c_U = UnitaryGate(U_matrix).control(1)
    
    #コントロールビット i から、ターゲットレジスタすべてへ適用
    qc.append(c_U, [i] + list(range(8, 12)))

qc.barrier()

# --- 逆量子フーリエ変換 (IQFT) ---
#ターゲットビットすべてに適用
for j in range(4):
    qc.swap(j, 8 - 1 - j)

#コントロールビットすべてに適用
for i in range(8):
    for m in range(i):
        qc.cp(-np.pi / (2**(i - m)), i, m)
    qc.h(i)

qc.barrier()

#計算
qc.measure_all()
backend = AerSimulator(method='matrix_product_state')
t_qc = transpile(qc, backend) #トランスパイルが必要
result = backend.run(t_qc, shots=1000).result().get_counts()

#解の確認
for r in result:
    print(f'{r[::-1]} | {result[r]}/1000')

#解のコントロールビットを集計して確認
all_pattern = {}
for r in result:
    pattern = r[::-1][:8][::-1] #上位ビットから並ぶようにひっくり返しておく
    if pattern in all_pattern:
        all_pattern += result[r]
    else:
        all_pattern = result[r]
print(all_pattern)

#ピークを確認するためのグラフ
x = range(2**8)
y = [0] * 2**8
for pattern in all_pattern:
    y[int(pattern, 2)] += all_pattern
plt.plot(x, y)
plt.xticks(range(0, 2**8 + 1, 32))
plt.show()

000000111000 | 66/1000
000000010001 | 57/1000
000000110100 | 64/1000
000000101000 | 69/1000
000000110001 | 56/1000
000000011000 | 68/1000
000000000010 | 63/1000
000000010100 | 64/1000
000000000100 | 58/1000
000000100100 | 68/1000
000000001000 | 66/1000
000000110010 | 67/1000
000000010010 | 69/1000
000000100010 | 47/1000
000000000001 | 63/1000
000000100001 | 55/1000
{'11000000': 253, '10000000': 258, '01000000': 239, '00000000': 250}