依赖于浮点相等的简单积分

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，那么另一个循环（或少一个循环）也是如此。有没有办法防止这种情况发生，或者大多数数学类型的计算器只能使用小数或“正常”浮点数的扩展？

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 还是从@987654342 漂移了-step/2 @。

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

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

Answer 2

可以轻松进行两项改进。

1. 使用递归是个坏主意。每个额外的调用都会创建一个新的堆栈帧。对于足够多的步骤，您将触发堆栈溢出。请改用循环。
2. 通常，您可以通过使用步数start、end 和n 来避免舍入问题。 kth 间隔的位置将在 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;
}