如何确定双精度浮点数的最小可察觉变化

Question

如何确定双精度浮点数的最小可察觉变化

7

我有一个问题，需要确定给定双精度变量 v 的最小值 eps，使得：

v+eps != v

请注意，这不是典型的问题清单任务，因为eps取决于任意数字v。

不应该通过在for循环中寻找此值来完成。是否有快速的方法来完成这个任务，例如通过位移？独立于编译器、优化标志和平台...

感谢您的回答。

- choosyg

可能是重复的问题：如何找到给定数字的最近的下一个/上一个双精度值（使用numeric_limits :: epsilon）。 - rubenvb

@BoBTFish，这完全没有回答问题。 - rubenvb

@rubenvb 是的，我开始写了，但是我还没有完全理解如何将其转换为 epsilon。所以我去读了一些资料。在我弄清楚自己想要表达什么之前，我会将其删除。 - BoBTFish

1

问题和标题要求不同。如果d是从（正的）v到下一个更大的可表示值的最小变化，e是最小的（正的）值，使得计算v+e不会产生v，则e是d的一半或稍微大一些，这取决于v的有效数字的低位。这是因为，在浮点运算中，数学结果只有到达下一个可表示值的一半就足以导致向上舍入。（如果恰好到一半，如果v的低位是偶数，则向下舍入，如果是奇数，则向上舍入，这就是为什么e可能略微大于d的原因。） - Eric Postpischil

你应该清楚地指定你想要的是d还是e，并更改标题和问题文本以匹配。（注意：这假设使用常见的四舍五入模式IEEE 754算术。） - Eric Postpischil

3个回答

1

我会使用类型转换（type punning）：

double
epsFor( double x )
{
    union
    {
        double d;
        unsigned long long i;
    } tmp;
    tmp.d = x;
    ++ tmp.i;
    double results = tmp.d - x;
    return results;
}

（从正式意义上讲，这是未定义行为，但在实践中，我不知道有现代编译器会失败。）

编辑：

请注意，C++允许中间表达式具有过度精度；由于我们关心的是精确结果，如果您直接在表达式中使用原始发布的函数而不是将其分配给double，则该函数可能会给出错误的结果。我已经添加了一个分配到函数中以避免这种情况，但请注意，很多编译器在这方面并不符合标准，至少默认情况下不符合标准。（g++是其中一个很好的例子，你需要特殊选项才能获得符合标准的行为，至少在打开优化时。如果您正在使用g++，则必须指定-ffloat-store选项，如果要获得正确的结果。）

- James Kanze

1

一定会因为未定义的行为而投反对票。您可以在GCC和Clang中设置优化标志，这将导致其出现问题。 - Puppy

@JonasWielicki 任何使用都明显依赖于实现。我使用了最常见的平台（Intel和AMD、Sparc、Power PC等）的适当类型。在这种情况下，我不确定 unint64_t 是否有用：它会导致在一些奇特的平台上（例如Unisys大型机）编译失败，这可能是一个优势。（或者不是。我不知道那些平台上的 long long 是什么，使用 uint64_t 可能会导致编译失败，尽管上述代码可以正常工作。） - James Kanze

@DeadMG G++过去一直保证了这一点。如果您愿意，可以使用memcpy来解决未定义的行为；标准的“意图”是reinterpret_cast应该与类型游戏一起工作，但实际上，g++只保证使用union（这是最常见的方法）可以实现它。由于上述内容显然仅适用于IEEE（和类似表示法，其中在尾数中使用隐式MSB），因此我只关心实际实践中发生的情况。 - James Kanze

"-ffloat-store"并不能完全解决高精度中间结果的问题，因为它会导致双重舍入。解决方案是采用更高精度，并将其纳入语言语义的定义之中。这就是C99编译器定义的FLT_EVAL_METHOD所用之处。最近的C++标准也定义了FLT_EVAL_METHOD，并具有与C中相同的含义。 - Pascal Cuoq

@PascalCuoq 我还没有验证。_文档_称-fexcess-precions仅适用于C，这意味着不适用于C ++。但在C++中使用它肯定也是有道理的。 - James Kanze

显示剩余6条评论

0

eps = std::numeric_limits<double>::epsilon() * v;

- Ivan Ishchenko

这里并不是正确的答案，因为存在四舍五入的问题...例如 v + veps(2/3)。 - Sneftel

这就是 1.0 和下一个可表示的 double 之间的差异。对于更大的 double，这个差距会更大。 - BoBTFish

1

@BoBTFish：epsilon()就是这样，所以你应该将它乘以v。 - Ivan Ishchenko

2

间隔并非线性增加。 - Puppy

1

考虑到OP不知道他想要ULP(v)还是它的一半，我认为这个答案已经足够好了。v * (1 + epsilon) - v可以计算出ULP的一半，而不使用nextafter()或类型转换，但有一半的时间会计算出2ULPs。也许像v + (v * epsilon / 2) - v这样的东西可以在所有情况下都起作用，但是按照现在的写法，我仍然担心它在某些情况下是错误的。 - Pascal Cuoq

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- rubenvb · Accepted Answer

C99函数nextafter是您需要的。或者使用Boost.Math的nextafter。根据定义，这是实现定义的（它依赖于内存中double的内部表示）。

有关此处提供的所有方法的比较，请参见实时演示以查看其他解决方案如何失败。

以下是测试代码作为参考，如果您想在自己的系统上运行它：

#include <cmath>
#include <cfloat>
#include <limits>
#include <iostream>
using std::cout;
#include <iomanip>
using std::setprecision;

#include <boost/math/special_functions/next.hpp>

double
epsFor( double x )
{
  union
  {
    double d;
    unsigned long long i;
  } tmp;
  tmp.d = x;
  ++ tmp.i;
  return tmp.d - x;
}

void test(double d)
{
  double d1 = std::nextafter(d,DBL_MAX);
  double d2 = d+std::numeric_limits<double>::epsilon() * d;
  double d3 = d+epsFor(d);
  double d4 = boost::math::nextafter(d, DBL_MAX);
  cout << setprecision(40)
       << "For value of d = " << d << '\n'
       << " std::nextafter: " << d1 << '\n'
       << " Boost solution: " << d4 << '\n'
       << " undefined beh.: " << d3 << '\n'
       << " numeric_limits: " << d2 << '\n';
}

int main()
{
  test(0.1);
  test(986546357654.354687);
}