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;

将输出：“他们不是......”这是错误的，这段代码：

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;

将输出：“他们是……”这是正确的。
为什么会这样？

4

5 回答 5

36

浮点精度问题

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

您无法精确比较浮点数。您不能精确地减去或除以它们。你不能精确地为他们计算任何东西。对它们进行的任何操作都可能（并且几乎总是会）给结果带来一些错误。甚至a=0.2f不是一个精确的操作。此处其他答案的作者很好地解释了其更深层次的原因。（我对此表示感谢并投票给他们。）

这是您的第一个也是更简单的错误。你不应该，从不，从不，从不，永远不要在它们上使用 == 或任何语言的等价物。

代替a==b，使用Abs(a-b)<HighestPossibleError代替。

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

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，是的

因此，浮点运算的下两个规则：尽量不要将较大的估值器除以较小的估值器，然后上帝保佑你不要减去接近的值。

有两种简单的方法可以避免这个问题。

  • 通过确定 x,y 的绝对值较大，然后将 1 除以较大的值，然后再减去较小的值。

  • 如果您想比较1/y against x，而您正在使用字母而不是值，并且您的操作没有错误，请将比较的两边乘以 y 并且您有1 against x*y。（通常你应该在那个操作中检查符号，但这里我们使用 abs 值，所以它是干净的。）结果比较根本没有除法。

用更短的方式：

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

我们已经知道应该这样比较1 against x*y：

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

就这些。

PS如果你真的需要它在一行，使用：

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

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

PPS 在您的第二个示例中，编译器会优化代码，以便在运行任何代码之前将 z 设置为 5。因此，即使是浮点数，检查 5 对 5 也有效。

于 2012-02-03T23:33:39.980 回答
13

问题是0.2不能用二进制精确表示，因为它的二进制扩展有无限位数：

 1/5: 0.0011001100110011001100110011001100110011...

这类似于1/3无法以十进制精确表示。由于x存储在float具有有限位数的 a 中，因此这些数字将在某些时候被截断，例如：

   x: 0.0011001100110011001100110011001

问题的出现是因为 CPU 内部经常使用更高的精度，所以当你刚刚计算1/y时，结果会有更多的数字，当你加载x比较它们时，x会被扩展以匹配 CPU 的内部精度。

 1/y: 0.0011001100110011001100110011001100110011001100110011
   x: 0.0011001100110011001100110011001000000000000000000000

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

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

   x: 0.0011001100110011001100110011001
   z: 0.0011001100110011001100110011001

许多编译器都有开关，您可以启用它来强制在每一步截断中间值以保持一致性，但是通常的建议是避免在浮点值之间进行直接比较，而是检查它们的差异是否小于某个 epsilon 值，即Gangnus 的建议是什么。

于 2012-02-03T23:52:13.207 回答
5

您必须精确定义两个近似值是乘法逆的含义。否则，您将不知道您应该测试的是什么。

0.2没有精确的二进制表示。如果您以有限的精度存储没有精确表示的数字，您将不会得到完全正确的答案。

同样的事情发生在十进制。例如，1/3没有精确的十进制表示。您可以将其存储为.333333. 但是你有一个问题。是3和.333333乘法逆元吗？如果将它们相乘，则得到.999999. 如果您希望答案为“是”，则必须为乘法逆元创建一个测试，而不是像乘法和测试等于 1 那样简单。

二进制也会发生同样的事情。

于 2012-02-03T23:42:14.143 回答
2

其他回复中的讨论很棒，所以我不会重复其中的任何一个，但没有代码。这里有一些代码来实际检查一对浮点数在相乘时是否正好给出 1.0。

代码做了一些假设/断言（通常在 x86 平台上遇到）：
-float是 32 位二进制（AKA single precision）IEEE-754
- 要么int'要么 'long是 32 位（我决定不依赖可用性of uint32_t)
-memcpy()将浮点数复制到整数/长整数，使得 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
于 2012-02-16T15:14:41.647 回答
0

令人吃惊的是，无论四舍五入规则是什么，您都希望两个版本的结果相同（两次错误或两次正确）！

最有可能的是，在第一种情况下，在计算 x==1/y 时会提升 FPU 寄存器的精度，而 z=1/y 真正存储的是单精度结果。

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

于 2012-02-14T11:49:07.627 回答