C中浮点运算中的下溢错误

Question

我是 C 新手，我的任务是创建一个函数

f(x) = sqrt[(x^2)+1]-1

它可以处理非常大的数字和非常小的数字。我正在在线界面上提交我的脚本以检查我的答案。

对于非常大的数字，我将表达式简化为：

f(x) = x-1

只需使用最高功率。这是正确的答案。

同样的逻辑不适用于较小的数字。对于小数（大约 1e-7），它们会很快被截断为零，甚至在平方之前。我怀疑这与 C 中的浮点精度有关。在我的教科书中，它说浮点类型的最小可能值为 1.17549e-38，精度为 6 位。所以虽然 1e-7 比 1.17e-38 大很多，但它的精度更高，因此四舍五入为零。这是我的猜测，如果我错了，请纠正我。

作为一种解决方案，我认为当 x

#include <math.h>
#include <stdio.h>

double feval(double x) {
    /* Insert your code here */
    if (x > 1e299) 
    {;
        return x-1;
    }
    if (x < 1e-6)
    {
        long double g;
        g = x;
        printf("x = %Lf\n", g);
        long double a;
        a = pow(x,2);
        printf("x squared = %Lf\n", a);
        return sqrt(g*g+1.)- 1.;
    }
    else
    { 
        printf("x = %f\n", x);
        printf("Used third \n");
        return sqrt(pow(x,2)+1.)-1;
    }
}

int main(void)
{
    double x;
    printf("Input: ");
    scanf("%lf", &x);
    double b;
    b = feval(x);
    printf("%f\n", b);
    return 0;
}

Answer 1

在这些情况下通常有用的一个技巧是基于身份

(a+1)*(a-1) = a*a-1

在这种情况下

sqrt(x*x+1)-1 = (sqrt(x*x+1)-1)*(sqrt(x*x+1)+1) 
                 /(sqrt(x*x+1)+1)
= (x*x+1-1) / (sqrt(x*x+1)+1)
= x*x/(sqrt(x*x+1)+1)

最后一个公式可以作为实现。对于 vwry small x sqrt(x*x+1)+1 将接近 2（对于足够小的 x 它将是 2）但我们不会放松评估它的精度。

Answer 2

以朴素的方式实现这里有两个问题：计算x * x时中间计算中的溢出或下溢，以及最终减1期间的减法取消。第二个问题是准确性问题。

ISO C 有一个标准数学函数hypot (x, y)，它可以准确地执行sqrt (x * x + y * y) 的计算，同时避免中间计算中的下溢和上溢。解决减法消除问题的一种常见方法是对计算进行代数转换，以便将其转换为乘法和/或除法。

结合这两个修复导致float 参数的以下实现。根据我的测试，它在所有可能的输入中的错误小于 3 ulps。

/* Compute sqrt(x*x+1)-1 accurately and without spurious overflow or underflow */
float func (float x)
{
    return (x / (1.0f + hypotf (x, 1.0f))) * x;
}

Answer 3

问题不在于跑到最小值，而在于精度。

正如您自己所说，您机器上的float 的精度约为 7 位数。所以让我们取 x = 1e-7，所以 x^2 = 1e-14。这仍然在float 的范围内，没有问题。但现在加 1。确切的答案是1.00000000000001。但如果我们只有 7 位精度，则将四舍五入为 1.0000000，即正好 1。所以你最终计算出的 sqrt(1.0)-1 正好是 0。

一种方法是使用sqrt 在x=1 和sqrt(x) ~ 1+0.5*(x-1) 周围的线性近似。这将导致近似f(x) ~ 0.5*x^2。

Answer 4

对于较小的输入，执行 1+x^2 时会出现截断错误。如果x=1e-7f，x*x 将很高兴地适合 32 位浮点数（由于1e-7 没有精确的浮点表示，但x*x 会很多小于 1 表示浮点精度不足以表示 1+x*x。

对 sqrt(1+x^2) 进行泰勒展开会更合适，最低阶为

sqrt(1+x^2) = 1 + 0.5*x^2 + O(x^4)

然后，你可以把你的结果写成

sqrt(1+x^2)-1 = 0.5*x^2 + O(x^4),

避免将非常小的数字添加到 1 的情况。

附带说明，您不应将pow 用于整数幂。对于 x^2，你应该只做x*x。任意整数幂的效率要高一些。例如，GNU 科学库具有高效计算任意整数幂的功能。