c - 如何在 C double 中表示无穷大？

Question

我从计算机系统:程序员的视角一书中了解到，IEEE 标准要求使用以下 64 位二进制格式表示 double float :

• s: 1位符号
• exp: 11位指数
• frac:52 位的分数

+infinity 表示为具有以下模式的特殊值:

• s = 0
• 所有的 exp 位都是 1
• 小数位全部为0

而且我认为 double 的完整 64 位应该按以下顺序排列:

(s)(exp)(分数)

所以我写了下面的C代码来验证一下:

//Check the infinity
double x1 = (double)0x7ff0000000000000;  // This should be the +infinity
double x2 = (double)0x7ff0000000000001; //  Note the extra ending 1, x2 should be NaN
printf("\nx1 = %f, x2 = %f sizeof(double) = %d", x1,x2, sizeof(x2));
if (x1 == x2)
    printf("\nx1 == x2");
else
    printf("\nx1 != x2");

但结果是:

x1 = 9218868437227405300.000000, x2 = 9218868437227405300.000000 sizeof(double) = 8
x1 == x2

为什么这个数字是有效数字而不是一些无穷大错误？

为什么 x1==x2？

(我使用的是 MinGW GCC 编译器。)

加1

我修改了如下代码并成功验证了 Infinity 和 NaN。

//Check the infinity and NaN
unsigned long long x1 = 0x7ff0000000000000ULL; // +infinity as double
unsigned long long x2 = 0xfff0000000000000ULL; // -infinity as double
unsigned long long x3 = 0x7ff0000000000001ULL; // NaN as double
double y1 =* ((double *)(&x1));
double y2 =* ((double *)(&x2));
double y3 =* ((double *)(&x3));

printf("\nsizeof(long long) = %d", sizeof(x1));
printf("\nx1 = %f, x2 = %f, x3 = %f", x1, x2, x3); // %f is good enough for output
printf("\ny1 = %f, y2 = %f, y3 = %f", y1, y2, y3);

结果是:

sizeof(long long) = 8
x1 = 1.#INF00, x2 = -1.#INF00, x3 = 1.#SNAN0
y1 = 1.#INF00, y2 = -1.#INF00, y3 = 1.#QNAN0

详细输出看起来有点奇怪，但我认为重点很明确。

PS.: 指针转换好像没必要。只需使用 %f 告诉 printf 函数以 double 格式解释 unsigned long long 变量。

添加 2

出于好奇，我用下面的代码检查了变量的位表示。

typedef unsigned char *byte_pointer;

void show_bytes(byte_pointer start, int len)
{
    int i;
    for (i = len-1; i>=0; i--)
    {
        printf("%.2x", start[i]);
    }
    printf("\n");
}

我尝试了下面的代码:

//check the infinity and NaN
unsigned long long x1 = 0x7ff0000000000000ULL; // +infinity as double
unsigned long long x2 = 0xfff0000000000000ULL; // -infinity as double
unsigned long long x3 = 0x7ff0000000000001ULL; // NaN as double
double y1 =* ((double *)(&x1));
double y2 =* ((double *)(&x2));
double y3 = *((double *)(&x3));

unsigned long long x4 = x1 + x2;  // I want to check (+infinity)+(-infinity)
double y4 = y1 + y2; // I want to check (+infinity)+(-infinity)

printf("\nx1: ");
show_bytes((byte_pointer)&x1, sizeof(x1));
printf("\nx2: ");
show_bytes((byte_pointer)&x2, sizeof(x2));
printf("\nx3: ");
show_bytes((byte_pointer)&x3, sizeof(x3));
printf("\nx4: ");
show_bytes((byte_pointer)&x4, sizeof(x4));

printf("\ny1: ");
show_bytes((byte_pointer)&y1, sizeof(y1));
printf("\ny2: ");
show_bytes((byte_pointer)&y2, sizeof(y2));
printf("\ny3: ");
show_bytes((byte_pointer)&y3, sizeof(y3));
printf("\ny4: ");
show_bytes((byte_pointer)&y4, sizeof(y4));

输出是:

x1: 7ff0000000000000

x2: fff0000000000000

x3: 7ff0000000000001

x4: 7fe0000000000000

y1: 7ff0000000000000

y2: fff0000000000000

y3: 7ff8000000000001

y4: fff8000000000000  // <== Different with x4

奇怪的是，虽然 x1 和 x2 与 y1 和 y2 具有相同的位模式，但 x4 的和与 y4 不同。

和

printf("\ny4=%f", y4);

给出这个:

y4=-1.#IND00  // What does it mean???

为什么不同？ y4是怎么得到的？

Answer 1

首先，0x7ff0000000000000 确实是双无穷大的位表示。但是强制转换没有设置位表示，它把 0x7ff0000000000000 的逻辑值解释为 64 位整数。因此，您需要使用一些其他方式来设置位模式。

设置位模式的直接方法是

uint64_t bits = 0x7ff0000000000000;
double infinity = *(double*)&bits;

但是，这是未定义的行为。C 标准禁止读取已存储为一种基本类型 (uint64_t) 的值作为另一种基本类型 (双)。这被称为严格的别名规则，允许编译器生成更好的代码，因为它可以假设一种类型的读取顺序和另一种类型的写入顺序是不相关的。

此规则的唯一异常(exception)是 char 类型:明确允许您将任何指针转换为 char* 并返回。所以你可以尝试使用这段代码:

char bits[] = {0x7f, 0xf0, 0, 0, 0, 0, 0, 0};
double infinity = *(double*)bits;

即使这不再是未定义的行为，它仍然是实现定义的行为:double 中字节的顺序取决于您的机器。给定的代码适用于像 ARM 和 Power 系列这样的大端机器，但不适用于 X86。对于 X86，您需要此版本:

char bits[] = {0, 0, 0, 0, 0, 0, 0xf0, 0x7f};
double infinity = *(double*)bits;

确实没有办法绕过这个实现定义的行为，因为不能保证机器会按照与整数值相同的顺序存储浮点值。甚至有机器使用这样的字节顺序:<1, 0, 3, 2> 我什至不想知道是谁想出了这个绝妙的主意，但它存在，我们必须接受它。

---

关于你的最后一个问题:浮点运算与整数运算有着本质上的不同。这些位具有特殊含义，浮点单元会考虑到这一点。特别是无穷大、NAN 和非规范化数字等特殊值以特殊方式处理。由于 +inf + -inf 被定义为生成 NAN，因此您的浮点单元会发出 NAN 的位模式。整数单元不知道无穷大或 NAN，因此它只是将位模式解释为一个巨大的整数并愉快地执行整数加法(在这种情况下恰好溢出)。生成的位模式不是 NAN 的位模式。它恰好是一个非常大的正 float (2^1023，准确地说)的位模式，但这没有任何意义。

---

实际上，有一种方法可以以可移植的方式设置除 NAN 之外的所有值的位模式:给定包含符号位、指数位和尾数位的三个变量，您可以这样做:

uint64_t sign = ..., exponent = ..., mantissa = ...;
double result;
assert(!(exponent == 0x7ff && mantissa));    //Can't set the bits of a NAN in this way.
if(exponent) {
    //This code does not work for denormalized numbers. And it won't honor the value of mantissa when the exponent signals NAN or infinity.
    result = mantissa + (1ull << 52);    //Add the implicit bit.
    result /= (1ull << 52);    //This makes sure that the exponent is logically zero (equals the bias), so that the next operation will work as expected.
    result *= pow(2, (double)((signed)exponent - 0x3ff));    //This sets the exponent.
} else {
    //This code works for denormalized numbers.
    result = mantissa;    //No implicit bit.
    result /= (1ull << 51);    //This ensures that the next operation works as expected.
    result *= pow(2, -0x3ff);    //Scale down to the denormalized range.
}
result *= (sign ? -1.0 : 1.0);    //This sets the sign.

这使用浮点单元本身将位移动到正确的位置。由于无法使用浮点运算与 NAN 的尾数位进行交互，因此无法在此代码中包含 NAN 的生成。好吧，您可以生成一个 NAN，但您无法控制其尾数位模式。