Numba cuda.jit与njit结果不一致，如何统一计算精度？

1. 统一浮点舍入模式 #

• CPU端：用Numba的njit时直接设fastmath=False禁用快速数学优化，同时在函数外调用numba.set_float_rounding('round_nearest')锁定舍入模式，避免编译器自动优化改变舍入行为。
• GPU端：编译CUDA函数时同样加fastmath=False，CUDA默认是向最近舍入，但快速数学选项会乱改规则，必须关掉。别用--fmad=true这类自动融合乘加的编译参数，除非你主动用FMA统一计算逻辑。

2. 强制完全一致的计算路径 #

• 既然已经用pyfma让CPU和GPU对齐，就得保证所有运算的顺序丝毫不差：比如数组累加不要依赖广播的自动优化，改成显式循环；GPU端避免让编译器自动重排运算顺序，手动固定每一步的计算流程。
• 乘加操作全用FMA：CPU用pyfma.fma(a,b,c)，GPU直接调用CUDA原生的fma(a,b,c)，确保两边的乘加运算完全等价。

3. 用高精度中间变量压减累积误差 #

• 别用单精度float32，直接换成双精度float64：Numba和CUDA都完美支持，现代GPU对double的性能已经足够用，能大幅降低单精度带来的误差。
• 累加操作改用Kahan求和算法，不管CPU还是GPU都实现一模一样的逻辑：

  @numba.njit(fastmath=False)
  def kahan_sum(arr):
      total = 0.0
      error = 0.0
      for x in arr:
          y = x - error
          t = total + y
          error = (t - total) - y
          total = t
      return total
  

  这个算法能有效减少累加时的舍入误差，让结果更接近高精度参考。

4. 对齐高精度参考的计算逻辑 #

• 用mpmath算参考结果时，别直接用它的内置函数，得完全复刻CPU/GPU的计算步骤：比如你的CPU代码是循环累加+FMA，那mpmath版本也要写一样的循环，用mpmath的FMA实现，这样得到的参考结果才是和你的计算路径对应的高精度值，而不是数学上的理想值——毕竟理想值和实际计算的误差不是一回事。

5. 定位误差源并修正 #

• 算出CPU/GPU结果和高精度参考的差值，找到误差最大的环节：是除法、开方这类单步运算的问题，还是累积累加的锅。
• 针对误差大的运算做优化：比如开方用牛顿迭代多跑一轮提升精度；把硬编码的常数（比如π）换成mpmath计算的高精度版本，别用默认的math.pi。

6. 验证结果 #

• 用numpy.testing.assert_allclose做验证，设置rtol=1e-5（对应5位小数）、atol=1e-6，确保三者的差异在允许范围内。
• GPU结果要先拷贝到CPU再对比，避免设备传输的额外误差。

内容的提问来源于stack exchange，提问作者slaw