C语言计算数学表达式时不同浮点数类型结果差异问题咨询

• 核心原因：灾难性抵消（Catastrophic Cancellation）
  从代数角度推导，原式可化简为：

  [((a+b)^2 - (a^2 + 2ab)) / b^2 = b^2 / b^2 = 1]
  但浮点数存储和计算存在固定精度限制：float仅支持约6_{7位十进制有效数字，double支持约15}17位十进制有效数字。
  计算a + b时，a为1000（1e3）、b为0.0001（1e-4），和为1000.0001，需要8位有效数字才能精确存储：

  • float精度不足，存储a + b时就已经产生了较大的截断误差
  • double精度足够，a + b的存储误差极小
    后续执行(a+b)² - (a² + 2ab)计算时，两个数值极为接近的大数相减会抵消全部高位有效数字，仅剩低位的误差部分，也就是灾难性抵消现象：
  • double场景下残余误差很小，除以b²（1e-8）放大后仅产生约千分之一的偏差，结果接近理论值1
  • float场景下a + b的初始误差就很大，抵消后残余误差远大于真实的b²值，再被1e8倍放大后就出现了6250000的异常结果
• 额外误差来源：pow/powf函数的精度损耗
  通用幂函数pow/powf是针对任意指数场景设计的，计算整数平方的精度远低于直接用乘法(a+b)*(a+b)，额外引入了不必要的计算误差。如果把代码中的pow调用全部替换为直接乘法，计算精度会有明显提升。
• 优化建议
  实际开发中遇到这类数值计算场景，要优先使用代数化简后的公式，避免出现两个接近值相减的计算逻辑，就能从根源上规避这类精度问题。

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