企业官网建设流程全解析

简介：面向量子计算入门者与MATLAB实践者的算法实现资料，围绕量子位叠加与纠缠、常见量子逻辑门及量子电路仿真展开，覆盖量子傅里叶变换、无序搜索、大整数因数分解等经典算法，适合课程实验、自学演练和小规模量子系统模拟。压缩包共5个文件，以4个.m脚本和1份.md说明为主，整体体量约5KB，脚本便于直接运行与改写，说明文档可帮助梳理项目入口、实验先后顺序和算法脉络。已有1883人学习，读者可据此搭建量子门与测量流程，复现搜索和因数分解示例，并理解目录结构与代码组织方式。对希望把量子算法从公式落到MATLAB代码、进一步连接Qiskit或真实硬件验证的读者，这份轻量资料提供了可操作的起点与排错参照。

1. 态向量加酉矩阵就够了：MATLAB 里跑量子计算算法的真实手感

不少人第一次听说用 MATLAB 做量子计算，第一反应是「那不是 Python 的地盘吗」。真动手会发现 MATLAB 在这件事上有它的便利：n 比特寄存器的态就是一个 2^n 维复向量，量子门就是酉矩阵，kron拼张量积、sparse压稀疏、eig求特征值、plot画概率分布，全是现成的工具箱函数，不用装任何量子 SDK。这份资源就是把量子计算里几套经典算法——Deutsch-Jozsa、Grover、量子傅里叶变换、相位估计——从态向量和酉矩阵这一层老老实实实现了一遍。它适合两类人：一类是要在课程设计或论文里复现算法、需要出概率分布图的人；另一类是想弄明白「门到底是怎么乘到态上去的」，不愿意停在调用现成库的人。

2. 工程骨架：core/algorithms/viz 三层怎么分，比特序怎么定

拿到一份算法代码，我先看的不是算法本身，是目录。因为算法写得再漂亮，只要比特序和下标约定不统一，跑出来的结果就是一团随机数，而且极难查。这一章先把骨架和约定定死，后面所有代码都在这套约定上跑。

2.1 目录分层：哪些文件是引擎，哪些只是脚本

一个能长期改的 MATLAB 量子算法工程，常见骨架大致是下面这样。分成core、algorithms、viz、tests四层，核心思想是让「造门」和「用门」彻底分开——core里的函数不认识任何具体算法，algorithms里的函数不自己拼矩阵。

quantum-matlab/ ├── core/ % 引擎层：只关心态和门，不含任何算法逻辑 │ ├── basisState.m % 构造计算基态 |x> │ ├── hadamard.m % 常用门的快捷构造 │ ├── singleQubitGate.m % 单比特门 -> 2^n 维算符 │ ├── controlledGate.m % 受控门 -> 2^n 维算符 │ ├── runCircuit.m % 按顺序把算符作用到态上 │ └── measureAll.m % 采样与态塌缩 ├── algorithms/ % 算法层：只调用 core，不碰 kron │ ├── deutschJozsa.m │ ├── grover.m │ ├── qftMatrix.m │ └── phaseEstimation.m ├── viz/ % 出图层 │ ├── plotProb.m │ └── plotBloch.m ├── tests/ % 断言层，改完必跑 │ └── testCore.m └── demo.m % 入口脚本，一次跑完所有例子

这么分层的好处很直接：要加一个新算法，只在algorithms里加一个文件；要换一个稀疏实现，只改core里那两三个造门函数，算法文件一行不动。血泪经验是很多人图省事，把kron直接写在算法脚本里，等到比特数一涨、想换稀疏路径时，得挨个脚本改。

2.2 为什么用 MATLAB：kron、sparse、eig 与画图链路

选 MATLAB 不是因为它跑得快——纯矩阵乘法 Python 的 NumPy 一点不慢。真正的差别在三个地方。第一是kron的语义和量子力学的张量积完全对得上，kron(A, B)就是把 A 作用在高位比特、B 作用在低位比特，写法几乎和公式一一对应。第二是sparse对受控门这类「单位阵加几个非零元」的算符极其友好，一个 4096×4096 的稠密复矩阵要 512 MB，稀疏存下来通常只有几万个非零元。第三是从算到画是一条链：算完概率幅直接bar出分布图，算完单比特态直接画 Bloch 球，中间不用导数据。

代价也得说清楚。MATLAB 没有原生的复数张量广播，多比特批量演化（比如一次跑 100 组参数）写得比 NumPy 啰嗦；for循环在旧版本里慢，得靠向量化。这套资源里的写法是「小比特数走稠密、大比特数走稀疏 + 分块」，下面 3.3 节会给出两套等价实现。

2.3 比特序与下标：big-endian 约定和 MATLAB 的 1 起下标

这是整个项目最容易翻车的地方，必须在写第一行代码前定死。约定如下：

只要这三条守住，后面所有门矩阵和态向量就是自洽的。反过来，只要有一处按 little-endian 写（第 n 位当最高位），结果不会报错，只会「看起来像随机」。我一般会在tests里放一条断言把这件事钉死：

% 断言：第 1 位是最高位。|10> 的整数下标应为 2，落在 psi(3) psi = basisState(2, 2); assert(abs(psi(3) - 1) < 1e-12 && nnz(psi) == 1, '比特序约定被破坏');

basisState(2, 2)的第二个参数是 0 起的整数下标，2 的二进制是10，最高位为 1，对应第 1 位为 1、第 2 位为 0。如果哪天有人把kron的顺序写反，这条断言会第一时间炸出来，比跑到 Grover 里去猜「为什么成功率只有 50%」省太多时间。

3. 核心模块拆解：一个门从矩阵乘到态向量上的完整链路

这一章把core里那几个文件逐个拆开。它们加起来不到一百行，但整个项目所有算法都建立在这几个函数上，写错一个，后面全错。

3.1 态准备：从比特串到 2^n 维列向量

量子态在 MATLAB 里就是一个列向量，长度 2^n，元素是复数概率幅。构造基态的函数只有三行，但下标换算必须和 2.3 节的约定严格对齐：

function psi = basisState(n, idx) % BASISSTATE 构造 n 比特计算基态 |idx>，idx 从 0 开始 % n : 比特数 % idx : 基态对应的整数，范围 0 .. 2^n-1，第 1 位为最高位 dim = 2^n; psi = zeros(dim, 1); psi(idx + 1) = 1; % 唯一的 +1 转换点，别在别处再转 end

逻辑说明：整个项目只在「整数下标 ↔ MATLAB 下标」这一处做加减一，其他函数一律接收 0 起的整数下标，内部再转。参数n决定维度，idx决定是哪个基矢。这么做的好处是算法层写basisState(n, 1)时心里想的就是|0...01>，不用记 MATLAB 的偏移。

初态搞定后，叠加态用一行就够：对全零态每位打一个 H，或者直接ones(dim,1)/sqrt(dim)。后者只有在 H 作用在|0...0>上时才等价，别拿它当通用公式用。

3.2 单比特门与受控门：kron 拼矩阵和稀疏受控写法

单比特门的构造靠kron把 2×2 矩阵扩到 2^n 维。核心思路是：在第target位的位置放这个 2×2 矩阵，其余位置放单位阵，从高位到低位依次kron过去。

function U = singleQubitGate(gate, n, target) % SINGLEQUBITGATE 把 2x2 的单比特门作用到 n 比特寄存器的第 target 位 % gate : 2x2 酉矩阵，例如 H、X、Rz(theta) % n : 总比特数 % target : 目标位编号，1 表示最高位 I2 = eye(2); U = 1; % kron 的单位元，第一轮 kron(1, x) 就是 x for k = 1:n if k == target U = kron(U, gate); else U = kron(U, I2); end end end function H = hadamard() % HADAMARD 返回 2x2 的 Hadamard 门 H = [1 1; 1 -1] / sqrt(2); end

逻辑说明：kron是左结合地从高位往低位拼，循环变量k从 1 走到 n，正好对应第 1 位到第 n 位。参数target一旦写错位置，kron不会报错，只会把门作用到别的比特上，所以这个函数必须配单测——singleQubitGate(X, 2, 1)作用在|00>上应该得到|10>，也就是psi(3)=1。

受控门不能这么拼，因为kron表达不了「控制位为 1 时才作用」这种条件。常见做法是直接按基矢逐行填：遍历所有 2^n 个整数下标，看控制位是不是 1，是就在目标位上把gate的四个元素填进去，不是就填一个 1。

function U = controlledGate(gate, n, ctrl, target) % CONTROLLEDGATE 把 2x2 的 gate 变成 n 比特上的受控门 % ctrl : 控制位编号（1 起），target : 目标位编号（1 起），两者不能相同 dim = 2^n; U = zeros(dim); % 先全零，避免单位阵残留 for s = 0:dim-1 if bitget(s, n - ctrl + 1) == 1 % 控制位为 1：在目标位上作用 gate i0 = bitset(s, n - target + 1, 0); % 目标位清零 step = 2^(n - target); % 目标位的权重 i1 = i0 + step; U(i0+1, i0+1) = gate(1,1); U(i0+1, i1+1) = gate(1,2); U(i1+1, i0+1) = gate(2,1); U(i1+1, i1+1) = gate(2,2); else U(s+1, s+1) = 1; % 控制位为 0，什么都不做 end end end

逻辑说明：bitget(s, n-ctrl+1)把「第 ctrl 位（1 起，最高位在前）」换算成 MATLABbitget从最低位数起的位数，这是全项目第二处需要小心换算的地方。bitset(s, ..., 0)把目标位清零，加2^(n-target)就得到目标位为 1 的另一个下标，这两个下标构成一对，正好把gate的四个元素写进去。因为每个s只会落进「控制位为 1」或「控制位为 0」其中一支，且每支覆盖全部行一次，所以不会出现覆盖冲突。

CNOT 就是controlledGate([0 1; 1 0], n, ctrl, target)。Toffoli 需要两个控制位，可以在这个函数基础上再加一层判断，或者直接用kron拼 Toffoli 的 8×8 矩阵再扩维，后者更快。

提示：controlledGate返回的是稠密阵，n=10 时就是 1024×1024，构造本身不慢，但如果你要在一百次循环里反复构造同一个门，把它挪到循环外面，否则耗时会莫名其妙翻几倍。

3.3 线路求值：逐门作用、批量矩阵与不构造大矩阵的分块法

有了门矩阵，跑线路就是一路左乘：

function psi = runCircuit(psi, ops) % RUNCIRCUIT 依次把 ops 里的算符作用到态向量上 % psi : 2^n 维列向量；ops : cell 数组，每个元素是 2^n x 2^n 的算符 for k = 1:numel(ops) psi = ops{k} * psi; end end

逻辑说明：ops{k} * psi是矩阵向量乘，MATLAB 会自动挑 BLAS 的zgemv，小规模下比自己写循环快一到两个数量级。但这条路的天花板很明显：每个算符都是 2^n × 2^n，n=12 时单个算符 4096×4096 复数双精度就是 512 MB，ops里存上二十个直接爆内存。这就是很多人跑到 10 比特以上就卡住的原因。

绕开的办法是根本不构造大矩阵，直接对态向量做分块运算。单比特门作用在第target位，等价于把态向量按 2^(n−target+1) 的长度切块，每块的前一半是目标位为|0>的分量、后一半是|1>的分量，然后对每块做一次 2×2 混合：

function psi = applySingleFast(psi, gate, n, target) % APPLYSINGLEFAST 不构造 2^n 维矩阵，直接把单比特门作用到态上 % 内存开销只有态向量本身，n=14 也能跑 blk = 2^(n - target + 1); % 每块包含目标位的 |0> 和 |1> 分量 M = reshape(psi, blk, []); % 每列是一个块 h = blk / 2; top = M(1:h, :); bot = M(h+1:end, :); M = [gate(1,1)*top + gate(1,2)*bot; gate(2,1)*top + gate(2,2)*bot]; psi = M(:); end

逻辑说明：reshape(psi, blk, [])把一维态向量摊成矩阵，列数由总维数自动推出，这就是「按位切块」的向量化写法。top是目标位为 0 的全部分量，bot是目标位为 1 的全部分量，两行赋值正好对应gate的两行。整个过程只额外申请一个和态向量同量级的内存，n=14（16384 维）在普通笔记本上几十毫秒就完。参数target的含义和singleQubitGate完全一致，所以这两套实现可以直接互换做交叉验证。

注意：applySingleFast只对单比特门成立。受控门要按控制位的取值做掩码，写法是把bot里控制位为 0 的那部分保留不动，只改控制位为 1 的部分，代码会多十行左右，思路一样。

3.4 测量：概率归一与塌缩采样

测量这一步看着简单，实际是最容易埋隐性 bug 的地方。概率幅取模平方得到概率，累加成分布函数，抽一个随机数落在哪段就塌缩到哪个基态：

function [idx, psiOut] = measureAll(psi) % MEASUREALL 对全寄存器做一次计算基测量 % 返回 0 起的测量结果 idx，以及塌缩后的态 psiOut p = abs(psi).^2; p = p / sum(p); % 抹掉浮点累积误差，保证和为 1 r = rand(); cdf = cumsum(p); idx = find(cdf >= r, 1, 'first') - 1; if isempty(idx), idx = numel(p) - 1; end % 兜底：r 落在末尾误差区 psiOut = zeros(size(psi)); psiOut(idx+1) = 1; end

逻辑说明：p / sum(p)这一行不能省，因为经过十几个门之后概率和通常会偏离 1 到 1e-15 量级，cumsum的最后一项可能小于 1，rand()恰好取到接近 1 的值时find会返回空，程序直接抛索引错误——这就是典型的「偶发崩溃，复现不了」。参数idx返回的是 0 起下标，和 2.3 节的约定保持一致。要固定随机结果做回归测试，在调用前加一句rng(42)。

4. 算法跑通：Deutsch-Jozsa、Grover 与相位估计的 MATLAB 写法

引擎层齐了，算法层其实很短。这一章挑三个最有代表性的：一个测「判定能力」，一个测「振幅放大」，一个测「相位读取」。三个跑通，基本能确认整套 core 没写错。

4.1 Deutsch-Jozsa：一个 oracle 矩阵判定常量还是平衡

Deutsch-Jozsa 要回答的问题是：给定一个把 n 位映射到 1 位的函数 f，它要么是常量（全 0 或全 1），要么是平衡的（一半 0 一半 1），用一次查询判定。经典做法最坏要 2^(n−1)+1 次。量子版本的关键是构造 oracle 矩阵Uf: |x>|y> → |x>|y ⊕ f(x)>，在 MATLAB 里就是一个置换矩阵：

function U = oracleDJ(n, fvec) % ORACLEDJ 构造 Deutsch-Jozsa 的 oracle 算符 % n : 输入比特数；fvec : 长度 2^n 的 0/1 向量，fvec(x+1) 即 f(x) dim = 2^(n + 1); U = zeros(dim); for x = 0:2^n-1 for y = 0:1 s = x * 2 + y; % 输入寄存器在高位 t = x * 2 + bitxor(y, fvec(x+1)); % 辅助位异或 f(x) U(t+1, s+1) = 1; end end end

逻辑说明：s = x*2 + y是把「输入寄存器 x」和「辅助位 y」拼成一个整数，输入寄存器在高位，这和 2.3 节的 big-endian 一致。U(t+1, s+1) = 1表示|s>映射到|t>。每个(x, y)对恰好写一次，所以U自动是置换矩阵，天然酉。参数fvec用向量而不是函数句柄传，是为了避免在循环里反复调用匿名函数拖慢速度。

线路部分是标准的四步：辅助位翻到|1>、全部打 H、跑 oracle、前 n 位再打 H，然后看前 n 位全零的概率。

n = 3; fvec = zeros(1, 2^n); % 常量函数 f(x)=0 % fvec = mod(0:2^n-1, 2); % 换成平衡函数，改这一行即可 ops = {}; ops{end+1} = singleQubitGate([0 1; 1 0], n+1, n+1); % X 打在辅助位 for k = 1:n+1 ops{end+1} = singleQubitGate(hadamard(), n+1, k); end ops{end+1} = oracleDJ(n, fvec); for k = 1:n ops{end+1} = singleQubitGate(hadamard(), n+1, k); end psi = runCircuit(basisState(n+1, 1), ops); p0 = norm(psi(1:2))^2; % 前 n 位全为 0 的概率 fprintf('P(输入寄存器全 0) = %.6f\n', p0);

逻辑说明：basisState(n+1, 1)构造的是|0...01>，最后一位是辅助位，整数下标 1 正好对应这个态。psi(1:2)取的是前 n 位全零、辅助位任意的两个分量，它们的模平方和就是「输入寄存器测得全零」的概率。常量函数下这个值应该是 1，平衡函数下应该是 0，中间值说明比特序或者 oracle 写错了。这个判据比看单个概率幅可靠得多。

4.2 Grover：两行代码的振幅放大与迭代次数

Grover 的整个算法可以浓缩成两行：把目标项的相位翻过来，再对所有项关于均值做一次反射。搜索 N=2^n 个元素，迭代约π/4·√N次后目标项概率接近 1。

function psi = grover(n, marked, iters) % GROVER 无测量版本的 Grover 迭代 % marked : 目标项的 0 起下标；iters : 迭代次数 N = 2^n; psi = ones(N, 1) / sqrt(N); % 均匀叠加，等价于每位打 H for it = 1:iters psi(marked+1) = -psi(marked+1); % 1) 翻转目标项相位 m = mean(psi); % 2) 关于均值反射 psi = 2 * m * ones(N, 1) - psi; end end

逻辑说明：mean(psi)算的是所有分量的算术平均，也就是<psi|s>，第二行就是2|s><s|psi> − |psi>的逐元素写法，不需要真的构造|s><s|这个秩一矩阵。参数iters一般取floor(pi/4*sqrt(N))，取多了概率会掉下来，这是 Grover 的固有振荡特性，不是 bug。

n = 4; N = 2^n; marked = 9; iters = floor(pi/4 * sqrt(N)); psi = grover(n, marked, iters); [~, idx] = max(abs(psi).^2); fprintf('目标 %d，测得 %d，成功概率 %.4f\n', marked, idx-1, abs(psi(marked+1))^2);

n=4时 N=16，迭代 3 次，目标项概率能到 0.96 以上。如果这个数字明显偏低，八成是marked的 0 起下标被写成了 1 起，或者mean那行被误写成sum/N之外的形式。

4.3 量子傅里叶变换与相位估计：把相位读出来

QFT 在 MATLAB 里就是一个范德蒙矩阵，一行写完：

function U = qftMatrix(n) % QFTMATRIX 返回 n 比特的量子傅里叶变换矩阵 N = 2^n; j = (0:N-1)'; k = 0:N-1; U = exp(2i*pi*(j*k)/N) / sqrt(N); end

逻辑说明：j*k是列向量乘行向量，得到 N×N 的相位矩阵，直接指数化再加归一化因子。这个矩阵是对称的，所以逆变换就是conj(U)，不用求逆，省一次inv的开销。工程上更快的写法是fft(eye(N))/sqrt(N)，两者数值上一致。

相位估计的用法是：给一个酉算符 U 和它的一个特征向量|u>，把对应的特征相位 φ 读出来。下面用最简单的相位门做例子，U = diag(1, e^(2πiφ))，特征向量取|1>：

n = 4; % 精度寄存器的比特数 phi = 0.25; % 待估相位，4 位精度正好能整除 psi = zeros(2^(n+1), 1); psi(2) = 1; % |0>^n ⊗ |1> ops = {}; for k = 1:n ops{end+1} = singleQubitGate(hadamard(), n+1, k); end for k = 1:n P = diag([1, exp(2i*pi*phi*2^(n-k))]); % 受控 U^(2^(n-k)) ops{end+1} = controlledGate(P, n+1, k, n+1); end ops{end+1} = kron(conj(qftMatrix(n)), eye(2)); % 逆 QFT 只作用前 n 位 psi = runCircuit(psi, ops); [~, idx] = max(abs(psi).^2); m = floor((idx-1) / 2); % 去掉末位辅助位 fprintf('估计相位 = %g（真值 %g）\n', m / 2^n, phi);

逻辑说明：第一个循环给精度寄存器打 H，第二个循环做受控演化，控制位是第 k 位、演化时长是2^(n-k)次方，这是 QPE 的标准结构。逆 QFT 用kron(conj(qftMatrix(n)), eye(2))只作用在前 n 位，末位辅助位保持不动。floor((idx-1)/2)是把「输入寄存器 + 辅助位」拼成的整数下标拆开，右移一位等价于除以 2。φ=0.25、n=4 时应该输出 4/16=0.25，一位不差。想测非整除的相位，比如 1/3，就得加大 n 看分布峰值的位置。

5. 避坑与排查：结果不对时先看这 5 个地方

算法跑出来不对，先别怀疑物理，百分之九十的坑集中在这五处。以上每一条都是我在不同项目里真撞过的，按「现象 → 原因 → 解决」记下来。

5.1 现象：恒定函数却以一定概率测出非零

Deutsch-Jozsa 用常量函数跑，p0应该严格等于 1，实际却出现 0.5 甚至 0.25 这类数字，而且每次跑都不一样。

原因：比特序不统一。最常见的是kron拼门时从低位往高位拼（for k = n:-1:1），但basisState和oracleDJ都按最高位在前算。两套约定混在一起，H 打在了错误的位上，此时结果既不报错也不完全随机，而是「有偏」。

解决：把 2.3 节那条断言放进tests/testCore.m，每次改完core就跑一遍。另外加一条更直接的：singleQubitGate(X, 3, 1)作用在|000>上应该得到下标为 4 的基态（二进制100），不是下标 1。这两条过了，比特序基本不会错。

5.2 现象：矩阵下标越界，或者基态总是偏移一位

Index exceeds matrix dimensions，或者明明写的是|00>却跑出|01>的结果，误差恒定位移 1。

原因：MATLAB 下标从 1 起，量子力学的基矢编号从 0 起。有人在算法层直接写psi(x)而不是psi(x+1)，或者bitset的位序没换算。

解决：只在basisState和measureAll这两个边界函数里做加减一，其他所有函数统一用 0 起整数下标。bitget(s, n-k+1)里的n-k+1是第二处必检点——bitget从最低位数起，第 k 位（最高位在前）对应的是第n-k+1位。这两处写对，位移问题就绝迹了。

5.3 现象：n 到 12 就 out of memory，或者跑得极慢

小比特数一切正常，比特数一加到 11、12，MATLAB 直接抛内存不足，或者跑到一半硬盘灯狂闪。

原因：把每个门都展开成 2^n × 2^n 的稠密复矩阵再存进ops。n=12 时单矩阵 512 MB，ops里存二十个就是 10 GB。而且大多数门其实只有一个 2×2 的实质内容。

解决：改用 3.3 节的applySingleFast，全程不构造大矩阵，内存只有态向量本身（n=14 才 256 KB）。受控门则用sparse：先speye(dim)起手，只改受控位为 1 的那些行，非零元数量级是 2^n 而不是 4^n。改完以后 n=14 的 Grover 在普通笔记本上几秒钟出结果。

5.4 现象：偶发索引错误，find返回空

measureAll偶尔抛Index exceeds，而且没法稳定复现，跑十次错一次。

原因：概率和没归一。经过十几个门之后，sum(abs(psi).^2)可能是 0.9999999999999998，cumsum的最后一项小于 1，而rand()恰好落在 1 附近，find(cdf >= r)就找不到任何下标。这不是玄学，是浮点累积误差的必然结果。

解决：p = p / sum(p)这一行必须留，再加一个isempty兜底。顺手把rng(42)固定住，让回归测试可复现。想更稳一点，可以在每次runCircuit结束后做一次psi = psi / norm(psi)，代价可忽略。

5.5 现象：相位估计的结果差一位，或者和真值互补

φ=0.25 应该输出 4/16，结果输出 3/16 或者 12/16。

原因：两处容易混。一是exp里的符号，exp(-2i*pi*...)和exp(+2i*pi*...)会得到互补的结果；二是U^(2^(n-k))的指数写反成2^(k-1)，这会得到一个位反转的答案，看起来像「二进制读反了」。

解决：固定用一个只有单个已知特征值的最小例子做单测——phi=0.25、n=4必须精确输出 0.25。改任何一行受控演化代码后，先跑这个例子。位反转的问题可以用一个自检排掉：把 φ 换成 0.5（二进制 1000），位反转后还是 0.5，看不出问题；换成 0.0625（二进制 0001），位反转后会变成 0.5，一眼就能认出来。

6. 二次开发：加一个自定义门，再把 12 比特跑进内存

改这套代码最常见的需求就两个：加自己的门、把比特数往上顶。

加自定义门有两条路。如果它是单比特门，写一个返回 2×2 矩阵的函数就够了，比如相位门：

function Rz = rzGate(theta) % RZGATE 绕 Z 轴旋转 theta 弧度 Rz = [exp(-1i*theta/2), 0; 0, exp(1i*theta/2)]; end

然后singleQubitGate(rzGate(pi/4), n, k)或applySingleFast(psi, rzGate(pi/4), n, k)直接就能用，不用改core。如果它是三比特以上的多比特门，别去手工填 2^n 维矩阵，用kron拼出小的 8×8 或 16×16 矩阵再扩维，或者拿controlledGate嵌套两层。

6.1 三条断言，改完必跑

我自己改这套代码时会固定跑三条断言，全过了才敢往下走：

% 1) H 连打两次等于单位阵 psi = basisState(3, 5); H = hadamard(); for k = 1:3, psi = applySingleFast(psi, H, 3, k); end for k = 1:3, psi = applySingleFast(psi, H, 3, k); end assert(norm(psi - basisState(3, 5)) < 1e-12, 'H 不是自逆'); % 2) 两套实现完全等价（稠密 vs 分块） a = singleQubitGate(rzGate(0.7), 4, 2) * basisState(4, 9); b = applySingleFast(basisState(4, 9), rzGate(0.7), 4, 2); assert(norm(a - b) < 1e-12, '分块路径和矩阵路径不一致'); % 3) Grover 在 n=2 时一次迭代成功率应为 1 psi = grover(2, 2, 1); assert(abs(psi(3))^2 > 1 - 1e-12, 'Grover 单次迭代未收敛');

第二条最值钱，它同时验证了比特序和门的展开顺序——只要两边有一处约定不一致，norm立刻就不为零。第三条用来卡 Grover 的迭代次数和marked下标。

6.2 上比特数：稀疏与分块的耗时对比

把core换成稀疏加量以后，同一台机器上的大致量级差别是这样：

差别的来源不是算法，是内存带宽。稠密路径每次矩阵向量乘都要把 512 MB 读一遍，分块路径只碰态向量本身。

从那以后我每次动core里任何一个造门函数，都强制先把那三条断言跑一遍，再去看算法层的输出——因为算法层出错的样子太好看了，概率分布依旧平滑、图依旧能画，只是数字全错，靠眼睛根本认不出来。这套代码不大，但它把一个道理讲得很清楚：量子算法在 MATLAB 里没有魔法，全是下标、张量积和归一化，把这三件事钉死，剩下的就是照着公式写。这份资源我放在下面，想上手的话直接拿demo.m跑一遍，再挑一个算法拆开改改门，比自己从零搭一套快得多。希望帮到你。

本文还有配套的精品资源，点击获取