复合辛普森法则取大n值时积分计算结果偏大问题咨询

我使用复合辛普森法则计算积分∫₂¹⁰⁰⁰(1/ln(x))dx，当n值很大（约500000）时，计算结果开始与参考值176.5644出现偏差，比如n=10000000时结果为184.1495。按理论n越大精度应越高，想请教该现象的原因。相关C++代码如下：

#include <iostream>
#include <cmath>

// 被积函数f(x)
float f(float x)
{
    return (float) 1 / std::log(x);
}

float my_simpson(float a, float b, long int n)
{
  if (n % 2 == 1) n += 1; // 复合辛普森法则要求n为偶数
  float area, h = (b-a)/n;

  float x, y, z;
  for (int i = 1; i <= n/2; i++)
  {
    x = a + (2*i - 2)*h;
    y = a + (2*i - 1)*h;
    z = a + 2*i*h;
    area += f(x) + 4*f(y) + f(z);
  }

  return area*h/3;
}

int main()
{
    std::cout.precision(20);
    int upperBound = 1'000;
    int subsplits = 1'000'000;

    float approx = my_simpson(2, upperBound, subsplits);

    std::cout << "Output: " << approx << std::endl;

    return 0;
}

原因分析 #

问题核心是单精度浮点数（float）的精度限制：

• float类型仅能提供约7位十进制有效数字，当n极大时，步长h=(b-a)/n会变得非常小（比如n=1e7时，h≈9.98e-5）。
• 循环中需要累加数百万次微小的计算项，float的舍入误差会随着累加次数的增加不断累积，最终完全掩盖了正确的计算结果，导致输出值偏离参考值。

解决验证 #

更新：将代码中所有float类型替换为double类型（双精度浮点数，提供约15-17位十进制有效数字）后，计算结果恢复正常，与参考值一致。

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