大规模量子模拟新突破：cuQuantum SDK v25.11技术详解

import cupy as cp
from cuquantum.bindings import cupauliprop
from cuquantum import cudaDataType

# 创建库句柄和工作空间描述符
handle = cupauliprop.create()
workspace = cupauliprop.create_workspace_descriptor(handle)

# 为工作空间分配GPU内存
ws_size = 1024 * 1024 * 64  # 示例：64 MiB
d_ws = cp.cuda.alloc(ws_size)
cupauliprop.workspace_set_memory(
    handle, workspace, cupauliprop.Memspace.DEVICE,
    cupauliprop.WorkspaceKind.WORKSPACE_SCRATCH, d_ws.ptr, ws_size)

# 辅助函数：将Pauli字符串编码为打包整数（每个量子位2位：X和Z掩码）
def encode_pauli(num_qubits, paulis, qubits):
    num_ints = cupauliprop.get_num_packed_integers(num_qubits)
    # 打包整数格式: [X_ints..., Z_ints...]
    packed = np.zeros(num_ints * 2, dtype=np.uint64)
    x_mask, z_mask = packed[:num_ints], packed[num_ints:]
    for p, q in zip(paulis, qubits):
        idx, bit = divmod(q, 64)
        if p in (cupauliprop.PauliKind.PAULI_X, cupauliprop.PauliKind.PAULI_Y):
            x_mask[idx] |= (1 << bit)
        if p in (cupauliprop.PauliKind.PAULI_Z, cupauliprop.PauliKind.PAULI_Y):
            z_mask[idx] |= (1 << bit)
    return packed

# 1. 分配设备缓冲区
# 定义容量（最大Pauli字符串数）并分配缓冲区
max_terms = 10000
num_packed_ints = cupauliprop.get_num_packed_integers(num_qubits)
d_pauli = cp.zeros((max_terms, 2 * num_packed_ints), dtype=cp.uint64, order="C")
d_coef = cp.zeros(max_terms, dtype=cp.float64, order="C")

# 2. 填充初始可观测量 (Z_62)
encoded_pauli = encode_pauli(num_qubits, [cupauliprop.PauliKind.PAULI_Z], [62])
# 分配第一项
d_pauli[0] = cp.array(encoded_pauli)
d_coef[0] = 1.0

# 3. 创建Pauli展开式
# 输入展开式：用可观测量预填充
expansion_in = cupauliprop.create_pauli_expansion(
    handle, num_qubits,
    d_pauli.data.ptr, d_pauli.nbytes,
    d_coef.data.ptr, d_coef.nbytes,
    cudaDataType.CUDA_R_64F,
    1, 1, 1  # num_terms=1, is_sorted=True, is_unique=True)

# 输出展开式：初始为空 (num_terms=0)，需要自己的缓冲区
d_pauli_out = cp.zeros_like(d_pauli)
d_coef_out = cp.zeros_like(d_coef)

expansion_out = cupauliprop.create_pauli_expansion(
    handle, num_qubits,
    d_pauli_out.data.ptr, d_pauli_out.nbytes,
    d_coef_out.data.ptr, d_coef_out.nbytes,
    cudaDataType.CUDA_R_64F,
    0, 0, 0)

# 在量子位0上创建Z旋转门
paulis = [cupauliprop.PauliKind.PAULI_Z]
qubits = [0]
gate = cupauliprop.create_pauli_rotation_gate_operator(
    handle, theta, 1, qubits, paulis)

# 获取输入展开式中当前项的视图
num_terms = cupauliprop.pauli_expansion_get_num_terms(handle, expansion_out)
view = cupauliprop.pauli_expansion_get_contiguous_range(
    handle, expansion_in, 0, num_terms)

# 应用门：in_expansion -> gate -> out_expansion
cupauliprop.pauli_expansion_view_compute_operator_application(
    handle, view, expansion_out, gate,
    True,         # adjoint?
    False, False, # make_sorted?, keep_duplicates?
    0, None,      # 截断策略（可选）
    workspace)

import numpy as np
result = np.zeros(1, dtype=np.float64)

# 计算迹
cupauliprop.pauli_expansion_view_compute_trace_with_zero_state(
    handle, view, result.ctypes.data, workspace)

import cuquantum.stabilizer as cust
# 电路信息
num_qubits = 5
num_shots = 10_000
num_measurements = 2

# 在GPU上创建电路
circ = cust.Circuit("""
H 0 1
X_ERROR(0.1) 1 2
DEPOLARIZE2(0.5) 2 3
CX 0 1 2 3
M 0 3
""")
sim = cust.FrameSimulator(
      num_qubits,
      num_shots,
      num_measurements)
sim.apply(circ)

circ2 = cust.Circuit("""
Z_ERROR(0.01) 1 4
""")
sim.apply(circ2)

# 获取测量翻转
m_table = sim.get_measurement_bits(bit_packed=False)
print(m_table.dtype)
# uint8
print(m_table.shape)
# (2, 10000)
print(m_table)
# [[0 0 0 ... 0 0 0]
#  [1 0 0 ... 0 1 1]]

x_table, z_table = sim.get_pauli_xz_bits(bit_packed=True)
print(x_table.dtype)
# uint8
print(x_table.shape)
# (5, 1252)

# 获取pauli表
pauli_table = sim.get_pauli_table()
num_frames_print = 5
for i in range(num_frames_print):
    print(pauli_table[i])
# ...XZ
# ZXX..
# ...Z.
# .....
# ...Z.

import stim

p = 0.001
circ_stim = stim.Circuit.generated(
    "surface_code:rotated_memory_z",
    distance=5,
    rounds=5,
    after_clifford_depolarization=p,
    after_reset_flip_probability=p,
    before_measure_flip_probability=p,
    before_round_data_depolarization=p,)
circ = cust.Circuit(circ_stim)
sim = cust.FrameSimulator(
    circ_stim.num_qubits,
    num_shots,
    circ_stim.num_measurements,
    num_detectors=circ_stim.num_detectors,)
sim.apply(circ)

pauli_table = sim.get_pauli_table()
for i in range(num_frames_print):
    print(pauli_table[i])