Demo 2：Abel HSP 离散对数

第 10 章 隐藏结构与群算法 · 本章演示 · 代码演示总索引

可视化：双图 —— N=15 拆 p1, p2 + CRT

对应教程：第 10 章 隐藏结构与群算法 中相关小节 · 依赖安装：pip install "unified-quantum[visualization,simulation]"

实现要点

问题：求 x 使 2^x ≡ 8 (mod 15)，期望 x=3。 Shor 把阶查找拆成对 N 的每个素因子 p 分别查找，再 CRT (lcm) 拼合：

• p=3: ord_3(2) = 2 - p=5: ord_5(2) = 4 - lcm(2, 4) = 4 = r ✓ QPE 估相位：U|a⟩ = |2a mod p⟩ 的本征相位是 j/r_p。 实操：用 numpy 模拟 QPE 输出（直接构造 U 所需门代价过高），CRT 与恢复 x 经典完成。

完整源码

demos/ch10-hidden-structures/02_abel_hsp.py

"""
ch10 演示 2：Abel HSP / 离散对数 + CRT period finding（N=15, a=2）

对应教程：chenzhaoyun.com/quantum-tutorial/ch10-hidden-structures/abelian-hidden-subgroup.html
依赖：pip install unified-quantum[visualization]
运行：python demos/ch10-hidden-structures/02_abel_hsp.py

问题：求 x 使 2^x ≡ 8 (mod 15)，期望 x=3。
Shor 把阶查找拆成对 N 的每个素因子 p 分别查找，再 CRT (lcm) 拼合：
  - p=3: ord_3(2) = 2    - p=5: ord_5(2) = 4    - lcm(2, 4) = 4 = r  ✓
QPE 估相位：U|a⟩ = |2a mod p⟩ 的本征相位是 j/r_p。
实操：用 numpy 模拟 QPE 输出（直接构造 U 所需门代价过高），CRT 与恢复 x 经典完成。
"""
from uniqc import Circuit, submit_task, wait_for_result
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from pathlib import Path
from math import gcd
from fractions import Fraction

OUT = Path(__file__).parent
N, A, H, SHOTS = 15, 2, 8, 1500


def simulate_qpe_phase(period: int, n_precision: int, shots: int, rng) -> np.ndarray:
    """模拟 QPE 估相位：j/period → k ≈ j·Q/period + 高斯噪声。"""
    Q = 1 << n_precision
    js = rng.integers(0, period, size=shots)
    sigma = max(0.5, np.sqrt(Q) / 4.0)
    ks = np.round(js * Q / period + rng.normal(0, sigma, size=shots)).astype(int) % Q
    return ks


def recover_period_via_cf(k: int, Q: int, max_period: int) -> int:
    """连分数：从 k/Q 恢复候选 period。"""
    if k == 0:
        return Q
    frac = Fraction(k, Q)
    for d in range(1, max_period + 1):
        if abs(frac - Fraction(round(frac * d), d)) < 1 / (2 * Q):
            return d
    return Q


def order_mod(a, p):
    for r in range(1, p + 1):
        if pow(a, r, p) == 1:
            return r
    return -1


def factor_via_order(a, N, r):
    """已知 a 的阶 r → gcd(a^(r/2) ± 1, N) 给出非平凡因子。"""
    if r % 2 != 0:
        return None
    for delta in (-1, 1):
        g = gcd(pow(a, r // 2, N) + delta, N)
        if 1 < g < N:
            return g, N // g
    return None


# === 1. 经典检查：列出 2^x mod 15 与部分周期 ===
r_full = order_mod(A, N)
r1 = order_mod(A, 3)
r2 = order_mod(A, 5)
print(f"target: {A}^x ≡ {H} (mod {N}); true x = 3, full period r = {r_full}")
print(f"CRT split: r_1 = ord_3({A}) = {r1}, r_2 = ord_5({A}) = {r2}")
print(f"classical powers: {[pow(A, x, N) for x in range(r_full + 1)]}")
print()

# === 2. 模拟 QPE 输出 ===
rng = np.random.default_rng(seed=2025)
samples_p1 = simulate_qpe_phase(r1, n_precision=3, shots=SHOTS, rng=rng)
samples_p2 = simulate_qpe_phase(r2, n_precision=3, shots=SHOTS, rng=rng)

# === 3. 连分数恢复 r_1, r_2，lcm 拼合 ===
def top_two_peaks(xs):
    vals, counts = np.unique(xs, return_counts=True)
    order = np.argsort(-counts)
    out = [int(vals[order[0]])]
    if len(order) > 1:
        out.append(int(vals[order[1]]))
    return out


k1 = next((k for k in top_two_peaks(samples_p1) if k != 0), 0)
k2 = next((k for k in top_two_peaks(samples_p2) if k != 0), 0)
recovered_r1 = recover_period_via_cf(k1, 8, max_period=N)
recovered_r2 = recover_period_via_cf(k2, 8, max_period=N)
r_recovered = recovered_r1 * recovered_r2 // gcd(recovered_r1, recovered_r2)

# === 4. 因子分解与离散对数恢复 ===
factors = factor_via_order(A, N, r_recovered)
x_recovered = next((x for x in range(r_recovered) if pow(A, x, N) == H), None)

# === 5. matplotlib 可视化 ===
fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# 左图：QPE 输出相位直方图（p=3 与 p=5 并列）
ax = axes[0]
bins = np.arange(0, 9)
h1, _ = np.histogram(samples_p1, bins=bins)
h2, _ = np.histogram(samples_p2, bins=bins)
w = 0.4
ax.bar(np.arange(8) - w / 2, h1 / SHOTS * 100, width=w, color="#1f77b4",
       edgecolor="black", label=f"p=3 (r={r1}): k = j·Q/r")
ax.bar(np.arange(8) + w / 2, h2 / SHOTS * 100, width=w, color="#d62728",
       edgecolor="black", label=f"p=5 (r={r2}): k = j·Q/r")
for r in (r1, r2):
    for j in range(r):
        k = round(j * 8 / r)
        ax.axvline(k, color="gray", linestyle=":", alpha=0.4)
        ax.annotate(f"j={j}", xy=(k, 50), xytext=(k, 52),
                    ha="center", fontsize=7, color="gray")
ax.set_xticks(np.arange(8))
ax.set_xticklabels([f"k={k}" for k in range(8)])
ax.set_xlabel("Measured precision register k (Q = 2^3 = 8)")
ax.set_ylabel("Probability (%)")
ax.set_title(
    f"Simulated QPE phase estimation for p=3 and p=5\n"
    f"({SHOTS} shots each, peaks at k = j·Q/r_p)",
    fontsize=10,
)
ax.set_ylim(0, 60)
ax.grid(axis="y", alpha=0.3)
ax.legend(loc="upper right", fontsize=9)

# 右图：CRT 拼合流程 + 离散对数恢复（纯文本）
ax = axes[1]
ax.axis("off")
sep = "─" * 48
txt = f"CRT combination of partial periods\n{sep}\n\n"
txt += f"  p_1 = 3:  ord_3({A}) = {r1}\n"
txt += f"  QPE sample k_1 = {k1}/Q → r_1 = {recovered_r1}  {'✓' if recovered_r1 == r1 else '✗'}\n\n"
txt += f"  p_2 = 5:  ord_5({A}) = {r2}\n"
txt += f"  QPE sample k_2 = {k2}/Q → r_2 = {recovered_r2}  {'✓' if recovered_r2 == r2 else '✗'}\n\n"
txt += f"  full period: lcm(r_1, r_2) = lcm({r1}, {r2}) = {r_recovered}\n\n"
txt += f"Factoring via order\n{sep}\n\n"
if factors is not None:
    f1, f2 = factors
    txt += f"  r = {r_recovered} (even) → gcd({A}^(r/2) ± 1, {N})\n"
    txt += f"  factors of {N} = ({f1}, {f2})  ✓\n\n"
else:
    txt += f"  r = {r_recovered} (odd or trivial) → try other branch\n\n"
txt += f"Discrete log recovery\n{sep}\n\n"
txt += f"  find x with {A}^x ≡ {H} (mod {N}); search x ∈ [0, r={r_recovered}):\n"
for x_try in range(r_recovered):
    val = pow(A, x_try, N)
    marker = "  ← MATCH" if val == H else ""
    txt += f"    {A}^{x_try} mod {N} = {val}{marker}\n"
txt += f"\n  recovered x = {x_recovered}  (true x = 3, match = {x_recovered == 3})"
ax.text(0.02, 0.98, txt, family="monospace", fontsize=9,
        verticalalignment="top", transform=ax.transAxes)

fig.suptitle(
    f"Abelian HSP discrete log on Z_{N}: a={A}, h={H}  "
    f"({SHOTS} QPE shots per prime, dummy:local:simulator)",
    fontsize=11, y=1.02,
)
plt.tight_layout()
plt.savefig(OUT / "02_abel_hsp.png", dpi=120, bbox_inches="tight")
plt.close()

# === 6. stdout 关键统计 ===
print("ch10-hidden-structures demo 2: Abel HSP discrete log + CRT period finding")
print(f"  QPE sample k_1 = {k1} → r_1 recovered = {recovered_r1} (true {r1})")
print(f"  QPE sample k_2 = {k2} → r_2 recovered = {recovered_r2} (true {r2})")
print(f"  combined r = lcm(r_1, r_2) = {r_recovered}")
print(f"  factorization of {N} = {factors}")
print(f"  discrete log: recovered x = {x_recovered} (true x = 3, match = {x_recovered == 3})")
print(f"  saved: {OUT / '02_abel_hsp.png'}")

运行结果

运行方式

# 仓库根目录执行；输出 PNG 到 demos/ch10-hidden-structures/02_abel_hsp.png，关键统计到 stdout
.venv/bin/python3 demos/ch10-hidden-structures/02_abel_hsp.py

---

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