Demo 3：CRT 恢复

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

可视化：3 子图（子群轨道+CRT 流程） —— Z_3^* + Z_5^* 拆解后合并

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

实现要点

模型: N=15 = 3 * 5, a=2, 目标 g=8. Shor DLP 在 N 为合数时, 可拆成在 p_i = 3, 5 上的子 DLP (量子周期查找), 然后用 CRT 拼回 x (mod ord_N(a)). 这里: Z_3^: a=2, ord_3(2)=2, 轨道 1 -> 2 -> 1; 8 mod 3 = 2, 故 x_1 = 1. Z_5^: a=2, ord_5(2)=4, 轨道 1 -> 2 -> 4 -> 3 -> 1; 8 mod 5 = 3, 故 x_2 = 3. CRT (n_1=2, n_2=4): x = 3 满足 x mod 2 = 1, x mod 4 = 3, 与期望 x=3 一致.

量子部分用 1-qubit S/T 门做"简化路径" eigenstate QPE, 对应各子群的阶: ord_3(2) = 2 -> S eigenstate phase 1/4 (对应 r=4 in S gate, 或换 T 适配 r=2) ord_5(2) = 4 -> S eigenstate phase 1/4.

完整源码

demos/ch09-algebra-number-theory/03_crt_recovery.py

"""
ch09 演示 3: 用中国剩余定理 (CRT) 拼接各素数子群上的 Shor 离散对数结果.

模型: N=15 = 3 * 5, a=2, 目标 g=8.
Shor DLP 在 N 为合数时, 可拆成在 p_i = 3, 5 上的子 DLP (量子周期查找),
然后用 CRT 拼回 x (mod ord_N(a)). 这里:
  Z_3^*: a=2, ord_3(2)=2, 轨道 1 -> 2 -> 1; 8 mod 3 = 2, 故 x_1 = 1.
  Z_5^*: a=2, ord_5(2)=4, 轨道 1 -> 2 -> 4 -> 3 -> 1; 8 mod 5 = 3, 故 x_2 = 3.
  CRT (n_1=2, n_2=4): x = 3 满足 x mod 2 = 1, x mod 4 = 3, 与期望 x=3 一致.

量子部分用 1-qubit S/T 门做"简化路径" eigenstate QPE, 对应各子群的阶:
  ord_3(2) = 2  -> S eigenstate phase 1/4 (对应 r=4 in S gate, 或换 T 适配 r=2)
  ord_5(2) = 4  -> S eigenstate phase 1/4.

对应教程: chenzhaoyun.com/quantum-tutorial/ch09-algebra-number-theory/finite-rings-ideals.html
依赖: pip install unified-quantum[visualization]
运行: python demos/ch09-algebra-number-theory/03_crt_recovery.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_PREC = 4


def orbit_mult(a, n_mod, start, max_steps=20):
    seen = {start: 0}
    state = start
    cycle = [state]
    while True:
        state = (a * state) % n_mod
        if state in seen:
            break
        seen[state] = len(cycle)
        cycle.append(state)
    return cycle, len(cycle) - seen[state]


def classical_dlog(a, g, n_mod):
    """返回最小 x 使 a^x ≡ g (mod n_mod), 找不到返回 -1."""
    for x in range(n_mod):
        if pow(a, x, n_mod) == g % n_mod:
            return x
    return -1


def crt_pair(x1, n1, x2, n2):
    """对两个未必互素的模 n_1, n_2, 解 x ≡ x_1 (mod n_1), x ≡ x_2 (mod n_2).
    返回 (x, n) 使 x 是该线性同余方程在模 n 下的唯一解.
    若 x_1, x_2 不一致 (mod gcd(n_1, n_2)) 则返回 None."""
    from math import gcd
    g = gcd(n1, n2)
    if (x1 - x2) % g != 0:
        return None
    n = n1 * n2 // g
    # 用扩展欧几里得合并
    def egcd(a, b):
        if b == 0:
            return (a, 1, 0)
        d, u, v = egcd(b, a % b)
        return (d, v, u - (a // b) * v)
    _, u1, _ = egcd(n1, n2 // g)
    x = (x1 * (n2 // g) * u1 + x2 * (n1 // g)) % n
    return x, n


# ---- 1. 经典前置 ----
N = 15
A = 2
G = 8  # 目标: 找 x 使 2^x ≡ 8 (mod 15)
P1, P2 = 3, 5

ord_3 = 2
ord_5 = 4
cycle_p1, r_p1 = orbit_mult(A, P1, 1)
cycle_p2, r_p2 = orbit_mult(A, P2, 1)

g1 = G % P1  # 8 mod 3 = 2
g2 = G % P2  # 8 mod 5 = 3
x1 = classical_dlog(A, g1, P1)  # log_2(2) mod 3 → 1
x2 = classical_dlog(A, g2, P2)  # log_2(3) mod 5 → 3

crt_res = crt_pair(x1, ord_3, x2, ord_5)
x_crt, mod_crt = crt_res  # x = 3 (mod lcm(2, 4)=4)

# 验证
x_true = classical_dlog(A, G, N)  # log_2(8) mod 15 → 3


# ---- 2. QPE on S eigenstate for both subgroups (period proxy) ----
def run_qpe_phase(unitary, prep, n_prec=N_PREC, shots=SHOTS):
    c = qpe_circuit(n_precision=n_prec, unitary_circuit=unitary,
                     state_prep=prep, measure=True)
    tid = submit_task(circuit=c, backend="dummy:local:simulator", shots=shots)
    res = wait_for_result(tid)
    counts = res.counts
    n_states = 2 ** n_prec
    dist = np.zeros(n_states, dtype=float)
    for bs, cnt in counts.items():
        # uniqc qpe bitstring length = n_prec + 1 (末位是未测量 eigenstate),
        # 末 n_prec 字符直接读为二进制整数.
        k = int(bs[-n_prec:], 2)
        dist[k] += cnt
    dist /= dist.sum()
    return dist


unitary_s = Circuit()
unitary_s.s(0)
prep_eigen = Circuit()
prep_eigen.x(0)
phase_dist = run_qpe_phase(unitary_s, prep_eigen, n_prec=N_PREC, shots=SHOTS)
peak_k = int(np.argmax(phase_dist))
peak_phase = peak_k / (2 ** N_PREC)


# ---- 3. matplotlib 可视化: 两个子群轨道 + QPE + CRT 流程 ----
fig, axes = plt.subplots(1, 3, figsize=(16.5, 4.8))

ax = axes[0]
labels_p1 = [f"|{x:02b}>" for x in cycle_p1]
ax.bar(range(len(cycle_p1)), [1.0] * len(cycle_p1), color="#1f77b4", edgecolor="black")
ax.set_xticks(range(len(cycle_p1)))
ax.set_xticklabels(labels_p1, fontsize=10)
for i in range(len(cycle_p1)):
    ax.text(i, 1.02, f"k={i}", ha="center", va="bottom", fontsize=9, color="dimgray")
ax.set_ylim(0, 1.3)
ax.set_xlabel(f"U_a^k |1> in Z_{P1}^* (a={A})")
ax.set_ylabel("probability")
ax.set_title(
    f"Subgroup Z_{P1}^*: order r={r_p1}\n"
    f"g mod {P1} = {g1},  log_{A}({g1}) = x_{P1} = {x1}"
)
ax.grid(axis="y", alpha=0.3)

ax = axes[1]
labels_p2 = [f"|{x:03b}>" for x in cycle_p2]
ax.bar(range(len(cycle_p2)), [1.0] * len(cycle_p2), color="#1f77b4", edgecolor="black")
ax.set_xticks(range(len(cycle_p2)))
ax.set_xticklabels(labels_p2, fontsize=10)
for i in range(len(cycle_p2)):
    ax.text(i, 1.02, f"k={i}\nx={i}", ha="center", va="bottom", fontsize=8, color="dimgray")
ax.set_ylim(0, 1.4)
ax.set_xlabel(f"U_a^k |1> in Z_{P2}^* (a={A})")
ax.set_ylabel("probability")
ax.set_title(
    f"Subgroup Z_{P2}^*: order r={r_p2}\n"
    f"g mod {P2} = {g2},  log_{A}({g2}) = x_{P2} = {x2}"
)
ax.grid(axis="y", alpha=0.3)

ax = axes[2]
xs = np.arange(2 ** N_PREC)
probs = phase_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 (k)  ->  phase = k / 16")
ax.set_ylabel("probability (%)")
ax.set_title(
    f"QPE on S eigenstate |1>  (proxy for r=4)\n"
    f"peak k={peak_k}, phase={peak_phase:.3f} -> r ~ 4 (matches ord_{P2})"
)
for k in (0, 4, 8, 12):
    ax.text(k, probs[k] + 1.5, f"k={k}\nphase={k/(2**N_PREC):.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"CRT recovery of discrete log:  N={N}={P1}*{P2}, a={A}, g={G}  "
    f"(QPE shots={SHOTS})",
    fontsize=12, y=1.02,
)
plt.tight_layout()
plt.savefig(OUT / "03_crt_recovery.png", dpi=120, bbox_inches="tight")
plt.close()


# ---- 4. stdout 关键统计 ----
print("ch09-algebra-number-theory demo 3: CRT recovery of Shor DLP")
print(f"  N = {N} = {P1} * {P2},  base a = {A},  target g = {G}")
print(f"  subgroup Z_{P1}^*: order r1 = {r_p1},  g mod {P1} = {g1},  x_{P1} = {x1}")
print(f"  subgroup Z_{P2}^*: order r2 = {r_p2},  g mod {P2} = {g2},  x_{P2} = {x2}")
print(f"  CRT (x ≡ {x1} mod {r_p1}, x ≡ {x2} mod {r_p2}): x = {x_crt}  (mod {mod_crt})")
print(f"  true x by brute force (mod {r_p1 * r_p2}): {x_true}")
print(f"  agreement (CRT vs classical brute-force): {x_crt == x_true}")
print(f"  verify: {A}^{x_crt} mod {N} = {pow(A, x_crt, N)}  (target g = {G})")
print(f"  QPE phase histogram peak: k = {peak_k}  (phase = {peak_phase:.4f})")
print(f"  QPE peak probability: {phase_dist[peak_k]*100:.2f}%")
print(f"  saved: {OUT / '03_crt_recovery.png'}")

运行结果

运行方式

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

---

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