Demo 1：Shor 离散对数

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

可视化：双图（轨道+直方图） —— N=15 离散对数骨架 + QPE 相位

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

实现要点

模型: 给定 N=15 (模合数), a=2, 寻找 x 使 g = a^x mod N. 乘法群 Z_15^* 的阶 phi(15)=8; a=2 在其中的阶为 4, 因为 2^1 = 2, 2^2 = 4, 2^3 = 8, 2^4 = 16 mod 15 = 1. 因此 x 仅在模 4 意义下唯一. Shor 算法的核心是把离散对数 归约到乘法算子 U_a |y> = |a*y mod N> 的周期查找: |QPE 估相位 r / 2^n> -> r = 4 -> x 在轨道 [1, 2, 4, 8] 中的位置.

量子部分用 1 比特 S 门做"简化路径" eigenstate (类似 ch04 demo 1): |1> 是 S 的本征态, 特征值 i = e^{2pii/4}, 用 4-qubit 精度寄存器 跑 QPE 应集中在 phase = 1/4 (即 k=4 of 16). 该相位对应 r = 4. 然后对每个目标 g in {1, 2, 4, 8}, 由 r 与轨道位置反推 x.

完整源码

demos/ch09-algebra-number-theory/01_discrete_log.py

"""
ch09 演示 1: Shor 离散对数小实例 N=15, a=2.

模型: 给定 N=15 (模合数), a=2, 寻找 x 使 g = a^x mod N.
乘法群 Z_15^* 的阶 phi(15)=8; a=2 在其中的阶为 4, 因为
2^1 = 2, 2^2 = 4, 2^3 = 8, 2^4 = 16 mod 15 = 1.
因此 x 仅在模 4 意义下唯一. Shor 算法的核心是把离散对数
归约到乘法算子 U_a |y> = |a*y mod N> 的周期查找:
|QPE 估相位 r / 2^n> -> r = 4 -> x 在轨道 [1, 2, 4, 8] 中的位置.

量子部分用 1 比特 S 门做"简化路径" eigenstate (类似 ch04 demo 1):
|1> 是 S 的本征态, 特征值 i = e^{2*pi*i/4}, 用 4-qubit 精度寄存器
跑 QPE 应集中在 phase = 1/4 (即 k=4 of 16). 该相位对应 r = 4.
然后对每个目标 g in {1, 2, 4, 8}, 由 r 与轨道位置反推 x.

对应教程: chenzhaoyun.com/quantum-tutorial/ch09-algebra-number-theory/finite-rings-ideals.html
依赖: pip install unified-quantum[visualization]
运行: python demos/ch09-algebra-number-theory/01_discrete_log.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_MOD = 15
A = 2
N_PREC = 4
N_STATES = 2 ** N_PREC  # 16


def orbit_mult(a, n_mod, start, max_steps=12):
    """返回从 |start> 出发, 应用 U_a (乘 a mod n_mod) 的轨道, 及周期 r.
    周期定义为轨道中首次出现重复时, 两次出现位置之差."""
    seen = {}
    state = start
    cycle = []
    period = -1
    for k in range(max_steps):
        if state in seen:
            period = k - seen[state]
            break
        seen[state] = k
        cycle.append(state)
        state = (a * state) % n_mod
    return cycle, period


def classical_dlog(a, g, n_mod, max_x=20):
    """经典暴力法: 找 x 使 a^x mod n_mod == g. 返回最小 x."""
    for x in range(max_x):
        if pow(a, x, n_mod) == g:
            return x
    return -1


def dlog_from_orbit(a, g, orbit):
    """由轨道与目标 g 反推 x (假定 g 在轨道中)."""
    rev = {val: idx for idx, val in enumerate(orbit)}
    return rev.get(g, -1)


# ---- 1. 经典前置: 算轨道、r、离散对数 ----
full_cycle, period_r = orbit_mult(A, N_MOD, 1)
targets = [1, 2, 4, 8]  # 轨道上的元素
known_x = [classical_dlog(A, g, N_MOD) for g in targets]

# ---- 2. uniqc: 1-qubit S gate eigenstate QPE, 模拟周期 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,
)
tid = submit_task(circuit=c_qpe, backend="dummy:local:simulator", shots=SHOTS)
qpe_result = wait_for_result(tid)
qpe_counts = qpe_result.counts


def precision_from_bitstring(bs, n_prec):
    """从 QPE 总 bitstring 中取精度寄存器, 转成整数.
    qpe_circuit 的 bitstring 长度为 n_prec + 1 (含 1 个未测量的 eigenstate qubit,
    显示在 bitstring 末尾固定为 '0'). 精度寄存器占据前 n_prec 个字符,
    按 q[1] 为 MSB 的二进制读法 (无需翻转)."""
    return int(bs[-n_prec:], 2)


qpe_dist = np.zeros(N_STATES, dtype=float)
for bs, cnt in qpe_counts.items():
    k = precision_from_bitstring(bs, N_PREC)
    qpe_dist[k] += cnt
qpe_dist /= qpe_dist.sum()


# ---- 3. 由 QPE 直方图推断 r, 然后由 r 反推 x ----
# QPE 峰位应为 4 (= N_STATES / period_r), period 4.
# 选 top-1 peak 推断周期
peak_k = int(np.argmax(qpe_dist))
estimated_period = N_STATES // peak_k if peak_k > 0 else N_STATES
# 简化: 用真实周期 4 做反推 (演示算法意图)
estimated_x = [dlog_from_orbit(A, g, full_cycle) for g in targets]

# ---- 4. matplotlib 可视化: 轨道柱图 + QPE 相位直方图 ----
fig, axes = plt.subplots(1, 2, figsize=(13.0, 4.6))

ax = axes[0]
orbit_labels = [f"|{x:04b}>" for x in full_cycle]
colors = ["#1f77b4"] * len(full_cycle)
ax.bar(range(len(full_cycle)), [1.0] * len(full_cycle), color=colors, edgecolor="black")
ax.set_xticks(range(len(full_cycle)))
ax.set_xticklabels(orbit_labels, fontsize=9)
for i in range(len(full_cycle)):
    ax.text(i, 1.02, f"k={i}\nx={i}", ha="center", va="bottom", fontsize=8, color="dimgray")
ax.set_ylim(0, 1.3)
ax.set_xlabel("U_a^k applied to |1>")
ax.set_ylabel("probability")
ax.set_title(
    f"Multiplication-by-{A} orbit in Z_{N_MOD}^*\n"
    f"|1> -> |{full_cycle[1]:04b}> -> |{full_cycle[2]:04b}> -> |{full_cycle[3]:04b}> -> |1>   (period r = {period_r})"
)
ax.grid(axis="y", alpha=0.3)

ax = axes[1]
xs = np.arange(N_STATES)
probs = qpe_dist * 100
bar_colors = ["#d62728" if (k in (0, 4, 8, 12)) else "#1f77b4" for k in xs]
ax.bar(xs, probs, color=bar_colors, edgecolor="black")
ax.set_xticks(xs)
ax.set_xticklabels([f"{k:04b}" for k in xs], rotation=0, fontsize=7)
ax.set_xlabel("counting register (4 bits) -> phase = k / 16")
ax.set_ylabel("probability (%)")
ax.set_title(
    f"QPE on eigenstate |1> of S gate (proxy for U_a)\n"
    f"shots={SHOTS}, peak at k=4 (phase=1/4) => 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=8, color="darkred")
ax.set_ylim(0, max(probs.max() * 1.25, 25))
ax.grid(axis="y", alpha=0.3)

fig.suptitle(
    f"Shor discrete-log skeleton: N={N_MOD}, a={A}, base order r={period_r}  "
    f"(QPE shots={SHOTS}, dummy:local:simulator)",
    fontsize=12, y=1.02,
)
plt.tight_layout()
plt.savefig(OUT / "01_discrete_log.png", dpi=120, bbox_inches="tight")
plt.close()


# ---- 5. stdout 关键统计 ----
peak_probs = qpe_dist
total_top = peak_probs[peak_k]
known_str = ", ".join(f"g={g} -> x={known_x[i]}" for i, g in enumerate(targets))
est_str = ", ".join(f"g={g} -> x={estimated_x[i]}" for i, g in enumerate(targets))

print("ch09-algebra-number-theory demo 1: Shor DLP N=15, a=2")
print(f"  shots: {SHOTS}")
print(f"  multiplicative orbit of |1> under *{A} mod {N_MOD}: {[f'{x:04b}' for x in full_cycle]}")
print(f"  order (period) r = {period_r}")
print(f"  QPE peak counting register: k = {peak_k}  (phase = {peak_k/N_STATES:.4f})")
print(f"  QPE peak probability: {total_top*100:.2f}%")
print(f"  phase 1/4 corresponds to r = {N_STATES // 4}")
print(f"  discrete log (classical brute-force): {known_str}")
print(f"  discrete log (estimated from orbit, r={period_r}): {est_str}")
all_match = all(known_x[i] == estimated_x[i] for i in range(len(targets)))
print(f"  agreement: {all_match}")
print(f"  saved: {OUT / '01_discrete_log.png'}")

运行结果

运行方式

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

---

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