logo标签列表

如何检查浮点数的依赖关系

c++algorithmfloating-pointfloating-accuracyinverse
12

我想在C++中确定一个浮点数是否是另一个浮点数的乘法逆元。问题是我必须使用第三个变量来完成。例如，这段代码：

float x=5,y=0.2;
if(x==(1/y)) cout<<"They are the multiplicative inverse of eachother"<<endl;
else cout<<"They are NOT the multiplicative inverse of eachother"<<endl;

将输出“they are not...”，这是错误的，而此代码：

float x=5,y=0.2,z;
z=1/y;
if(x==z) cout<<"They are the multiplicative inverse of eachother"<<endl;
else cout<<"They are NOT the multiplicative inverse of eachother"<<endl;

将输出："they are..."，这是正确的。
为什么会这样发生？

7
好的，我会尽可能地翻译这篇 Oracle 文章并使其更易懂。以下是您需要翻译的内容：http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html - Mysticial
3
((x*y) == 1)也不起作用吗？ - Vyktor
我在我的回答中添加了一些信息。再加上一分，感谢这个有成效的问题。 - Gangnus
5个回答
36

浮点数精度问题

    你有两个问题，但这两个问题都来自同一个根源。

浮点数无法精确比较，无法精确相减或相除，也无法对任何值进行精确计数。任何涉及浮点数的运算都可能（几乎总是）在结果中引入一些误差。即使 a=0.2f 也不是精确的操作。这些深层原因在其他答案的作者们已经很好地解释了。 （谢谢他们的回答和投票。）

这里出现了你的第一个，也是更简单的错误。永远不要、绝不要、再也不要 在它们上使用 == 或其等价物。

而应该使用 Abs(a-b)<HighestPossibleError 代替 a==b。

    但这并不是你任务中唯一的问题。

Abs(1/y-x)<HighestPossibleError 也不行。至少，它不够可靠。为什么？

假设 x=1000，y=0.001，并且 y 的起始相对误差为 10-6。

（相对误差=误差/值）。

值的相对误差在乘除法中是相加的。

1/y 约等于 1000。它的相对误差也是 10-6。（“1”没有误差）

这使得绝对误差=1000*10-6=0.001。当你稍后减去 x 时，该误差将是唯一剩下的。（在加减法中，绝对误差会相加，而 x 的误差则微不足道。）肯定，你不能指望如此大的误差，HighestPossibleError 一定设置得更低，否则你的程序将抛弃一个好的 x、y 组合。

因此，浮点数运算的下两个规则是：尽量不要除以较小的值，并且避免在此之后相减接近的值。

有两种简单方法可以解决这个问题。

  • 找出 x、y 中绝对值较大的一个，然后将 1 除以较大的那个数，最后再减去较小的数。

  • 如果你想比较 1/y 和 x，在你还在使用字母进行操作而非值，并且你的运算没有误差的情况下，将比较的两边都乘以 y，就可以得到 1 against x*y。（通常你应该检查这个运算中的符号，但在这里我们使用了绝对值，所以是可行的。） 这样比较的结果根本没有除法。

简单一点说：

1/y V x   <=>   y*(1/y) V x*y   <=>   1 V x*y

我们已经知道应该这样进行1对x*y的比较:

const float HighestPossibleError=1e-10;
if(Abs(x*y-1.0)<HighestPossibleError){...

就是这样。

顺便提一句，如果你真的需要把所有东西都放在一行里，可以使用：

if(Abs(x*y-1.0)<1e-10){...

但这是不好的风格。我不建议这样做。

P.P.S. 在您的第二个示例中，编译器会优化代码，使其在运行任何代码之前将z设置为5。因此，将5与5进行比较即使对于浮点数也可以工作。

13
问题在于 0.2 不能用二进制精确表示，因为它的二进制展开有无限多个数字:
 1/5: 0.0011001100110011001100110011001100110011...

这就像是在十进制中无法精确表示1/3一样。因为x是存储在一个具有有限比特数的float中，所以这些数字在某个时刻会被截断，例如:

   x: 0.0011001100110011001100110011001
问题的产生是因为CPU通常在内部使用更高的精度，因此当你刚刚计算出1/y时，结果将具有更多的数字，在加载x进行比较时，x会被扩展以匹配CPU的内部精度。
 1/y: 0.0011001100110011001100110011001100110011001100110011
   x: 0.0011001100110011001100110011001000000000000000000000

因此，当您进行逐位比较时，它们是不同的。

然而，在您的第二个示例中，将结果存储到一个变量中意味着在进行比较之前它被截断，因此在这个精度下比较它们是相等的：

   x: 0.0011001100110011001100110011001
   z: 0.0011001100110011001100110011001

许多编译器都有开关可以启用，以强制中间值在每一步被截断以保证精度的一致性，但通常建议避免直接比较浮点数，而是检查它们是否相差不到某个 epsilon 值，这就是Gangnus建议的内容。

5
你需要精确地定义两个近似值成为乘法逆的含义，否则你将不知道你应该测试什么。
0.2在二进制中没有精确的表示方式。如果你储存那些没有精确表示的数字，只使用有限的精度，你得到的答案将不是完全正确的。
在十进制中同样如此。例如，1/3没有精确的十进制表示方式。你可以将它储存为0.333333。但是这时会产生一个问题。3和0.333333是否为乘法逆？如果你将它们相乘，你得到0.999999。如果你希望答案是“是”，你必须创建一个测试乘法逆的方法，它不像简单的乘法和等于1一样简单。
在二进制中同样如此。
2
其他回复中的讨论很好，所以我不会重复任何内容，但没有代码。下面是一小段代码，用于检查一对浮点数相乘是否恰好等于1.0。
该代码做出了一些假设/断言（在x86平台上通常满足）：
- float是32位二进制（也称为单精度）IEEE-754
- int或long均为32位（我决定不依赖uint32_t的可用性）
- memcpy()将float复制到int/long，使得8873283.0f变成0x4B076543（即期望某些“字节序”）
另一个假设是：
- 它接收实际的浮点数，即乘法不会使用数学硬件/库可以在内部使用的更高精度值
#include <stdio.h>
#include <string.h>
#include <limits.h>
#include <assert.h>

#define C_ASSERT(expr) extern char CAssertExtern[(expr)?1:-1]

#if UINT_MAX >= 0xFFFFFFFF
typedef unsigned int uint32;
#else
typedef unsigned long uint32;
#endif
typedef unsigned long long uint64;

C_ASSERT(CHAR_BIT == 8);
C_ASSERT(sizeof(uint32) == 4);
C_ASSERT(sizeof(float) == 4);

int ProductIsOne(float f1, float f2)
{
  uint32 m1, m2;
  int e1, e2, s1, s2;
  int e;
  uint64 m;

  // Make sure floats are 32-bit IEE754 and
  // reinterpreted as integers as we expect
  {
    static const float testf = 8873283.0f;
    uint32 testi;
    memcpy(&testi, &testf, sizeof(testf));
    assert(testi == 0x4B076543);
  }

  memcpy(&m1, &f1, sizeof(f1));
  s1 = m1 >= 0x80000000;
  m1 &= 0x7FFFFFFF;
  e1 = m1 >> 23;
  m1 &= 0x7FFFFF;
  if (e1 > 0) m1 |= 0x800000;

  memcpy(&m2, &f2, sizeof(f2));
  s2 = m2 >= 0x80000000;
  m2 &= 0x7FFFFFFF;
  e2 = m2 >> 23;
  m2 &= 0x7FFFFF;
  if (e2 > 0) m2 |= 0x800000;

  if (e1 == 0xFF || e2 == 0xFF || s1 != s2) // Inf, NaN, different signs
    return 0;

  m = (uint64)m1 * m2;

  if (!m || (m & (m - 1))) // not a power of 2
    return 0;

  e = e1 + !e1 - 0x7F - 23 + e2 + !e2 - 0x7F - 23;
  while (m > 1) m >>= 1, e++;

  return e == 0;
}

const float testData[][2] =
{
  { .1f, 10.0f },
  { 0.5f, 2.0f },
  { 0.25f, 2.0f },
  { 4.0f, 0.25f },
  { 0.33333333f, 3.0f },
  { 0.00000762939453125f, 131072.0f }, // 2^-17 * 2^17
  { 1.26765060022822940E30f, 7.88860905221011805E-31f }, // 2^100 * 2^-100
  { 5.87747175411143754E-39f, 1.70141183460469232E38f }, // 2^-127 (denormalized) * 2^127
};

int main(void)
{
  int i;
  for (i = 0; i < sizeof(testData) / sizeof(testData[0]); i++)
    printf("%g * %g %c= 1\n",
           testData[i][0], testData[i][1],
           "!="[ProductIsOne(testData[i][0], testData[i][1])]);
  return 0;
}

输出结果（请查看ideone.com）：

0.1 * 10 != 1
0.5 * 2 == 1
0.25 * 2 != 1
4 * 0.25 == 1
0.333333 * 3 != 1
7.62939e-06 * 131072 == 1
1.26765e+30 * 7.88861e-31 == 1
5.87747e-39 * 1.70141e+38 == 1
+1。所以，二进制小数是精确的。你试过2^(-100)* 2^(+100)吗？ - Gangnus
@Gangnus：当然，如果是二进制的话，2的幂次方就是精确的。请参见ideone上更新的代码。我们甚至不需要在十进制中知道2^100或2^-100的所有有效数字。 - Alexey Frunze
我的意思是，在某些幂次以上，将2的幂放入浮点数的第二部分会出现问题。 - Gangnus
@Gangnus：超过最大指数，只有无穷大（代码在Inf和NaN上返回0）。低于最小指数，存在非规格化值（代码也处理它们）。请参见ideone上的另一个更新，演示了一种非规格化情况。 - Alexey Frunze
是的。我明白了。谢谢你。我在想操作那个边界附近应该设为 NaN 还是 0。 - Gangnus
0

令人惊讶的是，无论舍入规则如何，您都期望这两个版本的结果相同（两次错误或两次正确）！

很可能，在第一种情况下，当评估x==1/y时，在FPU寄存器中发生了向更高精度的提升，而z= 1/y实际上存储单精度结果。

其他贡献者已经解释了为什么5==1/0.2可能会失败，我不必重复。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接