量子傅里叶变换与相位估计深入

QFT|j⟩ = (1/√N) Σ_{k=0}^{N-1} exp(2πi·j·k/N) |k⟩，N = 2^n

等价的分量形式：
QFT|j⟩ = (1/√N) ⊗_{l=1}^{n} [ |0⟩ + exp(2πi·j/2^l)|1⟩ ]

→ 关键洞察：QFT 把「整数 j」变成「二进制小数 0.j_1 j_2 ... j_n」
  在相位上逐位展开

对每个比特 l：
  1. 施加 H 门到比特 l
  2. 对每个更低位 k > l，施加受控相位门 R_k：
     R_k = diag(1, exp(2πi/2^(k-l+1)))
→ 最后（可选）交换比特顺序（bit-reversal）

门数：
  n 个 H 门
  n(n-1)/2 个受控相位门
→ 总门数 O(n^2)，对比经典 FFT 的 O(N log N) = O(n·2^n)

经典 FFT：给定 N 个复数的完整向量，输出 N 个复数
  代价 O(N log N) = O(n·2^n)
  → 输入输出都是「N 个经典数」

量子 QFT：作用在量子态 |j⟩ 上，输出叠加态
  代价 O(n^2) 门
  → 但输出是「量子态」，无法直接读出全部 N 个振幅

QFT 本身并没有「比 FFT 快指数倍地算出傅里叶变换」
  经典：O(N log N) 算出全部 N 个系数
  量子：O(n^2) 得到叠加态，但读出任意一个振幅都需 O(N) 次测量

→ QFT 的价值不在于「单独使用」
  而在于「作为更大算法的一个子程序」
  只要后续操作能「在傅里叶基上」完成而不需读出
  指数加速才成立

U|u⟩ = exp(2πi·φ)|u⟩，φ ∈ [0, 1) 未知
目标：估计 φ 的 n 位二进制近似 0.φ_1 φ_2 ... φ_n

寄存器：
  R1：n 个「计数比特」，初态 |0⟩^⊗n
  R2：特征态寄存器，初态 |u⟩

步骤：
  1. H^⊗n 作用于 R1 → 均匀叠加 (1/√N) Σ_j |j⟩
  2. 受控 U^(2^j)：R1 的第 j 位控制 U 的 2^j 次幂
     在 |u⟩ 上累积相位 exp(2πi·φ·2^j)
  3. R1 的态变为 (1/√N) Σ_j exp(2πi·φ·j)|j⟩
  4. 对 R1 做逆 QFT（QFT†）
  5. 测量 R1 → 得到 φ 的 n 位二进制近似

步骤 3 的态恰是 QFT|φ_approx⟩ 的形式
  （相位 2πiφj 对应 QFT 定义中的 2πijk/N，k = φ·N）
→ 逆 QFT 把它「反变换」回计算基，集中到 |φ_approx⟩
→ 测量即读出 φ 的二进制近似

朴素做法：对第 j 位，连续施加 2^j 次受控 U → 指数级门数
高效做法：反复「平方」U
  U^2 = U·U，U^4 = U^2·U^2，...
  只需 j 次平方即可得到 U^(2^j)
→ 总代价 O(n) 次受控 U（每次是一份 U 的线路）
→ 这是 QPE 高效的关键：平方代替重复

设 φ = φ* + δ，φ* 为 n 位近似，0 ≤ δ < 2^(-n)
测量得到 φ* 的概率：
  P(φ*) = (1/N^2) · |Σ_{j=0}^{N-1} exp(2πiδj)|^2
        = (1/N^2) · sin^2(πδN) / sin^2(πδ)
  最小概率出现在 δ = 1/2（恰在两点中间）：
  P ≥ 4/π^2 ≈ 0.405

1. 增加计数比特数 n：
   精度 ε = 2^(-n)，每加一个比特精度翻倍
   但代价是线路深度与比特数线性增长

2. 额外增加 t 个比特（n → n + t）：
   测量后「四舍五入」到 n 位
   成功概率 ≥ 1 - 1/(2·2^t) = 1 - 2^(-t-1)
   → t = 3 时成功概率 ≥ 0.9375
   → t = 10 时成功概率 ≥ 0.9995

核心思想：每次只保留 1 个计数比特，逐位从低位到高位估计

流程（第 k 轮，从最低位开始）：
  1. 制备单个计数比特为 |+⟩
  2. 受控 U^(2^(n-k)) 施加到特征态
  3. 相位反馈：用已知的高位相位做一个 R_z 修正
  4. H 门 + 测量 → 得到第 k 位
→ 比特数从 n 降到 1（加特征态寄存器）
   代价：需要「经典反馈」与逐轮串行，线路深度仍 O(n^2) 但比特省

思路：不显式做逆 QFT，而是「逐位测量相位」
  对 k = 1, 2, 4, 8, ...（2 的幂）分别估计 φ 的小数部分
  利用 Hadamard 测试读出 cos(2π·2^k·φ)
  再用经典后处理重建 φ 的二进制展开
特点：对「单次测量」友好，适合早期硬件
缺点：噪声鲁棒性差，需要重复采样

Shor 因数分解：
  把「找周期 r」转化为「相位估计」问题
  U|y⟩ = |a·y mod N⟩，其特征相位 φ ≈ s/r
  用连分数从 φ 恢复 r → 分解 N

HHL 解线性方程组：
  对 A 做相位估计，把特征值 λ_j 提取到寄存器
  再用受控旋转做 1/λ_j 缩放

量子化学（能量估计）：
  对 e^(-iHt) 做相位估计，提取基态能量 E_0
  → 精度 ε 需计数比特 n ≈ log(1/ε)

量子模拟（振幅估计）：
  用 QPE 或振幅放大估计振幅

总比特数 = n（计数寄存器）+ n_sys（特征态/系统寄存器）+ 辅助比特

计数比特 n 由精度决定：n = log2(1/ε) + t（t 为额外比特）
系统比特 n_sys 由问题规模决定：
  Shor 分解 N：n_sys ≈ 2·log2(N)
  量子化学：n_sys = 轨道数 × 2（自旋轨道）
→ 总比特数是「精度 + 问题规模」的叠加

主要开销：受控 U^(2^j) 需要 U 的线路重复 O(n) 次（平方序列）
若 U 的线路深度为 D_U，则 QPE 深度 ≈ n · D_U + O(n^2)（逆 QFT）
→ 逆 QFT 的 O(n^2) 通常不是瓶颈，瓶颈是 n · D_U

例（Shor 分解 2048 位 RSA）：
  需要约 4000 ~ 6000 个逻辑比特
  需要约 10^9 个 T 门（表面码下折算到物理比特 10^6 量级）
→ 这就是「量子优势门槛」的具体数字来源

1. 深度：n · D_U + O(n^2) 远超当前相干时间
   当前设备 T2 ~ 100 μs，门时间 ~ 100 ns → 相干窗口 ~ 10^3 门
   而 QPE 需要 10^6 门以上

2. 精度：噪声使相位信息随机化，估计值分布变宽
   要压制噪声需重复测量 O(1/ε^2) 次

3. 比特数：容错所需的物理比特数远超当前设备

理想 QPE：测量分布集中在 φ_approx
含噪 QPE：分布被展宽，出现「卫星峰」
  - 退相干把主峰压低、尾部抬高
  - 读出误差使峰位偏移
→ 用重复测量 + 中位数/众数估计可部分恢复
  但代价随噪声指数增长

from qiskit import QuantumCircuit
import numpy as np

def qft(n, inverse=False):
    """构造 n 比特 QFT（inverse=True 时构造逆 QFT）"""
    qc = QuantumCircuit(n, name="IQFT" if inverse else "QFT")
    for j in range(n):
        qc.h(j)
        for k in range(j + 1, n):
            angle = np.pi / (2 ** (k - j))
            if inverse:
                angle = -angle
            qc.cp(angle, k, j)
    if not inverse:                 # 正 QFT 末尾做比特反序
        for j in range(n // 2):
            qc.swap(j, n - 1 - j)
    return qc

from qiskit import QuantumCircuit, transpile
from qiskit_aer import AerSimulator
import numpy as np

def qpe_circuit(phi, n_count=4):
    """用 n_count 个计数比特估计单比特相位 phi"""
    qc = QuantumCircuit(n_count + 1, n_count)
    qc.x(n_count)                          # 特征态 |1⟩
    for q in range(n_count):
        qc.h(q)                            # 计数寄存器均匀叠加
    # 受控 U^(2^j)：U = p(2*pi*phi)，U^(2^j) = p(2*pi*phi*2^j)
    for q in range(n_count):
        qc.cp(2 * np.pi * phi * (2 ** q), q, n_count)
    qc.append(qft(n_count, inverse=True), range(n_count))
    qc.measure(range(n_count), range(n_count))
    return qc

qc = qpe_circuit(0.25, n_count=4)
sim = AerSimulator()
result = sim.run(transpile(qc, sim), shots=4096).result()
counts = result.get_counts()
print(counts)      # 预期集中在 |0100⟩（0.25 = 二进制 0.0100）