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

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

#コントロール10ビット + ターゲット5ビット
qc = QuantumCircuit(15)

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

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

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

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

qc.barrier()

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

#コントロールビットすべてに適用
for i in range(10):
    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][:10][::-1] #上位ビットから並ぶようにひっくり返しておく
    if pattern in all_pattern:
        all_pattern += result[r]
    else:
        all_pattern = result[r]
print(all_pattern)

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

010101010000010 | 5/1000
110101010011010 | 15/1000
111010101001000 | 1/1000
001010101101000 | 3/1000
（中略）
110010101100001 | 1/1000
001010101000001 | 1/1000
001010101110000 | 1/1000
001111001000010 | 1/1000
{'0101010101': 120, '1101010110': 28, （中略）, '1010010101': 1, '0010000001': 1}

#aを決め、a^(2^i) mod N を計算
a_power = pow(4, 2**i, 21)

000000000010000 | 104/1000
100101010110000 | 1/1000
010111101000100 | 1/1000
110010101010000 | 3/1000
（中略）
010110101000001 | 1/1000
000011010100001 | 1/1000
011100101000100 | 1/1000
000101010100001 | 1/1000
{'0000000000': 347, '1010101001': 8, （中略）, '0101001110': 1, '1010101000': 1}