企业官网建设流程全解析

1. 为什么“同样一份数学代码”，别人的库能快十倍

做数学库这件事，最早我是很不以为然的。上学那会儿写矩阵乘法，三层循环套上去，能跑出结果就算完事。直到第一次把一份 N=2048 的矩阵乘法丢进项目里，发现耗时整整两秒多，而同门的程序用同样的算法只花了不到一百毫秒。那一刻我才意识到一个问题：一个数学库的性能，和你“会用公式”之间，隔着一整座冰山。

1.1 先厘清“高性能数学库”到底在解决什么问题

所谓高性能数学库，它不是指把 sin、cos、sqrt 这些函数重新发明一遍，而是指在一个特定硬件平台上，把数学运算的吞吐量、延迟、精度和资源占用压到极致。常见的形态包括 BLAS 级别的向量与矩阵运算、LAPACK 级别的线性方程组求解、FFT、稀疏矩阵运算，以及统计与随机数生成。它们共同的特点是：计算模式高度规则、数据规模大、可并行性强，因此非常吃“工程优化”，而不只是算法复杂度。

拿矩阵乘法举例，经典的 i-k-j 三重循环，时间复杂度都是 O(n^3)，这是算法层面的天花板。但同样的 O(n^3)，不同写法在真实机器上的耗时差距能拉到 20 倍以上。原因是：现代 CPU 的算力远远超过内存系统能供给的数据速度，你的代码能不能让计算单元一直有数据可算，决定了最终效率。高性能数学库的本质工作，就是“喂养”计算单元，而不是炫数学技巧。

1.2 性能瓶颈的第一性原理：访存 > 计算

很多刚接触优化的同学会把注意力放在“减少乘法次数”上，比如用斯特拉森算法。但真实项目中，斯特拉森很少有人直接用，原因很简单：它的常数因子大、数值稳定性差，更重要的是，对现代 CPU 来说，乘法本身已经是极廉价的指令了，真正贵的是把数据从内存搬进寄存器的那一下。一次 L1 cache miss 大概是 4 个时钟周期，L2 miss 约 12 个，L3 miss 约 40 个，而一次主存访问要 200 个时钟周期以上。相比之下，一次浮点乘法在 Haswell 以后的 CPU 上只需要 0.5 个周期左右的吞吐。算一下就知道，省一个乘法指令，远不如省一次 cache miss 来得实在。

所以我在优化任何数学函数之前，第一件事就是计算访存强度，也就是“每访问一次内存，能支撑多少次浮点运算”。对矩阵乘法这种访存密集型任务，如果能把数据块塞进 L2/L1 cache，让同一块数据被反复复用，性能自然就上来了。这个思路贯穿了后面所有优化步骤。

1.3 一次真实对比：三种矩阵乘法实现的时间差异

为了让你快速感知差距，我放一张我在自己机器上测过的数据（Intel i5-12400，单线程，N=1024）：

同一个算法，同一个精度（float），只是让数据的流动方式和计算指令的组织方式不同，结果差了 34 倍。这不是魔法，而是下面每一层优化叠加起来的效果。接下来我会把每一步拆开讲，并且把为什么这样做、踩过什么坑都交代清楚。

2. 数据布局与内存子系统：数学库最快的那条腿

如果你只记住一句话，那就是：数学库性能优化的第一步是数据布局，不是 CPU 指令。很多人在这一层还没做好就急着上 AVX、上多线程，结果收益很有限，因为瓶颈在内存搬运，不在计算。

2.1 行主序与列主序：一次cache miss就足以改变命运

在 C/C++ 里，二维数组默认是行主序存储，也就是 a[i][j] 在内存里是连续的一行一行排下来的。而 Fortran 和 MATLAB 的默认存储是列主序，a[i][j] 的 i 变化最快，j 变化最慢。

这带来一个直接影响：访问 a[i][k] 时按行读是连续的、高速的；但如果你在三重循环里对 b[k][j] 按列去读，那每一次读取都大概率是 cache miss。我最早写朴素矩阵乘法的时候就是这样，内层循环对 b 做列访问，每一列都跨越整个矩阵的内存范围，几乎全 miss。这 310ms 就是这么来的。

解决办法有两个方向：一是把 b 矩阵转置成行主序的临时矩阵，二是交换循环顺序，让内层循环对 a、b 都是行访问。我实际项目中两个都用，优先用转置临时矩阵，因为它不需要改动计算逻辑，只是多一次 O(n^2) 的拷贝，对 n=1024 来说，这个拷贝成本远低于列访问多出来的 cache miss。

2.2 AoS/SoA选择：线性代数库如何为SIMD铺路

矩阵是天生适合行主序的，但很多数学库不只是处理矩阵，还需要处理“一批点”的坐标变换、批量向量的运算。这时候就涉及到 AoS（Array of Structures，结构体数组）和 SoA（Structure of Arrays，数组结构体）的选择。

举个例子，处理 10000 个三维点的旋转，如果定义成：

struct Point { float x, y, z; }; Point points[10000];

这样声明，points[0].x、points[0].y、points[0].z 在内存里是紧挨着的。你要对每个坐标做乘法，访问模式是“隔着 12 字节取一个 float”，这对 SIMD 非常不友好，想要一次加载 4 个 float，必须做 gather 或者手动拆包。

改成 SoA 就清爽了：

struct Points { float x[10000]; float y[10000]; float z[10000]; };

这样 x 数组是一整块连续内存，一次可以 load 4 个甚至 8 个 float，直接进入向量寄存器做运算，效率提升非常明显。代价是代码可读性略差，抽象层级要设计好。我在自己的小库里对外暴露面向对象的接口，内部存储全部用 SoA。

2.3 对齐分配与跨平台陷阱

数据布局定了，还有一个容易忽略的细节：内存对齐。

AVX2 一次可以加载 32 字节，AVX-512 是 64 字节。如果你用 malloc 动态分配一个 float 数组，它默认只保证 16 字节对齐，运气好落在 32 字节边界，运气不好就落到一个尴尬的位置。当未对齐地址传给对齐 load 指令时，轻则性能回退，重则直接段错误。

我吃过一次亏：在调试版里一切正常，开了 -O2 用了对齐 load intrinsic 后，程序随机崩。后来用 addr2line 才发现是数组首地址没对齐。解决方式有三种：

• posix_memalign或_aligned_malloc，明确分配对齐内存。
• C++17 的std::aligned_alloc，注意它要求 size 是对齐值的整数倍，容易踩坑。
• 或者更简单，声明一个 64 字节对齐的结构体数组，内部放一个足够大的alignas(64) float data[]。

对齐这件事，同一个代码在不同编译器下的行为也会不一样，尤其是 MSVC 的__declspec(align)和 GCC 的__attribute__((aligned))写法不同，封装一层宏是值得的。

2.4 Bandwidth vs Compute Intensity：理解Roofline模型给优化划边界

“Roofline”听起来很高深，说白了就是画两条线：一条叫内存带宽上限，一条叫计算峰值上限。横轴是计算强度，也就是每读入一个字节数据能做多少次浮点运算；纵轴是实际能达到的性能。如果你的计算强度很低，性能就会被带宽压着走，再怎么优化计算指令都没用，因为数据喂不进来。只有把计算强度推进到超过某个临界点，性能才会跑到计算上限那一段。

做数学库优化，我习惯在动手前先算这个临界点。例如，我的机器内存带宽实测约 30GB/s，AVX2 单线程浮点峰值约 50 GFLOPs，那么临界计算强度就是 50/30，约等于 1.7 FLOP/Byte。如果某个运算的计算强度不足 1.7，优化的重心应该放在缓存友好和数据复用上；如果远超 1.7，那就放心去压榨计算指令和 ILP。

这解释了一个常见困惑：为什么同样一台机器，矩阵乘法能做得很高效，而一个简单的逐元素加和却永远跑不上高 GFLOPs。逐元素加和每个元素只读一次，写一次，几乎没有复用，计算强度极低，上限就是内存带宽。这类运算，再牛的库也只是把带宽用到极限而已。

3. 让编译器把你的循环变成向量：微架构级优化实操

数据布局搞定之后，下一步就是让你的循环尽量贴近 CPU 的向量执行单元。这一步是“榨干硬件”的关键，也是最容易写出“看着优化了实际上更慢”代码的地方。

3.1 自动矢量化：编译器不是万能的

很多编译器确实有自动向量化能力，GCC 的 -O3、Clang 的 -O3 以及 Intel 编译器的 -xHost 都能把简单循环向量化。但这里有个前提：循环结构必须足够简单、内存访问必须连续、循环次数最好能在编译期或运行期确定且不被复杂的别名分析拦下。

我遇到最典型的自动向量化失败场景是这样的：

void vector_add(const float* a, const float* b, float* c, int n) { for (int i = 0; i < n; i++) { c[i] = a[i] + b[i]; } }

逻辑再简单不过了，可如果调用方传入了c和a指向同一块内存的别名，编译器为了安全，就没法放心重排指令，向量化就直接放弃。解决办法是加__restrict__告诉编译器这些指针没有重叠，或者提前拷贝到临时缓冲区。GCC 的-fopt-info-vec能输出向量化报告，建议读者开起来看一眼，很多时候你会发现编译器“卡”在哪儿。

3.2 依赖链与延迟：为什么你的累加循环跑不满

还有一个常被忽略的坑：循环携带依赖。看这段代码：

float sum = 0.0f; for (int i = 0; i < n; i++) { sum += a[i]; }

每次sum += a[i]都依赖上一次sum的值，下一次加法必须等上一次的 FADD 指令完成。FADD 的延迟通常是 4 个时钟周期左右，这意味着这个循环的吞吐被限制在每 4 个周期加一个元素，和 CPU 实际上每周期能做 2 次浮点加法的能力相去甚远。

解决办法是“多路累加”，把累加器拆成 4 个独立的变量，每次循环交替加：

float sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0; for (int i = 0; i < n; i += 4) { sum0 += a[i]; sum1 += a[i+1]; sum2 += a[i+2]; sum3 += a[i+3]; } float sum = (sum0 + sum1) + (sum2 + sum3);

这样四条依赖链并行执行，吞吐立刻翻 3-4 倍。更优雅的写法是让编译器在自动向量化时自己展开，但显式写多路累加在大多数编译器上表现更稳定。拆成 4 路还是 8 路，取决于目标微架构拥有几个加法执行单元，一般 4 到 8 路比较合适，再多也不会更快，反而容易寄存器溢出。

3.3 FMA、reduce与分支消除：一个累加器的完整改造

真正到库级别，累加还会用到 FMA（Fused Multiply-Add），也就是一条指令同时完成乘法和加法，既省时间又减少中间舍入误差。对多项式求值、点积这类操作，FMA 几乎是性能神器。

看一个多项式求值的例子。普通写法是霍纳法，逐项累乘加，依赖链很长。优化后我可以把它拆成偶次项和奇次项两条独立链，分别用 FMA 计算，最后合并。这样依赖链深度减半，吞吐自然上去。不过要注意，FMA 改变了浮点运算的舍入顺序，如果你的库对结果精度有逐位匹配的硬性要求，这里需要做版本控制。

另外，能不用分支就不用分支。现代 CPU 的分支预测在简单规则下很准，但面对数据相关的分支（比如if (x > 0) y = 1; else y = 0;），预测失败代价大概是 15-20 个周期。这种场景可以直接用比较指令生成掩码，再和常数做位运算，做到无分支。数学库里的符号函数、绝对值、clamp 都可以这样优化。

3.4 编译优化选项的组合拳

单文件里写了很多优化代码，但如果编译选项不对，前面的辛苦全部白费。我常用的组合是：

• GCC/Clang：-O3 -march=native -ffast-math -fno-math-errno -fomit-frame-pointer -funroll-loops
• 用-ffast-math之前要想清楚，它会假设 NaNs/Infs 不会出现、允许重排浮点运算、关闭 errno 设置，对安全性要求高的代码慎用。我自己一般会单独为内部性能敏感模块开，对外 API 保持保守。
• 大项目里开-flto（链接时优化）收益也很大，它能让跨编译单元的常量传播和内联共享变成可能，但会显著增加编译时间和内存占用。

我在实际项目中还发现一个细节：-march=native会启用编译器对当前 CPU 的自动检测，如果你把编译好的库放到另一台老机器上，可能会出现 illegal instruction 崩溃。分发库的时候应该用-march=x86-64-v3这类更保守的指令集级别，而不是图省事直接 native。

4. 多线程并行：线程数不等于性能

单核优化做到位之后，接下来必然是并行。但并行不是“开几个线程然后等奇迹”，这一层踩坑最多，也最容易做出负优化。

4.1 两种并行切分策略：数据分片 vs 管道分片

数学库的并行通常用数据分片，因为每个数据元素的计算是独立的，天然适合把矩阵按行块或列块切给不同线程。管道分片在数学库里用得少，因为每一步的输出要被下一步消费，线程间经常需要同步，同步成本很容易吃掉并行收益。

举个例子，矩阵乘法把输出矩阵的 i 维按行切成 4 段，每段交给一个线程去算对应的行块，这是一种典型的数据分片。注意切分要尽量让每段运算量相等，避免一个线程算得飞快，另一个还在跑长尾。如果矩阵规模不确定，最简单的办法是动态调度，比如用 OpenMP 的 schedule(dynamic, chunk_size)。

4.2 伪共享：你不敢信的“负优化”案例

伪共享是我见过最隐蔽的并行性能杀手。它的原理是：CPU 的缓存一致性是以 cache line（通常 64 字节）为单位的。如果两个线程分别修改不同变量，但这两个变量凑巧落在同一个 cache line 里，缓存一致性协议会让两个核心不断互相通知“这个缓存行被改了”，导致明明互不相干的写入变成串行化，性能暴跌。

我做过一个实验，开 4 个线程，每个线程对数组的某一段独立求和，然后把结果写进四个不同的全局变量。仅仅因为这四个 float 变量紧挨着存放在一起，性能比把它们隔开 64 字节对齐慢了三倍多。解决办法是给不同线程的结果变量加 padding，让它们落在不同 cache line：

struct alignas(64) PerThreadResult { float value; char padding[60]; };

这样就把每个线程的写入隔离到独立缓存行。类似的还有循环切块大小，最好也按 cache line 的整数倍去切，避免两个线程把同一行数据拆开处理。

4.3 时间分片与调度开销：OpenMP schedule的选择

用 OpenMP 时，schedule 的选择直接影响负载均衡。static调度是编译期均分，开销最小，但如果数据分布导致不同区域计算量不等，线程会空等。dynamic调度是运行时按块领取，能自动均衡，但每领一块都有锁开销。guided是介于两者之间，块大小递减，适合循环内部计算量差异大的情况。

我自己的经验是：如果循环内计算量均匀，优先static；如果计算量难以预估，选guided而非dynamic，因为它的调度次数更少。有一回矩阵维度正好让 8 线程的 static 切分把一个线程分到几乎为空的行区间，另一个线程却要算整个最重的列块，性能直接掉了一半。换成 guided 之后才恢复正常。

4.4 同步原语与原子操作在数学库里合理边界

数学库内部的同步越少越好。多线程累加场景里，最常见的错误是直接把sum += local_value放进临界区，或者频繁使用原子加。这类操作在低竞争下勉强能用，高竞争时锁和原子操作的开销甚至比计算本身还高。

正确的做法是每线程维护局部结果，最后用一个 O(n_threads) 的归约合并结果。这个归约阶段即使加锁，竞争也只发生一次，开销忽略不计。如果并发级别不高，也可以直接用无锁的原子加做最终归约，但要注意原子操作对 cache line 的影响，和上面伪共享是同一个问题。

5. GEMM三层优化实战：从朴素实现到接近硬件峰值

泛泛谈优化容易陷入空洞，下面我用最经典的 GEMM（通用矩阵乘法，C = A * B）完整走一遍优化流程。GEMM 是整个 BLAS 体系的基石，也是几乎所有线性代数库评测里的重头戏。

5.1 版本0：朴素三重循环，性能基准

void matmul_naive(const float* a, const float* b, float* c, int n) { for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { float acc = 0.0f; for (int k = 0; k < n; ++k) { acc += a[i * n + k] * b[k * n + j]; } c[i * n + j] = acc; } } }

这个版本的问题前面已经讲过：b 是按列访问的，内层循环每次b[k][j]都在跳跃地址，cache miss 率极高。在 n=1024 时，实测单线程耗时 310ms 左右。作为基准，它的价值是告诉我们在没有任何优化时的代价。

5.2 版本1：转置+缓存分块+向量化

第一步，先把 b 转置一遍，让内层循环里两个矩阵都变成按行访问：

void matmul_t1(const float* a, const float* b, float* c, int n) { float* bt = new float[n * n]; for (int i = 0; i < n; ++i) for (int j = 0; j < n; ++j) bt[j * n + i] = b[i * n + j]; for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { float acc = 0.0f; for (int k = 0; k < n; ++k) { acc += a[i * n + k] * bt[j * n + k]; } c[i * n + j] = acc; } } delete[] bt; }

这一步已经能降到 76ms 左右。但还不够，因为每次内层循环取 a 的一行和 bt 的一行，这两行都来自主存，没有充分利用 cache。下一步是分块，让计算的最小单元在一个块内反复复用 L1/L2 缓存：

void matmul_t2(const float* a, const float* b, float* c, int n, int block) { for (int i0 = 0; i0 < n; i0 += block) { for (int j0 = 0; j0 < n; j0 += block) { for (int k0 = 0; k0 < n; k0 += block) { for (int i = i0; i < i0 + block; ++i) { for (int k = k0; k < k0 + block; ++k) { float a_tmp = a[i * n + k]; for (int j = j0; j < j0 + block; ++j) { c[i * n + j] += a_tmp * b[k * n + j]; } } } } } } }

block 大小一般选 32 或 64。选小了块太小，每次加载的 a、b 片段复用次数不够；选大了，块内数据超过 L2 容量，又会频繁换出。这个版本在 32 块大小时，耗时能再降到 30ms 上下。最后叠加 AVX2 向量化，内层 j 循环一次算 8 个 float，耗时进一步降到 15ms 左右。

5.3 版本2：多线程+进一步调优

把分块后的 i0 层直接抛给 OpenMP：

#pragma omp parallel for schedule(guided) for (int i0 = 0; i0 < n; i0 += block) { ... }

注意我选择对 i0 并行，而不是内层 j0，因为并行粒度越粗，线程切分和合并的开销就越小。在 4 线程实测下，这个版本能跑到 9ms 左右。

进一步调优还有两个方向：一是把 block 的尺寸按矩阵维度动态调整，二是让内层循环用restrict修饰指针、避免别名误判。还有一个容易忽视的点：c[i * n + j]在内层循环里被反复读写，如果 block 不太大，c 的对应片段会留在 L2 里，但为了保险，可以为每个输出块开一个 block×block 的临时累加数组，算完再写回。这样 c 的读改写从 main memory 变成了 cache 内的局部操作，效率更高。

5.4 性能数据对比与提升路径检验

把这几个版本放一起，就能看到每一层优化的贡献：

从数据看，转置和分块带来的收益甚至超过向量化。这验证了之前说的“访存优先”原则：先解决数据怎么流动，再考虑计算指令怎么排。向量化的收益是在缓存优化做完之后才显出来，如果一开始就上 AVX 而没有分块，效果会大打折扣。

6. 验证与基准测试：跑得快的前提是算得对

优化做得越多，正确性验证就越不能省。尤其是做数学库，服务对象可能是别的工程师在使用，一个错误的结果比慢十倍还要糟糕。

6.1 测试数据设计：三组必测用例

我一般会设计三组测试：

第一组是“已知答案”用例。比如让矩阵元素全是 1，或者单位矩阵乘任意矩阵，结果必定一眼能看出对错。这类用例跑一遍，基本功能错误立刻暴露。

第二组是“随机数据 + 高精度参考”。用 double 甚至 long double 实现一个相同运算的朴素版本，作为参考，然后对 float 的优化版本做相对误差对比。随机数据能覆盖到各种数值区间，避免某个特定输入下精度刚好够好。

第三组是“边界条件”。包括 n=0、n=1、n=非 4 的倍数、NaN/Inf 输入、非常大和非常小的数值。很多优化代码在 n 不是 SIMD 宽度倍数时，会写出错误的后处理逻辑，这组用例就是抓这个的。

6.2 ULP误差控制：浮点数学库的“误差牌照”

数学库和普通业务代码不同，它不能只说“结果大概差不多”，而是要以 ULP（Unit in the Last Place）为单位定义误差上限。ULP 是指两个相邻浮点数之间的间隔，1 ULP 误差就是你结果的最后一位可能有偏差。

对不同函数，业界对不同函数有不同要求：加减乘除这类四则运算，IEEE 754 要求正确舍入，误差不超过 0.5 ULP；exp、log、sin、cos 这类超越函数，通常允许 1-4 ULP。如果你的库要跟 BLAS/LAPACK 对标，误差也得定在相近水平。

所以我在优化每个函数时，都会先写一个“慢但正确”的参考实现，再用随机输入批量对比，要求误差在预设 ULP 范围内。这里有个经验：浮点优化里常见的错误不是算错，而是重排运算顺序后误差被放大。比如前面提到的多路累加和 FMA，改变累加顺序会改变舍入误差，虽然通常在一个 ULP 以内，但如果库的接口契约里写死了“逐位匹配”，就需要专门跑全量黄金样本。

6.3 基准测试方法论：高频抖动、预热与统计口径

基准测试也不能随便跑一遍就看数字。现代 CPU 的频率是动态的，同一个循环，冷启动和热启动能差 30% 以上。我做基准的标准流程是：

1. 用taskset或numactl把进程绑到固定核，避免操作系统调度导致核心漂移。
2. 先在循环外跑一轮预热，把缓存填起来，触发频率拉升。
3. 正式测量时，每个尺寸跑 5-10 次，取中位数而不是平均值。中位数能有效剔除系统中断和偶发抖动的影响。
4. 同时记录功耗和温度，防止机器过热降频影响结论。

另外，我刚入行时犯过一个错：优化完一个函数后，随手用系统自带的 chrono 计时，结果把打印、内存分配也一起算进去了，导致“优化效果”被其他开销稀释。正确的做法是只对被测试的核心计算循环计时，内存分配、初始化、打印都放到计时区间之外。

7. 真实项目的额外复杂度：不是只有矩阵乘法

GEMM 是最典型的例子，但一个真正可用的数学库，还要面对大量“不那么规整”的问题。这些才是拉开库之间差距的地方。

7.1 特殊函数、向量运算与数值稳定性

比如指数、对数、三角函数这些特殊函数。它们的优化不完全是“算得快”，而是“在快的同时保持数值稳定”。比方说计算log(1+x)，如果直接实现log(1.0 + x)，当 x 很小时，1.0 + x 会先被舍入到 1.0，导致结果直接变成 0。正确做法是用专门的log1p展开式，对很小和很大的 x 分别处理。

同样地，exp(x)在 x 接近边界时要防止溢出到 Inf，tan(x)在接近π/2时会爆炸，这些边界行为必须仔细定义。数学库的每一个函数背后都有一本“数值食谱”，有的来自经典参考书，有的来自平台原厂库的对照。

7.2 稀疏矩阵与数据结构取舍

真实场景里，很多矩阵是稀疏的，比如图计算、有限元模拟中的矩阵，非零元素占比可能只有 1%。这时如果你还用稠密矩阵的存储方式去优化，即使 GEMM 再快也是巨大的浪费。

稀疏矩阵的优化重点变成“只访问非零元素”。常用的存储格式有 CSR（Compressed Sparse Row）和 CSC（Compressed Sparse Column）。CSR 用一个数组存所有非零值，再用两个数组记录每一行的起始位置和每行非零元素的列号。这种方式省内存，但访问模式变得不规则，很难像稠密矩阵那样用 SIMD 吃到收益。

我做过一个稀疏矩阵向量乘，普通 CSR 实现遇到 GPU 上的分支发散严重，换了 ELLPACK 格式（把最多非零元素的那行定成宽度，其他行补零）之后，SIMD 友好度大幅提升，速度翻了 3 倍。代价是内存略多一点，非零分布极不均匀时会浪费严重。

7.3 精度与性能的矛盾：什么时候该用float，什么时候必须double

很多人以为 float 一定比 double 快，这是误解。在支持双精度 FMA 的芯片上，double 的吞吐和 float 差不多（比如某些服务器 CPU 的 FP64 单元是完整宽度的），只是内存带宽消耗减半。而有些消费级 CPU 的 FP64 单元是被阉割的，跑 double 就明显慢。所以“用 float 还是 double”取决于目标硬件和数值需求。

我的经验是：如果业务数值范围很大、需要高精度累加，先用 double 做原型，确认误差水平后再决定是否下沉到 float。如果确实需要 float 性能，但担心累加误差，可以考虑用 Kahan 补偿求和——它只增加少量指令，却能把累加误差降到一个很低的水平。

8. 我在多次重写数学库后最想分享的几点经验

最后说点个人感受，可能比前面的技术细节更有价值。

第一，优化是迭代式的，不是一蹴而就。不要试图“一次写出最优代码”，而是先写好正确版本，再逐步按访存、指令级并行、多线程的顺序一层层加优化。每做一步就测一次，既能定位瓶颈，也能防止后面的改动破坏前面的成果。

第二，用性能剖析器验证你的直觉，不要靠猜。我经常打开perf stat看 cache-misses、branch-misses、IPC 这些指标。有一次我以为瓶颈是分支预测失败，数据出来却是 L1 cache miss 占了大头，方向完全反了。

第三，版本管理里记录每个 commit 的性能增量。数学库的重构风险极高，没有性能回归记录的优化很容易在几个月后被一个看似无害的改动悄悄毁掉。我会在 CI 里挂一个简单的性能基准作业，每次 merge 前自动跑一遍关键 case，超过阈值就报警。

第四，不要迷信“某个库就是最快”的传说。每个库都有自己的适用场景。MKL 在 Intel 平台上是标杆，但换到 AMD 或者 ARM 上未必最优。真正可靠的做法，是拿你自己的业务矩阵尺寸和访问模式去实测，跑出一个可复现的对比数据再下结论。

数学库这条路，做久了你会发现，真正难的不是数学公式，而是让公式在真实的硬件上跑得又快又稳。每一次优化，都是和硬件的一次对话。希望这篇文章能帮你把第一句“对话”说清楚。