双精度奇怪的行为.需要解释

Question

这是代码:

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

static double const x = 665857;
static double const y = 470832;

int main(){
    double z = x*x*x*x -(y*y*y*y*4+y*y*4);
    printf("%f \n",z);
    return 0;
}

神秘地(对我来说)这个代码打印"0.0"如果在32位机器上编译(或者在64位机器上使用-m32标志,就像我的情况一样)使用GCC 4.6.据我所知,浮点运算可能会溢出/下溢它们或者它们会失去精度,但是...... 0？怎么样？

提前致谢.

Answer 1

问题不在于数字溢出.问题是双精度没有足够的精度来区分减法的两个操作数.

值为x*x*x*x196573006004558194713601.

值为y*y*y*y*4+y*y*4196573006004558194713600.

这些数字有78位,只有最后一位不同.双精度数字只有53位.其他数字四舍五入到只有53位.

在您的情况下,两个操作数四舍五入到相同的数字,因此它们的差异为0.

如果你稍微重写z的表达式,就会发生更奇怪的事情:

double z = x * x * x * x - ((y * y + 1) * y * y * 4);

通过此更改,您将获得33554432!为什么？因为舍入中间结果的方式导致右操作数的最后一位不同.最后一位的值是2 ^(78-53)= 2 ^ 25.

Answer 2

使用任意精度整数计算表达式:

Prelude> 665857^4 - 4*(470832^4 + 470832^2)
1

由于double通常只有53位精度而中间结果有78位,因此精度不足以精确计算结果,因此它是四舍五入的,最后一位在某些时候被遗忘.