Python 与量子计算初级教程

1Python 环境搭建

开始量子计算编程之前，我们需要搭建好 Python 环境。本教程使用 NumPy 进行矩阵运算，用纯 Python 模拟量子系统。

安装依赖

bash# 安装核心依赖
pip install numpy matplotlib

# 可选：安装 Qiskit（IBM 量子计算框架）
pip install qiskit qiskit-aer

核心工具

NumPy：用于复数数组和矩阵运算，是量子计算模拟的基础。

Matplotlib：用于绘制量子态、概率分布等可视化图表。

Qiskit（可选）：IBM 开源的量子计算 SDK，可在真实量子计算机上运行电路。

pythonimport numpy as np

# 验证环境
print(f"NumPy 版本: {np.__version__}")
print("环境就绪！开始量子计算之旅 🚀")

2量子比特与 Python

经典比特只能是 0 或 1，而量子比特（qubit）可以同时处于 0 和 1 的叠加态。我们用 Python 中的复数数组来表示量子态。

量子态的数学表示

pythonimport numpy as np

# 基态 |0⟩ 和 |1⟩
ket_0 = np.array([1, 0], dtype=complex)
ket_1 = np.array([0, 1], dtype=complex)

# 叠加态 |+⟩ = (|0⟩ + |1⟩) / √2
ket_plus = np.array([1, 1], dtype=complex) / np.sqrt(2)

print("|0⟩ =", ket_0)
print("|1⟩ =", ket_1)
print("|+⟩ =", np.round(ket_plus, 4))

# 验证归一化
print(f"|+⟩ 归一化: {np.sum(np.abs(ket_plus)**2):.4f}")

3量子门与矩阵运算

量子门是对量子态的操作，数学上是酉矩阵（Unitary Matrix）。用 Python 实现量子门就是做矩阵乘法。

基本量子门

pythonimport numpy as np

# Pauli-X 门（量子 NOT 门）
X = np.array([[0,1],[1,0]], dtype=complex)

# Hadamard 门（创建叠加态）
H = np.array([[1,1],[1,-1]], dtype=complex) / np.sqrt(2)

# Pauli-Y, Pauli-Z 门
Y = np.array([[0,-1j],[1j,0]], dtype=complex)
Z = np.array([[1,0],[0,-1]], dtype=complex)

# 用 Hadamard 门创建叠加态
ket_0 = np.array([1, 0], dtype=complex)
ket_plus = H @ ket_0  # 矩阵乘法
print("H|0⟩ =", np.round(ket_plus, 4))

# X 门翻转 |0⟩ → |1⟩
ket_1 = X @ ket_0
print("X|0⟩ =", ket_1)

4布洛赫球可视化

布洛赫球是表示单量子比特状态的几何工具。球面上的每个点对应一个合法的量子态。用 Python 可以计算并绘制布洛赫球上的状态向量。

布洛赫向量计算

pythondef bloch_vector(alpha, beta):
    """计算量子态 [alpha, beta] 的布洛赫向量"""
    x = 2 * np.real(np.conj(alpha) * beta)
    y = 2 * np.imag(np.conj(alpha) * beta)
    z = np.abs(alpha)**2 - np.abs(beta)**2
    return x, y, z

# |0⟩ → 北极 (0, 0, 1)
print("|0⟩:", bloch_vector(1, 0))
# |+⟩ → 赤道 (1, 0, 0)
print("|+⟩:", bloch_vector(1/np.sqrt(2), 1/np.sqrt(2)))

5量子纠缠编程

量子纠缠是多量子比特系统中最神奇的现象。用 Python 的张量积（tensor product）可以构建多量子比特态。

创建 Bell 态

pythondef tensor(a, b):
    """张量积"""
    return np.kron(a, b)

# Hadamard 门 + CNOT → Bell 态
H = np.array([[1,1],[1,-1]]) / np.sqrt(2)
CNOT = np.array([[1,0,0,0],
                 [0,1,0,0],
                 [0,0,0,1],
                 [0,0,1,0]], dtype=complex)

# |Φ+⟩ = (|00⟩ + |11⟩) / √2
ket_00 = tensor(np.array([1,0]), np.array([1,0]))
state = CNOT @ tensor(H @ np.array([1,0]), np.array([1,0]))
print("Bell态 |Φ+⟩ =", np.round(state, 4))
print("P(00) =", abs(state[0])**2)
print("P(11) =", abs(state[3])**2)

6量子测量与统计

测量使量子态坍缩到某个基态，概率由振幅的模方决定。用 Python 的 random 模块可以模拟量子测量的随机性。

模拟量子测量

pythonimport numpy as np

def measure(state, shots=1000):
    """模拟量子测量 shots 次"""
    probs = np.abs(state)**2
    results = np.random.choice(len(state), size=shots, p=probs)
    counts = {}
    for r in results:
        key = format(r, f'0{int(np.log2(len(state)))}b')
        counts[key] = counts.get(key, 0) + 1
    return counts

# 测量叠加态 |+⟩
ket_plus = np.array([1, 1]) / np.sqrt(2)
counts = measure(ket_plus, shots=1000)
print("测量结果:", counts)

7Deutsch 算法 Python 实现

Deutsch 算法是第一个展示量子计算优势的算法——只需1次查询就能判断函数是常数还是平衡的，而经典计算需要2次。

pythondef deutsch_algorithm(oracle):
    """Deutsch 算法：判断 f 是常数还是平衡函数"""
    # 初始态 |01⟩
    state = np.kron(np.array([1,0]), np.array([0,1]))
    # 对两个量子比特都加 H 门
    HH = np.kron(H, H)
    state = HH @ state
    # 应用 oracle（酉矩阵）
    state = oracle @ state
    # 对第一个量子比特加 H 门
    HI = np.kron(H, np.eye(2))
    state = HI @ state
    # 测量第一个量子比特
    prob_0 = np.abs(state[0])**2 + np.abs(state[1])**2
    return "常数函数" if prob_0 > 0.9 else "平衡函数"

8Grover 搜索 Python 实现

Grover 算法可以在 O(√N) 次查询内从 N 个元素中找到目标，比经典的 O(N) 快得多。我们用 Python 实现 2-qubit Grover（从4个元素中搜索）。

pythondef grover_search(target, n_qubits=2):
    """Grover 搜索算法"""
    N = 2**n_qubits
    # 初始均匀叠加态
    state = np.ones(N, dtype=complex) / np.sqrt(N)
    # Oracle: 翻转目标态的相位
    oracle = np.eye(N, dtype=complex)
    oracle[target, target] = -1
    # 扩散算子
    s = np.ones((N,N), dtype=complex) / N
    diffusion = 2 * s - np.eye(N)
    # 迭代 ≈ π√N/4 次
    iterations = int(np.pi * np.sqrt(N) / 4)
    for _ in range(iterations):
        state = oracle @ state
        state = diffusion @ state
    return np.abs(state)**2