依赖于浮点相等性的简单集成

Question

我有以下非常粗略的积分计算器：

// definite integrate on one variable
// using basic trapezoid approach
float integrate(float start, float end, float step, float (*func)(float x))
{
    if (start >= (end-step))
        return 0;
    else {
        float x = start; // make it a bit more math-like
        float segment = step * (func(x) + func(x+step))/2;
        return segment + integrate(x+step, end, step, func);
    }
}

以及用法示例：

static float square(float x) {return x*x;}
int main(void)
{
    // Integral x^2 from 0->2 should be ~ 2.6
    float start=0.0, end=2.0, step=0.01;
    float answer = integrate(start, end, step, square);
    printf("The integral from %.2f to %.2f for X^2 = %.2f\n", start, end, answer );
}

$ run
The integral from 0.00 to 2.00 for X^2 = 2.67

如果 start >= (end-step) 处的相等性检查不起作用会怎样？例如，如果它的计算结果为 2.99997 而不是 3，那么另一个循环（或少一个循环）也是如此。有没有办法防止这种情况发生，或者大多数数学类型的计算器是否只使用小数或 'normal' 浮点数的扩展？

Answer 1

如果给定 step，一种编写循环的方法（您应该为此使用循环，而不是递归）是：

float x;
for (float i = 0; (x = start + i*step) < end - step/2; ++i)
    …

关于此的一些要点：

我们用 i 保持整数计数。只要步数合理，这里面就不会出现浮点数舍入错误。（我们可以使 i 和 int，但是 float 可以很好地计算整数值，并且使用 float 可以避免 int-to-float 转换 i*step.)
我们不是重复递增 x（或递归传递的 start），而是每次都将其重新计算为 start + i*step。这只有两个可能的舍入误差，在乘法和加法中，因此它避免了通过重复加法累积误差。
我们使用end - step/2作为阈值。即使计算出的 x 偏离 end 与 end - step/2 一样远，这也使我们能够捕捉到所需的端点。这是我们能做的最好的，因为如果它漂移距离理想间隔点超过半个 step，我们无法判断它是否从 end-step 漂移 +step/2 或-step/2 来自 end.

这假定 step 是 end-start 的整数除法，或者非常接近它，因此循环中有整数步。如果不是，循环应该重新设计一点，提前停止一步，然后在最后计算部分宽度的一步。

开头提到了被给step。另一种方法是，您可能会得到一些要使用的步长，然后从中计算步长。在那种情况下，我们将使用整数步数来控制循环。循环终止条件根本不涉及浮点舍入。我们可以将 x 计算为 (float) i / NumberOfSteps * (end-start) + start.

Answer 2

可以轻松进行两项改进。

使用递归是个坏主意。每个额外的调用都会创建一个新的堆栈框架。对于足够多的步骤，您将触发堆栈溢出。请改用循环。
通常，您可以使用 start、end 和 n 步数来避免舍入问题。第 k 个间隔的位置将在 start + k * (end - start) / n;

因此您可以将函数重写为

float integrate(float start, float end, int n, float (*func)(float x))
{
    float next = start;
    float sum = 0.0f;
    for(int k = 0; k < n; k++) {
        float x = next;
        next = start + k * (end - start) / n; 
        sum += 0.5f * (next - x) * (func(x) + func(next));
    }
    return sum;
}

依赖于浮点相等性的简单集成

Simple integration that depends on floating point equality

c

math

floating-point

decimal