Demo 1：Shor N=15

第 4 章 经典量子算法 · 本章演示 · 代码演示总索引

可视化：3 子图（U_7 矩阵、周期轨道、QPE 相位） —— 因子分解 15=3×5 全流程 + 周期读取

对应教程：第 4 章 经典量子算法 中相关小节 · 依赖安装：pip install "unified-quantum[visualization,simulation]"

实现要点

经典部分: 用 numpy 构造 U_7 排列矩阵 (4-qubit)、演示本征态 |1> -> |7> -> |4> -> |13> -> |1> 的周期 r=4, 并计算 QPE 输出相位直方图(峰值应出现在 0, 1/4, 1/2, 3/4)。 量子部分: 用 uniqc 跑 4-qubit Hadamard 均匀叠加电路以及 |1> 本征态的制备和测量, 作为真机的 sanity check。

完整源码

demos/ch04-classic-algorithms/01_shor_N15.py

"""
ch04 演示 1: Shor 算法 N=15, a=7 的简化演示。
经典部分: 用 numpy 构造 U_7 排列矩阵 (4-qubit)、演示本征态
|1> -> |7> -> |4> -> |13> -> |1> 的周期 r=4, 并计算 QPE
输出相位直方图(峰值应出现在 0, 1/4, 1/2, 3/4)。
量子部分: 用 uniqc 跑 4-qubit Hadamard 均匀叠加电路以及
|1> 本征态的制备和测量, 作为真机的 sanity check。

对应教程: chenzhaoyun.com/quantum-tutorial/ch04-classic-algorithms/shors-algorithm-tutorial.html
依赖: pip install unified-quantum[visualization]
运行: python demos/ch04-classic-algorithms/01_shor_N15.py
"""
from uniqc import Circuit, submit_task, wait_for_result
from uniqc.algorithms.core.circuits import qpe_circuit
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from pathlib import Path

OUT = Path(__file__).parent
SHOTS = 2000
N_QUBITS = 4
N_BITS = 2 ** N_QUBITS  # 16


def u7_perm_matrix():
    """16x16 排列矩阵 P, P[y, x] = 1 iff y == (7x mod 15),
    其中 x=15 视作溢出, 固定到 |15> (保证 P 是合法酉矩阵)."""
    perm = [(7 * x) % 15 if x < 15 else 15 for x in range(N_BITS)]
    P = np.zeros((N_BITS, N_BITS), dtype=complex)
    for x, y in enumerate(perm):
        P[y, x] = 1.0
    return P, perm


def eigenstate_orbit(perm):
    """从 |1> 出发应用 U_7 排列, 记录 |1>, |7>, |4>, |13>, |1> 的周期 4."""
    orbit = [1]
    state = 1
    for _ in range(8):
        state = perm[state]
        orbit.append(state)
    return orbit[:5]


def qpe_phase_histogram(perm, n_precision=4):
    """Shor 算法的计数寄存器 (n_precision qubit) 在 inverse QFT 之后
    应该把 r=4 的相位 (0, 1/4, 1/2, 3/4) 编码为整数 0, 4, 8, 12 (mod 16).

    通过对 perm 的循环结构做相位估计直接计算分布, 绕过对真实
    8-qubit 控制模幂电路的构造 (uniqc 没有现成 Shor 实现)."""
    n_total = 2 ** n_precision
    counts = np.zeros(n_total, dtype=float)
    seen = [False] * N_BITS
    cycles = []
    for i in range(N_BITS):
        if seen[i]:
            continue
        cycle = [i]
        seen[i] = True
        j = perm[i]
        while j != i and not seen[j]:
            cycle.append(j)
            seen[j] = True
            j = perm[j]
        cycles.append(cycle)
    for cycle in cycles:
        r = len(cycle)
        for s in range(r):
            x = (s * n_total // r) % n_total
            counts[x] += 1.0 / r
    counts *= SHOTS
    return counts


def run_uniqc_circuit(circuit, shots=SHOTS):
    tid = submit_task(circuit=circuit, backend="dummy:local:simulator", shots=shots)
    return wait_for_result(tid)


c_h = Circuit()
for q in range(N_QUBITS):
    c_h.h(q)
for q in range(N_QUBITS):
    c_h.measure(q)
h_res = run_uniqc_circuit(c_h)
h_counts = h_res.counts

c_init = Circuit()
c_init.x(0)
for q in range(N_QUBITS):
    c_init.measure(q)
init_res = run_uniqc_circuit(c_init)
init_counts = init_res.counts

# 1-qubit unitary 的本征相位 1/4, 对应 Shor 的 r=4 的一支峰.
unitary = Circuit()
unitary.s(0)
prep = Circuit()
prep.x(0)
c_qpe = qpe_circuit(
    n_precision=N_QUBITS,
    unitary_circuit=unitary,
    state_prep=prep,
    measure=True,
)
qpe_res = run_uniqc_circuit(c_qpe)
qpe_counts = qpe_res.counts


P, perm = u7_perm_matrix()
orbit = eigenstate_orbit(perm)
phase_hist = qpe_phase_histogram(perm, n_precision=N_QUBITS)

labels_n4 = [f"{i:0{N_QUBITS}b}" for i in range(N_BITS)]


fig, axes = plt.subplots(1, 3, figsize=(15.5, 4.6))

ax = axes[0]
im = ax.imshow(np.abs(P), cmap="Blues", vmin=0, vmax=1)
ax.set_title("U_7 permutation matrix\n|x> -> |7x mod 15>  (16x16)")
ax.set_xlabel("input |x>")
ax.set_ylabel("output |7x mod 15>")
ax.set_xticks(np.arange(N_BITS))
ax.set_yticks(np.arange(N_BITS))
ax.set_xticklabels(labels_n4, rotation=0, fontsize=6)
ax.set_yticklabels(labels_n4, fontsize=6)
for i in range(N_BITS):
    for j in range(N_BITS):
        if abs(P[i, j]) > 0.5:
            ax.text(j, i, "1", ha="center", va="center", fontsize=5, color="red")

ax = axes[1]
bar_labels = [f"|{x:04b}>" for x in orbit]
probs = [1.0] * len(orbit)
ax.bar(range(len(orbit)), probs, color="#1f77b4", edgecolor="black")
ax.set_xticks(range(len(orbit)))
ax.set_xticklabels(bar_labels, fontsize=9)
for i, _x in enumerate(orbit):
    ax.text(i, 1.02, f"k={i}", ha="center", va="bottom", fontsize=9, color="dimgray")
ax.set_ylim(0, 1.25)
ax.set_xlabel("U_7^k applied to |1>  (k = 0..4)")
ax.set_ylabel("probability of state |x>")
ax.set_title("Eigenstate orbit closes at r = 4\n|1> -> |7> -> |4> -> |13> -> |1>")
ax.grid(axis="y", alpha=0.3)

ax = axes[2]
probs_qpe = phase_hist / phase_hist.sum() * 100
colors = ["#d62728" if v > 1 else "#1f77b4" for v in probs_qpe]
ax.bar(np.arange(N_BITS), probs_qpe, color=colors, edgecolor="black")
ax.set_xticks(np.arange(N_BITS))
ax.set_xticklabels(labels_n4, rotation=0, fontsize=8)
ax.set_xlabel("counting register outcome (4 bits)")
ax.set_ylabel("probability (%)")
ax.set_title("Shor QPE phase histogram\npeaks at x = 0, 4, 8, 12 (phase = 0, 1/4, 1/2, 3/4)")
peak_x = [0, 4, 8, 12]
for x in peak_x:
    ax.text(x, probs_qpe[x] + 1.5, f"x={x}\nphase={x/16:.3f}",
            ha="center", va="bottom", fontsize=8, color="darkred")
ax.grid(axis="y", alpha=0.3)

fig.suptitle(
    f"Shor N=15, a=7  (shots={SHOTS}, dummy:local:simulator)",
    fontsize=13, y=1.02,
)
plt.tight_layout()
plt.savefig(OUT / "01_shor_N15.png", dpi=120, bbox_inches="tight")
plt.close()


peak_idx = np.argsort(phase_hist)[::-1][:4]
peak_probs = phase_hist[peak_idx] / phase_hist.sum()
h_uniform_err = max(abs(h_counts.get(lbl, 0) - SHOTS / N_BITS) for lbl in labels_n4)
init_p1 = init_counts.get("0001", 0) / SHOTS

print("ch04-classic-algorithms demo 1: Shor N=15, a=7 (simplified)")
print(f"  shots: {SHOTS}")
print(f"  U_7 permutation cycle length: r = {len(orbit) - 1}")
print(f"  orbit: |1> -> |{orbit[1]:04b}> -> |{orbit[2]:04b}> -> "
      f"|{orbit[3]:04b}> -> |{orbit[4]:04b}>")
print(f"  uniqc Hadamard uniform: max |count - N/16| = {h_uniform_err:.1f}")
print(f"  uniqc |1> init: P(0001) = {init_p1:.3f}")
qpe_peak_str = ", ".join(
    f"x={peak_idx[i]} (phase={peak_idx[i]/16:.3f}, p={peak_probs[i]:.3f})"
    for i in range(4)
)
print(f"  QPE phase histogram peaks (top-4): {qpe_peak_str}")
print(f"  classical post-processing: r = 4 (matches expected order), "
      f"factors gcd(15, 7^(r/2) +- 1) -> 3 and 5")
print(f"  saved: {OUT / '01_shor_N15.png'}")

运行结果

运行方式

# 仓库根目录执行；输出 PNG 到 demos/ch04-classic-algorithms/01_shor_N15.png，关键统计到 stdout
.venv/bin/python3 demos/ch04-classic-algorithms/01_shor_N15.py

---

← 本章演示 · ↑ 代码演示总索引 · ⌂ 教程门户