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

Question

我想确定（在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; 

将输出：“他们是......”，这是正确的。
这是为什么发生？

Answer 1

浮动精度问题

您这里有两个问题，但都来自同根

不能精确地进行比较花车。你不能精确地减去或分割它们。你无法准确地为他们计算什么。对它们的任何操作都可能（并且几乎总是）给结果带来一些错误。即使a=0.2f不是一个精确的操作。这里的其他答案的作者很好地解释了更深层次的原因。 _{（我的感谢和投票给他们。）}

这里是你的第一个也是更简单的错误。你永远不应该，从未，从未，从未，NEVER使用它们==或等值的任何语言。代替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，y

因此，接下来的两个浮动操作规则：尽量不要分大于由较小的一个评估师和上帝保存你减去接近值之后。

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

• 通过创始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){... 

这是所有。

---

P.S.如果你真的需要它所有在一行，使用方法：

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

但它是坏的风格。我不会建议。

P.P.S.在第二个例子中，编译器优化代码，以便在运行任何代码之前将z设置为5。因此，即使对于浮标，检查5对5也是如此。

Answer 2

您将不得不准确定义两个逼近是乘法逆的含义。否则，你将不知道你应该测试什么。

0.2没有确切的二进制表示。如果你存储的数字没有精确的表示，精确度有限，你将无法得到完全正确的答案。

同样的事情发生在十进制。例如，1/3没有确切的十进制表示形式。您可以将其存储为.333333。但是，你有一个问题。 3和.333333乘法反转？如果你乘以它们，你会得到.999999。如果你想让答案为“是”，你必须创建一个乘法倒数的测试，它不像乘法和测试的等于1那么简单。

同样的事情发生在二进制文件中。

Answer 3

的问题是，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 is suggesting。

Answer 4

引人注目的是，无论舍入规则是什么，你都期望两个版本的结果是相同的（无论是两次错误还是两次）！

很可能，在第一种情况下，在评估x == 1/y时，FPU寄存器中的精度更高，而z = 1/y确实存储单精度结果。

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

Answer 5

其他回复中的讨论非常棒，所以我不会重复其中的任何内容，但是没有代码。下面是一些代码，用于实际检查一对浮点数是否在乘法时精确到1.0。

该代码使得几个假设/断言（其通常满足x86平台上）：
- float的是32位二进制（AKA single precision）IEEE-754
- 要么int的或long的有32位（我决定不依靠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