Demo 2：模幂 U_7

第 9 章 高等代数与数论量子算法 · 本章演示 · 代码演示总索引

可视化：3 子图 —— 16 输入态全测 + QPE

对应教程：第 9 章 高等代数与数论量子算法 中相关小节 · 依赖安装：pip install "unified-quantum[visualization,simulation]"

实现要点

模型: 取 N=15, a=7. U_7 在 Z_15 上的轨道结构 (用于 Shor 周期查找): 固定点: 0, 5, 10, 15 (其中 15 -> 0, 不是真正的"不动"而是溢出折叠). 3 个不相交 4-cycle: (1,7,4,13), (2,14,8,11), (3,6,12,9). 所以周期 r = 4.

量子部分: 因为 uniqc 不支持 add_gate('U', matrix=...), 这里用"单输入目标" 的等价实现: 对每个 x in [0, 15], 构造多受控 X 门, 只在输入为 |x> 时将 diff 位置翻转, 把 |x> 映射到 |7x mod 15>. 控制只用"非 diff"位, 避免互相干扰. 对其它输入态, 电路作用为 identity (按基态正交, 不会污染 |x> 那一支). 16 次独立测量验证全部 a*x mod N 映射.

QPE 部分: 用 1-qubit S 门 eigenstate |1> (与 ch04 demo 1 一致的简化路径), 4-qubit 精度寄存器测相位 1/4 -> r=4.

完整源码

demos/ch09-algebra-number-theory/02_modular_exp.py

"""
ch09 演示 2: Shor 算法的模幂子程序 U_a |x> = |a*x mod N>.

模型: 取 N=15, a=7. U_7 在 Z_15 上的轨道结构 (用于 Shor 周期查找):
  固定点: 0, 5, 10, 15  (其中 15 -> 0, 不是真正的"不动"而是溢出折叠).
  3 个不相交 4-cycle: (1,7,4,13), (2,14,8,11), (3,6,12,9).
  所以周期 r = 4.

量子部分: 因为 uniqc 不支持 add_gate('U', matrix=...), 这里用"单输入目标"
的等价实现: 对每个 x in [0, 15], 构造多受控 X 门, 只在输入为 |x> 时将 diff
位置翻转, 把 |x> 映射到 |7x mod 15>. 控制只用"非 diff"位, 避免互相干扰.
对其它输入态, 电路作用为 identity (按基态正交, 不会污染 |x> 那一支).
16 次独立测量验证全部 a*x mod N 映射.

QPE 部分: 用 1-qubit S 门 eigenstate |1> (与 ch04 demo 1 一致的简化路径),
4-qubit 精度寄存器测相位 1/4 -> r=4.

对应教程: chenzhaoyun.com/quantum-tutorial/ch04-classic-algorithms/shors-algorithm-tutorial.html
依赖: pip install unified-quantum[visualization]
运行: python demos/ch09-algebra-number-theory/02_modular_exp.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_PER_X = 60
N_Q = 4
N_STATES = 2 ** N_Q
N_MOD = 15
A = 7
N_PREC = 4


def u7_map():
    """返回 16 元素数组 perm, perm[x] = 7x mod 15, 并把 15 折叠到 0."""
    perm = [(A * x) % N_MOD for x in range(N_STATES)]
    perm[-1] = 0
    return perm


def orbit(start):
    seen = {start: 0}
    state, cycle = start, [start]
    while True:
        state = (A * state) % N_MOD
        if state in seen:
            return cycle, len(cycle) - seen[state]
        seen[state] = len(cycle)
        cycle.append(state)


def cx_neg(c, target, controls):
    """X(target) 受控于 (qubit, value) 列表 (0/1), 用 pre/post-flip 处理 0."""
    for q, v in controls:
        if v == 0:
            c.x(q)
    c.add_gate('X', qubits=target, control_qubits=[q for q, _ in controls])
    for q, v in controls:
        if v == 0:
            c.x(q)


def build_u7_single_target(x_val):
    """构造对单一输入 |x_val> -> |7x_val mod 15> 的电路. 控制只用非 diff 位."""
    y_val = (A * x_val) % N_MOD if x_val < 15 else 0
    xb = [(x_val >> i) & 1 for i in range(N_Q)]
    yb = [(y_val >> i) & 1 for i in range(N_Q)]
    diff = [i for i in range(N_Q) if xb[i] != yb[i]]
    nd = [i for i in range(N_Q) if xb[i] == yb[i]]
    if not diff:
        return None
    c = Circuit()
    for i in range(N_Q):
        if xb[i]:
            c.x(i)
    for tgt in diff:
        ctrls = [(j, xb[j]) for j in nd if j != tgt]
        if ctrls:
            cx_neg(c, tgt, ctrls)
        else:
            c.x(tgt)
    for q in range(N_Q):
        c.measure(q)
    return c


def decode(bs):
    """uniqc bitstring 直接当二进制读即为整数值."""
    return int(bs[-N_Q:], 2)


# ---- 1. 经典前置 ----
perm = u7_map()
M = np.zeros((N_STATES, N_STATES), dtype=int)
for x, y in enumerate(perm):
    M[y, x] = 1
cycle1, period_r = orbit(1)
cycle2, _ = orbit(2)
cycle3, _ = orbit(3)

# 16 次测量验证 7x mod 15 映射
measured_y, expected_y = [], []
for x in range(N_STATES):
    y_exp = perm[x]
    expected_y.append(y_exp)
    if x == y_exp:
        measured_y.append(y_exp)
        continue
    c = build_u7_single_target(x)
    res = wait_for_result(submit_task(circuit=c, backend="dummy:local:simulator",
                                      shots=SHOTS_PER_X))
    measured_y.append(decode(max(res.counts, key=res.counts.get)))

# ---- 2. QPE on S eigenstate (proxy for r=4) ----
unitary = Circuit(); unitary.s(0)
prep = Circuit(); prep.x(0)
c_qpe = qpe_circuit(n_precision=N_PREC, unitary_circuit=unitary,
                     state_prep=prep, measure=True)
qpe_counts = wait_for_result(submit_task(circuit=c_qpe, backend="dummy:local:simulator",
                                         shots=2000)).counts
qpe_dist = np.zeros(N_STATES, dtype=float)
for bs, cnt in qpe_counts.items():
    qpe_dist[int(bs[-N_PREC:], 2)] += cnt
qpe_dist /= qpe_dist.sum()


# ---- 3. 可视化: 排列热图 + 验证柱图 + QPE 相位直方图 ----
fig, axes = plt.subplots(1, 3, figsize=(16.5, 4.8))
labels_n4 = [f"{i:04b}" for i in range(N_STATES)]

ax = axes[0]
ax.imshow(M, cmap="Blues", vmin=0, vmax=1)
for i in range(N_STATES):
    for j in range(N_STATES):
        if M[i, j] == 1:
            ax.text(j, i, "1", ha="center", va="center", fontsize=7, color="red")
ax.set_title(f"U_7 map  |x> -> |{A}x mod {N_MOD}>  (16x16, 15->0 fold, r={period_r})")
ax.set_xlabel("input  |x>")
ax.set_ylabel("output  |7x mod 15>")
ax.set_xticks(np.arange(N_STATES))
ax.set_yticks(np.arange(N_STATES))
ax.set_xticklabels(labels_n4, rotation=0, fontsize=7)
ax.set_yticklabels(labels_n4, fontsize=7)

ax = axes[1]
xs = np.arange(N_STATES); w = 0.4
agree = sum(1 for i in range(N_STATES) if expected_y[i] == measured_y[i])
ax.bar(xs - w / 2, expected_y, width=w, color="#1f77b4", edgecolor="black",
       label=f"expected  {A}x mod {N_MOD}")
ax.bar(xs + w / 2, measured_y, width=w, color="#d62728", edgecolor="black",
       label=f"measured (shots={SHOTS_PER_X})")
ax.set_xticks(xs)
ax.set_xticklabels(labels_n4, rotation=0, fontsize=7)
ax.set_xlabel("input  |x>"); ax.set_ylabel("output index")
ax.set_title(f"Verification: x -> {A}x mod {N_MOD} for all 16 inputs\nagreement = {agree}/{N_STATES}")
ax.set_ylim(0, N_STATES); ax.legend(fontsize=8, loc="upper left"); ax.grid(axis="y", alpha=0.3)

ax = axes[2]
probs = qpe_dist * 100
ax.bar(np.arange(N_STATES), probs,
       color=["#d62728" if k in (0, 4, 8, 12) else "#1f77b4" for k in range(N_STATES)],
       edgecolor="black")
ax.set_xticks(np.arange(N_STATES))
ax.set_xticklabels([f"{k:04b}" for k in range(N_STATES)], rotation=0, fontsize=7)
ax.set_xlabel("counting register (k)  ->  phase = k / 16"); ax.set_ylabel("probability (%)")
ax.set_title("QPE on eigenstate |1> of S gate (period-4 proxy)\npeak at k=4 (phase 1/4) confirms r = 4")
for k in (0, 4, 8, 12):
    ax.text(k, probs[k] + 1.5, f"k={k}\nphase={k/N_STATES:.3f}",
            ha="center", va="bottom", fontsize=7, color="darkred")
ax.set_ylim(0, max(probs.max() * 1.25, 30))
ax.grid(axis="y", alpha=0.3)
fig.suptitle(f"Modular exponentiation U_{A} on Z_{N_MOD}: "
             f"permutation matrix + measured mapping + QPE phase histogram",
             fontsize=12, y=1.02)
plt.tight_layout()
plt.savefig(OUT / "02_modular_exp.png", dpi=120, bbox_inches="tight")
plt.close()


# ---- 4. stdout 关键统计 ----
peak_k = int(np.argmax(qpe_dist))
print("ch09-algebra-number-theory demo 2: modular exponentiation U_7 on N=15")
print(f"  per-x shots: {SHOTS_PER_X}  (16 runs total = {16 * SHOTS_PER_X} shots)")
print(f"  orbit structure:")
print(f"    cycle 1: {[f'{x:04b}' for x in cycle1]}  (r = {period_r})")
print(f"    cycle 2: {[f'{x:04b}' for x in cycle2]}")
print(f"    cycle 3: {[f'{x:04b}' for x in cycle3]}")
print(f"  agreement (16 inputs): {agree}/{N_STATES}")
disc = [(i, measured_y[i], expected_y[i])
        for i in range(N_STATES) if measured_y[i] != expected_y[i]]
print(f"  discrepancies: {disc}")
print(f"  QPE peak: k = {peak_k}  (phase = {peak_k/N_STATES:.4f})")
print(f"  QPE peak probability: {qpe_dist[peak_k]*100:.2f}%")
print(f"  saved: {OUT / '02_modular_exp.png'}")

运行结果

运行方式

# 仓库根目录执行；输出 PNG 到 demos/ch09-algebra-number-theory/02_modular_exp.png，关键统计到 stdout
.venv/bin/python3 demos/ch09-algebra-number-theory/02_modular_exp.py

---

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