第四节：GPU/CPU浮点运算的差异（加法、乘法、融合乘加）
说白了，GPU 和 CPU 的浮点运算结果不一样，这几乎是所有做异构计算的人都会遇到的坑。我记得刚接触 CUDA 那会儿，总觉得都是 x86 的机器，GPU 也是当代硬件，算个浮点数还能差到哪去？结果项目联调的时候，脚本一跑，数值对不上，排查了整整一周才找到根因——就是浮点运算的中间精度和归约顺序搞的鬼。
核心结论先摆出来：CPU 和 GPU 在浮点运算上的差异，不是“谁对谁错”的问题，而是“谁怎么算”的问题。标准一样（IEEE 754），但实现细节不一样，结果就会有微小偏差。
4.1 浮点运算标准——大家读同一本教材
CPU 和 GPU 都遵循 IEEE 754 标准。单精度（float）32 位，双精度（double）64 位。格式一样，但中间计算精度和舍入方式不同。
给你看个表，一目了然：
| 特性 | CPU（x86/x64） | GPU（NVIDIA/AMD） |
|---|---|---|
| 单精度（float） | 1 sign + 8 exp + 23 mantissa | 1 sign + 8 exp + 23 mantissa |
| 双精度（double） | 1 sign + 11 exp + 52 mantissa | 1 sign + 11 exp + 52 mantissa |
| 中间精度 | 默认 80 位扩展精度（x87），SSE/AVX 可配置 | 通常与目标精度一致（float 运算按 float 算） |
| FMA | 可启用/禁用（编译器选项控制） | 默认启用，硬件原生支持 |
| 归约顺序 | 确定性强（串行加） | 不确定性强（并行树形加） |
数据格式一样，这没毛病。但中间精度一不一样，直接决定了最终结果的差异。我习惯把浮点运算的差异分成三类来说：
- 一元运算：加法、减法，误差主要来自舍入
- 二元运算：乘除，误差来自尾数乘法和舍入
- 混合运算：融合乘加（FMA），这个是重灾区
4.2 加法——看似简单，细节不少
两个浮点数相加，CPU 和 GPU 的套路是一样的：对阶->尾数相加->规格化->舍入。但微妙的地方在于：
- 中间精度不一样。x87 的 FPU 默认用 80 位中间寄存器，加了两个 float，结果先存在 80 位里，最后再截断回 float。GPU 呢？大多数 GPU 直接按 32 位算完就截断。差就差在这一步。
- 加法顺序影响结果。1e8 + 1e-8 + (-1e8) 这个经典例子，CPU 按顺序算还是 GPU 并行算，结果完全不同。
我曾经在排查一个气候模型的联调问题时，发现 CPU 和 GPU 的差异达到 1e-5 量级，逐拍对比后发现就是 80 位中间精度导致最后一位不同。解决方案很简单：让 CPU 也用 /fp:precise 控制中间精度，或者干脆强制都用 double 算。
4.3 乘法——尾数乘法的精度差异
浮点乘法看起来没什么花头：指数相加，尾数相乘，然后规格化+舍入。但问题出在尾数相乘的精度上。
CPU 的乘法器实现很成熟，特别是 AVX-512 里的 FMA 单元，精度是有保证的。GPU 的乘法器面积受限，有些老架构的 GPU 乘法尾数截断会丢掉 1-2 个 ULP（Unit in the Last Place）。
举个例子，假设你跑个简单的程序：
float a = 1.23456789f;
float b = 8.76543210f;
float c = a * b;CPU 上可能算出来是 10.823097...，GPU 上可能是 10.823096...。最后一位差 1，这就是 ULP 差异。
千万别小看 1 个 ULP 的差异。在迭代数值算法（如求解偏微分方程）中，每一步多一个 ULP 误差，1000 步迭代后误差就可能放大到不可接受。我见过一个金融风险评估的项目，就因为 GPU 单精度乘法比 CPU 多偏了 2 个 ULP，跑了 10 万次蒙特卡洛后结果整体漂移了 0.01%。
4.4 融合乘加（FMA）——最大的差异来源
FMA，全称 Fused Multiply-Add，就是一条指令完成 a * b + c。
这里面的关键点：FMA 不截断中间乘积的尾数，整个乘加过程只做一次舍入。而普通乘法+加法，乘法结果先舍入一次，加法结果再舍入一次，两轮舍入，精度自然不如 FMA。
来看组对比：
| 表达式 | CPU（不启用 FMA） | GPU（默认 FMA） | 绝对差异 |
|---|---|---|---|
| a*b+c | mul 舍入 -> add 舍入 | FMA 一次舍入 | 最多 1 ULP |
| a*b-c*d | 两次 mul + 一次 add | 两次 FMA | 可达 2 ULP |
| 点积 a·b | 串行乘加，多次舍入 | 树形 FMA，低舍入次数 | 随长度增加 |
你想想看，如果一个矩阵乘法的核心循环里，GPU 用了 FMA 而 CPU 没启用，那最后的结果差异会在哪里？说白了就是 每个浮点运算的结果都差那么一点点，累计起来差异就明显了。
嗯，这里要注意：启用 FMA 并不意味着精度一定更高或更低。FMA 舍入次数少，但在某些场景下反而会因为“更精确”而产生和 CPU 结果不同的数值路径变化，导致结果偏离。
经验之谈：我建议在 GPU 一致性测试中，把 FMA 当作一个独立的变量来控制。对比测试时至少跑三组：CPU（无 FMA）、GPU（FMA 打开）、CPU（FMA 打开），这样才能明确差异到底是谁产生的。
4.5 编译器优化带来的差异
说到编译器，其实 CPU 端也得背锅。同样的 C++ 代码，gcc 和 MSVC 编译出来的浮点运算顺序都可能不同。比如：
float sum = a + b + c + d;编译器会把它规约成 ((a+b)+c)+d 还是 (a+b)+(c+d)？没有标准答案。不同优化等级下，CPU 端的偏差就已经存在了。
GPU 端呢？NVCC 编译器默认会做大量的指令重排和 FMA 融合。尤其是 -fmad=true 选项，几乎是默认打开的。如果不加控制，GPU 代码里的乘加运算会被大量替换成 FMA 指令。
怎么控制？
// NVCC 禁用 FMA
// 编译时加 --fmad=false
// MSVC 禁用 FMA
// #pragma float_control(precise, on)
// GCC/Clang 禁用 FMA
// -ffp-contract=off我在项目中遇到过这样一个 case：CPU 端用 MSVC 的 /fp:fast 编译，GPU 端默认开 FMA，结果两个结果差异大到 1e-3 量级。后来把 CPU 端改成 /fp:strict 并 在 GPU 端加上 --fmad=false，差异才收敛到 1e-6 以内。所以，做一致性检查时，两边编译器选项必须对齐。
4.6 总结差异点
把上面的内容浓缩成一张图，方便你快速回顾：
看完这张图你应该清楚了：加法的差异在中间精度和加法顺序，乘法的差异在尾数截断，FMA 的差异在舍入次数和编译器自动融合。
做一致性检查时，我会单独列一个清单，逐项确认：
- 中间精度：CPU 端是否用了 /fp:precise 或 /fp:strict？
- FMA 开关：两端是否一致？GPU 端加没加 --fmad=false？
- 编译器优化：-O2 还是 -O3？有没有自动向量化？
- 归约方式：GPU 的并行规约用了什么树？是否采用 Kahan 补偿求和？