floating-point - 是否存在 x 的浮点值，其中 xx == 0 为假？

Question

在大多数情况下，我知道应该使用一系列值（abs(xy) < epsilon）来实现浮点比较测试，但是自减法是否意味着结果为零？

// can the assertion be triggered?
float x = //?;
assert( x-x == 0 )

我的猜测是 nan/inf 可能是特殊情况，但我更感兴趣的是简单值会发生什么。

编辑：

如果有人可以引用参考（IEEE 浮点标准），我很乐意选择答案？

Answer 1

正如您所暗示的那样，inf - infis NaN，它不等于零。同样，NaN - NaN是NaN。然而，对于任何有限浮点数来说，确实如此x（x - x == 0.0取决于舍入模式，结果可能是负零，但在浮点算术中x - x，负零比较等于）。0.0

编辑：给出明确的标准参考有点棘手，因为这是 IEEE-754 标准中规定的规则的一个新兴属性。具体来说，它遵循第 5 条中定义的运算正确舍入的要求。减法就是这样的运算（第 5.4.1 节“算术运算”），正确舍入的结果x - x是适当符号的零（第 6.3 节，第 3 段）：

当两个符号相反的操作数之和（或两个符号相同的操作数之差）恰好为零时，该和（或差）的符号在除roundTowardNegative之外的所有舍入方向属性中应为+0；在该属性下，精确零和（或差）的符号应为 -0。

所以x - xmust 的结果是+/- 0，因此必须比较等于0.0（第 5.11 节，第 2 段）：

比较应忽略零符号。

进一步编辑：这并不是说有缺陷的编译器不会导致该断言触发。您的问题模棱两可；没有有限的浮点数x是x - x == 0错误的。但是，这不是您发布的代码检查的内容；它检查 C 风格语言中的某个表达式是否可以计算为非零值；特别是，在某些平台上，通过某些（考虑不周的）编译器优化，该x表达式中变量的两个实例可能具有不同的值，导致断言失败（特别是如果x是某些计算的结果，而不是常量，可表示的值）。这是这些平台上数字模型中的一个错误，但这并不意味着它不会发生。

Answer 2

如果表示被转换（例如从 64 位内存格式到 x86 上的 80 位内部寄存器格式），我希望断言在某些情况下可能会触发。

Answer 3

是的，除了特殊情况x-x将始终为 0。但x*(1/x)不会始终为 1 ;-)

Answer 4

是的，除特殊情况外，自减应始终导致零。

在调整指数和尾数的比较之前加、减、乘或除会出现问题。当指数相同时，减去尾数，如果它们相同，则一切都为零。

http://grouper.ieee.org/groups/754/

Answer 5

我对主要问题的回答：“是否存在 x 的浮点值，其中 xx == 0 为假？” 是：至少英特尔处理器上的浮点实现不会在“+”和“-”操作中产生算术下溢，因此您将无法找到 xx == 0 为假的 x。对于所有支持 IEEE 754-2008 的处理器也是如此（请参阅下面的参考资料）。

我对另一个你的问题的简短回答：如果 (xy == 0) 和 (x == y) 一样安全，所以 assert(xx == 0) 是可以的，因为在xx 或 ( xy)。

原因如下。浮点数/双精度数将以尾数和二进制指数的形式保存在内存中。在标准情况下，尾数被归一化：它 >= 0.5 并且 < 1。在<float.h>您可以从 IEEE 浮点标准中找到一些常量。现在对我们来说有趣的是只关注

#define DBL_MIN         2.2250738585072014e-308 /* min positive value */
#define DBL_MIN_10_EXP  (-307)                  /* min decimal exponent */
#define DBL_MIN_EXP     (-1021)                 /* min binary exponent */

但不是每个人都知道，你可以有比DBL_MIN少的双倍数字。如果您对 DBL_MIN 下的数字进行算术运算，则该数字将不会被规范化，因此您可以像使用整数一样使用该数字（仅使用尾数进行操作），而不会出现任何“舍入错误”。

备注：我个人尽量不要使用“round errors”这个词，因为在算术计算机操作中没有错误。这些操作仅与具有相同计算机数（如浮点数）的 +、-、* 和 / 操作不同。浮点数子集上有确定性操作，可以以（尾数，指数）的形式保存，每个浮点数都有明确定义的位数。我们可以将这种浮点子集命名为计算机浮点数。所以经典浮点运算的结果将被投影回计算机浮点数集。这种投影操作是确定性的，并且具有很多特征，例如 if x1 >= x2 then x1*y >= x2*y。

对不起，长话，回到我们的主题。

为了准确显示如果我们使用小于 DBL_MIN 的数字进行操作，我用 C 语言编写了一个小程序：

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

void DumpDouble(double d)
{
    unsigned char *b = (unsigned char *)&d;
    int i;

    for (i=1; i<=sizeof(d); i++) {
        printf ("%02X", b[sizeof(d)-i]);
    }
    printf ("\n");
}

int main()
{
    double x, m, y, z;
    int exp;

    printf ("DBL_MAX=%.16e\n", DBL_MAX);
    printf ("DBL_MAX in binary form: ");
    DumpDouble(DBL_MAX);

    printf ("DBL_MIN=%.16e\n", DBL_MIN);
    printf ("DBL_MIN in binary form: ");
    DumpDouble(DBL_MIN);

    // Breaks the floating point number x into its binary significand
    // (a floating point value between 0.5(included) and 1.0(excluded))
    // and an integral exponent for 2
    x = DBL_MIN;
    m = frexp (x, &exp);
    printf ("DBL_MIN has mantissa=%.16e and exponent=%d\n", m, exp);
    printf ("mantissa of DBL_MIN in binary form: ");
    DumpDouble(m);

    // ldexp() returns the resulting floating point value from
    // multiplying x (the significand) by 2
    // raised to the power of exp (the exponent).
    x = ldexp (0.5, DBL_MIN_EXP);   // -1021
    printf ("the number (x) constructed from mantissa 0.5 and exponent=DBL_MIN_EXP (%d) in binary form: ", DBL_MIN_EXP);
    DumpDouble(x);

    y = ldexp (0.5000000000000001, DBL_MIN_EXP);
    m = frexp (y, &exp);
    printf ("the number (y) constructed from mantissa 0.5000000000000001 and exponent=DBL_MIN_EXP (%d) in binary form: ", DBL_MIN_EXP);
    DumpDouble(y);
    printf ("mantissa of this number saved as double will be displayed by printf(%%.16e) as %.16e and exponent=%d\n", m, exp);

    y = ldexp ((1 + DBL_EPSILON)/2, DBL_MIN_EXP);
    m = frexp (y, &exp);
    printf ("the number (y) constructed from mantissa (1+DBL_EPSILON)/2 and exponent=DBL_MIN_EXP (%d) in binary form: ", DBL_MIN_EXP);
    DumpDouble(y);
    printf ("mantissa of this number saved as double will be displayed by printf(%%.16e) as %.16e and exponent=%d\n", m, exp);

    z = y - x;
    m = frexp (z, &exp);
    printf ("z=y-x in binary form: ");
    DumpDouble(z);
    printf ("z will be displayed by printf(%%.16e) as %.16e\n", z);
    printf ("z has mantissa=%.16e and exponent=%d\n", m, exp);

    if (x == y)
        printf ("\"if (x == y)\" say x == y\n");
    else
        printf ("\"if (x == y)\" say x != y\n");

    if ((x-y) == 0)
        printf ("\"if ((x-y) == 0)\" say \"(x-y) == 0\"\n");
    else
        printf ("\"if ((x-y) == 0)\" say \"(x-y) != 0\"\n");
}

此代码产生以下输出：

DBL_MAX=1.7976931348623157e+308
DBL_MAX in binary form: 7FEFFFFFFFFFFFFF
DBL_MIN=2.2250738585072014e-308
DBL_MIN in binary form: 0010000000000000
DBL_MIN has mantissa=5.0000000000000000e-001 and exponent=-1021
mantissa of DBL_MIN in binary form: 3FE0000000000000
the number (x) constructed from mantissa 0.5 and exponent=DBL_MIN_EXP (-1021) in binary form: 0010000000000000
the number (y) constructed from mantissa 0.5000000000000001 and exponent=DBL_MIN_EXP (-1021) in binary form: 0010000000000001
mantissa of this number saved as double will be displayed by printf(%.16e) as 5.0000000000000011e-001 and exponent=-1021
the number (y) constructed from mantissa (1+DBL_EPSILON)/2 and exponent=DBL_MIN_EXP (-1021) in binary form: 0010000000000001
mantissa of this number saved as double will be displayed by printf(%.16e) as 5.0000000000000011e-001 and exponent=-1021
z=y-x in binary form: 0000000000000001
z will be displayed by printf(%.16e) as 4.9406564584124654e-324
z has mantissa=5.0000000000000000e-001 and exponent=-1073
"if (x == y)" say x != y
"if ((x-y) == 0)" say "(x-y) != 0"

所以我们可以看到，如果我们使用小于 DBL_MIN 的数字，它们将不会被规范化（参见 参考资料0000000000000001）。我们正在使用这些数字，就像使用整数一样，并且没有任何“错误”。因此，如果我们指定y=xthenif (x-y == 0)就和 一样安全if (x == y)，并且assert(x-x == 0)可以正常工作。在此示例中，z = 0.5 * 2 ^(-1073) = 1 * 2 ^(-1072)。这个数字实际上是我们可以双倍保存的最小数字。所有数字小于 DBL_MIN 的算术运算都像整数乘以 2 ^(-1072) 一样工作。

因此，我在配备 Intel 处理器的 Windows 7 计算机上没有下溢问题。如果有人有另一个处理器，比较我们的结果会很有趣。

有人知道如何使用 - 或 + 操作产生算术下溢吗？我的实验看起来像这样，那是不可能的。

已编辑：我稍微修改了代码以提高代码和消息的可读性。

附加链接：我的实验表明，http ://grouper.ieee.org/groups/754/faq.html#underflow 在我的 Intel Core 2 CPU 上是绝对正确的。它的计算方式不会在“+”和“-”浮点运算中产生下溢。我的结果独立于 Strict (/fp:strict) 或 Precise (/fp:precise) Microsoft Visual C 编译器开关（请参阅http://msdn.microsoft.com/en-us/library/e7s85ffb%28VS.80%29 .aspx和http://msdn.microsoft.com/en-us/library/Aa289157）

另一个（可能是最后一个）链接和我的最后一句话：我找到了一个很好的参考http://en.wikipedia.org/wiki/Subnormal_numbers，其中的描述与我之前写的相同。包括非正规数或非正规数（现在通常称为次正规数，例如在IEEE 754-2008 中）遵循以下声明：

“非正规数保证浮点数的加法和减法永远不会下溢；两个相邻的浮点数总是有一个可表示的非零差。如果没有逐渐下溢，即使值不相等，减法 a−b 也会下溢并产生零。”</p>

所以我所有的结果在任何支持 IEEE 754-2008 的处理器上都必须是正确的。

Answer 6

关于马克所说的——查看这个链接http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.18。（但不确定它是否适用于您的情况。）