Demo 1：Shor 破 RSA

第 16 章 译码与量子密码分析 · 本章演示 · 代码演示总索引

可视化：单子图 —— 3+3 qubit QPE 估 r=4 + GCD

对应教程：第 16 章 译码与量子密码分析 中相关小节 · 依赖安装：pip install "unified-quantum[visualization,simulation]"

实现要点

核心：Shor 算法把整数分解归约为周期寻找。给定 N = 15，随机选 a = 7， f(x) = 7^x mod 15 以 r = 4 为周期。量子部分用相位估计 (QPE) 提取 s/r 的足够精确二进近似；经典后处理通过连分数展开恢复 r，再用 gcd(a^(r/2) ± 1, N) 得到非平凡因子。本 demo 用 uniqc 的 qpe_circuit helper 构造周期 4 的酉算子 U（在 |000⟩-|011⟩ 子空间上的本征值为 e^{2πi s/4}, s=0,1,2,3），跑 3-比特 QPE，测出 phase s/4，然后 连分数得到 r = 4，最后 gcd(7^2 ± 1, 15) = 3, 5 恢复因子。

完整源码

demos/ch16-cryptanalysis-coding/01_shor_rsa.py

"""
ch16-cryptanalysis-coding demo 1: Shor's algorithm factoring N = 15 = 3 x 5

对应教程：chenzhaoyun.com/quantum-tutorial/ch16-cryptanalysis-coding/quantum-cryptanalysis.html
核心：Shor 算法把整数分解归约为周期寻找。给定 N = 15，随机选 a = 7，
     f(x) = 7^x mod 15 以 r = 4 为周期。量子部分用相位估计 (QPE) 提取
     s/r 的足够精确二进近似；经典后处理通过连分数展开恢复 r，再用
     gcd(a^(r/2) ± 1, N) 得到非平凡因子。本 demo 用 uniqc 的 qpe_circuit
     helper 构造周期 4 的酉算子 U（在 |000⟩-|011⟩ 子空间上的本征值为
     e^{2πi s/4}, s=0,1,2,3），跑 3-比特 QPE，测出 phase s/4，然后
     连分数得到 r = 4，最后 gcd(7^2 ± 1, 15) = 3, 5 恢复因子。

依赖：pip install unified-quantum[visualization]
运行：python demos/ch16-cryptanalysis-coding/01_shor_rsa.py
"""
from __future__ import annotations

from fractions import Fraction
from math import gcd
from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

from uniqc import Circuit, submit_task, wait_for_result
from uniqc.algorithms.core.circuits.qpe import qpe_circuit

OUT = Path(__file__).parent
N_TARGET = 15
A_BASE = 7
N_PRECISION = 3  # 8 buckets; enough to resolve r = 4
N_SYSTEM = 3     # 3-qubit system with |000>,|001>,|010>,|011> as the r=4 eigenspace
SHOTS = 1500


# === 1. Classical pre-check: verify period of f(x) = 7^x mod 15 ===
def classical_period(a: int, n: int, max_x: int = 32) -> int:
    for r in range(1, max_x + 1):
        if pow(a, r, n) == 1:
            return r
    return -1


def classical_powers(a: int, n: int, count: int) -> list[int]:
    return [pow(a, x, n) for x in range(count)]


powers = classical_powers(A_BASE, N_TARGET, 9)
true_r = classical_period(A_BASE, N_TARGET)
print("ch16-cryptanalysis-coding-demo-1: Shor for N = 15 = 3 x 5")
print(f"  classical pre-check: f(x) = {A_BASE}^x mod {N_TARGET}")
print(f"    first 9 powers: {powers}")
print(f"    detected period r = {true_r}")


# === 2. Build QPE circuit for the modular-exponentiation unitary U ===
# A textbook Shor with N = 15, a = 7 needs U_7 acting on 4 qubits. We model it
# with a 3-qubit diagonal unitary U whose |000>,|001>,|010>,|011> subspace
# has eigenvalues exactly e^{2πi s/4} for s = 0, 1, 2, 3. This gives the same
# period-4 phase structure as the real 7^x mod 15 modular exponentiation on
# its 4-cycle eigenspace.
#   RZ(π/2) on q0  ->  contributes α to |q0=1> states
#   RZ(π) on q1    ->  contributes β to |q1=1> states
#   RZ(-3π/2) on q2 ->  sets the global phase so the 4 target phases match {0, π/2, π, 3π/2}
unitary = Circuit()
unitary.rz(0, np.pi / 2.0)
unitary.rz(1, np.pi)
unitary.rz(2, -3.0 * np.pi / 2.0)

state_prep = Circuit()
state_prep.h(0); state_prep.h(1); state_prep.h(2)  # |+++>


# === 3. Submit QPE and collect phase-register samples ===
qpe = qpe_circuit(
    n_precision=N_PRECISION,
    unitary_circuit=unitary,
    state_prep=state_prep,
    measure=True,
)

task_id = submit_task(circuit=qpe, backend="dummy:local:simulator", shots=SHOTS)
result = wait_for_result(task_id)
counts = result.counts
total = sum(int(v) for v in counts.values())


# === 4. Convert phase-register samples to integer s in [0, 2^N) ===
def bs_to_int(bs: str, n: int) -> int:
    return int(bs, 2) if len(bs) == n else int(bs[-n:], 2)


n_buckets = 1 << N_PRECISION
s_hist = np.zeros(n_buckets, dtype=np.int64)
for bs, cnt in counts.items():
    s_hist[bs_to_int(bs, N_PRECISION)] += int(cnt)
s_probs = s_hist / max(total, 1)


# === 5. Continued-fraction post-processing: each s -> fraction s / 2^N -> r ===
def fraction_from_phase(s: int, n: int) -> Fraction:
    return Fraction(s, 1 << n).limit_denominator(64)


# Expected theoretical phases s/r for r = 4 with N_PRECISION bits
expected_s = [0, 2, 4, 6]   # s/8 = 0, 1/4, 1/2, 3/4


# === 6. Classical post-processing: gcd(a^(r/2) ± 1, N) ===
def factor_from_period(a: int, n: int, r: int) -> tuple[int, int]:
    if r <= 0 or r % 2 != 0:
        return -1, -1
    half = pow(a, r // 2, n)
    if half == n - 1:
        return -1, -1
    g1 = gcd(half - 1, n)
    g2 = gcd(half + 1, n)
    if 1 < g1 < n:
        return g1, n // g1
    if 1 < g2 < n:
        return g2, n // g2
    return -1, -1


# Use the most-observed NON-ZERO s as the canonical measurement
# (s = 0 is uninformative: it would correspond to phase 0, which makes r ambiguous)
nonzero_s = [s for s in range(n_buckets) if s != 0 and s_hist[s] > 0]
top_s = max(nonzero_s, key=lambda s: s_hist[s])
top_phase = Fraction(top_s, n_buckets)
top_r = top_phase.denominator if top_phase.numerator != 0 else 0
factor_1, factor_2 = factor_from_period(A_BASE, N_TARGET, top_r)


# === 7. matplotlib visualisation: phase histogram + post-processing chain ===
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13.0, 5.2))

# Left: phase-register histogram
buckets = np.arange(n_buckets)
bar_colors = ["#d62728" if int(s) in expected_s else "#1f77b4" for s in buckets]
ax1.bar(buckets, s_probs, color=bar_colors, ec="black", lw=0.5)
for s in expected_s:
    ax1.axvline(s, color="#d62728", ls="--", lw=0.9, alpha=0.4)
ax1.text(0.5, s_probs.max() * 0.92,
         f"expected s for r=4: {expected_s}",
         ha="center", fontsize=10, color="#7a1f1f", style="italic",
         bbox=dict(facecolor="white", alpha=0.85, ec="none"))
ax1.set_xlabel(f"phase-register integer s (over {n_buckets} buckets, "
               fr"phase = $s/2^{N_PRECISION}$)")
ax1.set_ylabel("measured probability")
ax1.set_xticks(buckets)
ax1.set_xticklabels([format(s, f"0{N_PRECISION}b") for s in buckets], fontsize=9)
ax1.set_ylim(0, s_probs.max() * 1.15 + 0.005)
ax1.set_title(f"QPE on U (period-4 modular exponentiation)\n"
              f"{N_PRECISION} precision qubits, {total} shots")
ax1.grid(alpha=0.3, axis="y")

# Right: factor-extraction flow diagram (text panels)
ax2.axis("off")
half_pow = pow(A_BASE, top_r // 2, N_TARGET)
flow_lines = [
    (0.96, f"N = {N_TARGET}, pick a = {A_BASE} (gcd(a,N) = {gcd(A_BASE, N_TARGET)})"),
    (0.86, f"compute 7^x mod 15 -> period r = {true_r}"),
    (0.75, f"build QPE with period-4 unitary U (3 system qubits)"),
    (0.64, f"measure phase register: most-likely s = {top_s:0{N_PRECISION}b}"),
    (0.53, f"continued fraction: {top_s}/{n_buckets} = {top_phase} -> r = {top_r}"),
    (0.41, f"check r even?  r = {top_r} -> "
           f"{'YES' if top_r % 2 == 0 else 'NO'}"),
    (0.29, f"compute 7^(r/2) mod 15 = 7^{top_r // 2} mod 15 = {half_pow}"),
    (0.17, f"gcd(7^(r/2) - 1, 15) = {gcd(half_pow - 1, N_TARGET)}"),
    (0.05, f"gcd(7^(r/2) + 1, 15) = {gcd(half_pow + 1, N_TARGET)}"),
    (-0.06, f"=>  {N_TARGET} = {factor_1} x {factor_2}  "
            f"({'SUCCESS' if factor_1 * factor_2 == N_TARGET else 'FAIL'})"),
]
for y, txt in flow_lines:
    ax2.text(0.02, y, txt, ha="left", va="top", fontsize=10, family="monospace")
ax2.set_xlim(0, 1)
ax2.set_ylim(-0.18, 1.05)
ax2.set_title("Shor's algorithm: classical post-processing chain\n"
              "(from phase measurement to non-trivial factor)")

plt.suptitle(f"ch16-cryptanalysis-coding demo 1: Shor for N = {N_TARGET} = "
             f"{factor_1} x {factor_2}", y=1.02)
plt.tight_layout()
out_path = OUT / "01_shor_rsa.png"
plt.savefig(out_path, dpi=120, bbox_inches="tight")
plt.close()


# === 8. stdout key statistics ===
print(f"\n  N_PRECISION = {N_PRECISION} ({n_buckets} buckets), shots: {total}")
print(f"  distribution at expected s values: "
      f"{{0: {s_hist[0]/total:.3f}, 2: {s_hist[2]/total:.3f}, "
      f"4: {s_hist[4]/total:.3f}, 6: {s_hist[6]/total:.3f}}}")
print(f"  continued fraction {top_s}/{n_buckets} -> {top_phase}, recovered r = {top_r}")
print(f"  a^(r/2) mod N = {half_pow}, gcd(a^(r/2)-1,N)={gcd(half_pow-1, N_TARGET)}, "
      f"gcd(a^(r/2)+1,N)={gcd(half_pow+1, N_TARGET)}")
print(f"  FACTORED: {N_TARGET} = {factor_1} x {factor_2}")
print(f"  PNG -> {out_path}")

运行结果

运行方式

# 仓库根目录执行；输出 PNG 到 demos/ch16-cryptanalysis-coding/01_shor_rsa.png，关键统计到 stdout
.venv/bin/python3 demos/ch16-cryptanalysis-coding/01_shor_rsa.py

---

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