使用 epsilon 将双精度数与零进行比较

Posted 2023-03-11

【中文标题】使用 epsilon 将双精度数与零进行比较

【英文标题】：Compare double to zero using epsilon

【发布时间】：2012-11-21 20:15:36

【问题描述】：

今天，我翻阅了一些 C++ 代码（由其他人编写），发现了这个部分：

double someValue = ...
if (someValue <  std::numeric_limits<double>::epsilon() && 
    someValue > -std::numeric_limits<double>::epsilon()) 
  someValue = 0.0;

我正在尝试弄清楚这是否有意义。

epsilon() 的文档说：

该函数返回 1 和大于 1 的最小值之间的差值，该值可以 [用双精度数] 表示。

这是否也适用于 0，即epsilon() 是大于 0 的最小值？或者0 和0 + epsilon 之间是否有可以用double 表示的数字？

如果不是，那么比较不等于someValue == 0.0？

【问题讨论】：

1 附近的 epsilon 很可能比 0 附近的要高得多，因此可能会有 0 到 0+epsilon_at_1 之间的值。我猜这部分的作者想使用一些小的东西，但他不想使用魔法常数，所以他只是使用了这个本质上是任意值。

比较浮点数很困难，甚至鼓励使用 epsilon 或阈值。请参考：cs.princeton.edu/introcs/91float和cygnus-software.com/papers/comparingfloats/comparingfloats.htm

第一个链接是 403.99999999

epsilon() 是最小的正值。因此，如果我们假设 epsilon() 是 e，我们得到 1+e != 1，所以是的，epsilon 是大于 0 的最小值，并且在 0 和 0 + e 之间没有数字

IMO，在这种情况下，numeric_limits<>::epsilon 的使用具有误导性且无关紧要。如果实际值与 0 的差异不超过某个 ε，我们想要假设 0。并且 ε 应该根据问题规范而不是与机器相关的值来选择。我怀疑当前的 epsilon 是无用的，因为即使是几个 FP 操作也可能累积比这更大的错误。

【参考方案1】：

假设 64 位 IEEE double，则有 52 位尾数和 11 位指数。让我们分解一下：

1.0000 00000000 00000000 00000000 00000000 00000000 00000000 × 2^0 = 1

大于1的最小可表示数：

1.0000 00000000 00000000 00000000 00000000 00000000 00000001 × 2^0 = 1 + 2^-52

因此：

epsilon = (1 + 2^-52) - 1 = 2^-52

0 和 epsilon 之间有数字吗？很多...例如最小正可表示（正常）数是：

1.0000 00000000 00000000 00000000 00000000 00000000 00000000 × 2^-1022 = 2^-1022

实际上在 0 和 epsilon 之间有 (1022 - 52 + 1)×2^52 = 4372995238176751616 数，占所有可表示的正数的 47%...

【讨论】：

太奇怪了，你可以说“47% 的正数”:)

@configurator：不，你不能这么说（不存在“自然”有限度量）。但你可以说“47% 的正可表示数字”。

@ybungalobill 我想不通。指数有 11 位：1 个符号位和 10 个值位。为什么 2^-1022 而不是 2^-1024 是最小的正数？

@PavloDyban：仅仅因为指数没有有符号位。它们被编码为偏移量：如果编码的指数是0 <= e < 2048，则尾数乘以 2 的 e - 1023 的幂。例如。 2^0 的指数编码为e=1023，2^1 编码为e=1024，2^-1022 编码为e=1。 e=0 的值是为次正规和实零保留的。

@PavloDyban：2^-1022 也是最小的正常数。最小的数字实际上是0.0000 00000000 00000000 00000000 00000000 00000000 00000001 × 2^-1022 = 2^-1074。这是次正规的，意味着尾数部分小于 1，所以用指数e=0 编码。

【参考方案2】：

X 和X 的下一个值之间的差值根据X 而不同。epsilon() 只是1 和1 的下一个值之间的差值。0 与0 的下一个值之间的差异不是epsilon()。

相反，您可以使用std::nextafter 将双精度值与0 进行比较，如下所示：

bool same(double a, double b)

  return std::nextafter(a, std::numeric_limits<double>::lowest()) <= b
    && std::nextafter(a, std::numeric_limits<double>::max()) >= b;


double someValue = ...
if (same (someValue, 0.0)) 
  someValue = 0.0;

【讨论】：

+1 用于提及nextafter；但请注意，这种用法不太可能符合程序员的意图。假设 64 位 IEEE 754，在您的示例中 same(0, 1e-100) 返回 false，这可能不是程序员想要的。程序员可能宁愿想要一些小的阈值来测试相等性，例如+/-1e-6 或 +/-1e-9，而不是 +/-nextafter。

【参考方案3】：

假设我们正在使用适合 16 位寄存器的玩具浮点数。有一个符号位、一个 5 位指数和一个 10 位尾数。

这个浮点数的值是尾数，解释为二进制十进制值，乘以 2 的指数次方。

在 1 附近，指数等于 0。所以尾数的最小位数是1024的一部分。

接近 1/2 的指数是负一，所以尾数的最小部分是原来的一半。使用 5 位指数可以达到负 16，此时尾数的最小部分相当于 32m 的一部分。在负 16 指数处，该值大约是 32k 的一部分，比我们上面计算的大约 1 的 epsilon 更接近于零！

现在这是一个玩具浮点模型，它不能反映真实浮点系统的所有怪癖，但反映小于 epsilon 的值的能力与真实浮点值相当相似。

【参考方案4】：

此外，拥有这样一个函数的一个很好的原因是删除“非正规”（那些不能再使用隐含的前导“1”并具有特殊 FP 表示的非常小的数字）。你为什么想做这个？因为有些机器（特别是一些较旧的 Pentium 4s）在处理非规范化时会变得非常非常慢。其他人只是变得有点慢。如果您的应用程序并不真正需要这些非常小的数字，那么将它们清零是一个很好的解决方案。考虑这一点的好地方是任何 IIR 滤波器或衰减函数的最后一步。

另请参阅：Why does changing 0.1f to 0 slow down performance by 10x?

和http://en.wikipedia.org/wiki/Denormal_number

【讨论】：

这会删除更多的数字，而不仅仅是非规范化的数字。它将普朗克常数或电子质量更改为零，如果您使用这些数字，将会得到非常非常错误的结果。

【参考方案5】：

对于 IEEE 浮点，在最小的非零正值和最小的非零负值之间，存在两个值：正零和负零。测试一个值是否在最小的非零值之间，相当于测试是否与零相等；然而，赋值可能会产生影响，因为它会将负零变为正零。

可以想象，浮点格式可能具有介于最小有限正负值之间的三个值：正无穷小、无符号零和负无穷小。我不熟悉实际上以这种方式工作的任何浮点格式，但这样的行为将是完全合理的，并且可以说比 IEEE 的更好（也许不够好到值得添加额外的硬件来支持它，但在数学上 1 /(1/INF)、1/(-1/INF) 和 1/(1-1) 应该代表三种不同的情况，说明三个不同的零）。我不知道是否有任何 C 标准会要求有符号的无穷小（如果存在）必须比较为零。如果他们不这样做，上面的代码可以有效地确保例如将一个数字重复除以 2 最终会产生零，而不是停留在“无穷小”上。

【讨论】：

“1/(1-1)”（来自您的示例）不是无穷大而不是零吗？

数量 (1-1)、(1/INF) 和 (-1/INF) 都表示零，但是将正数除以它们中的每一个理论上应该会产生三个不同的结果 ( IEEE 数学认为前两个相同）。

【参考方案6】：

可以使用以下程序打印数字 (1.0, 0.0, ...) 周围的 epsilon 近似值（可能的最小差异）。它打印以下输出：epsilon for 0.0 is 4.940656e-324epsilon for 1.0 is 2.220446e-16 稍微思考一下就清楚了，我们用于查看其 epsilon 值的数字越小，epsilon 就越小，因为指数可以调整为该数字的大小。

#include <stdio.h>
#include <assert.h>
double getEps (double m) 
  double approx=1.0;
  double lastApprox=0.0;
  while (m+approx!=m) 
    lastApprox=approx;
    approx/=2.0;
  
  assert (lastApprox!=0);
  return lastApprox;

int main () 
  printf ("epsilon for 0.0 is %e\n", getEps (0.0));
  printf ("epsilon for 1.0 is %e\n", getEps (1.0));
  return 0;

【讨论】：

你检查了哪些实现？ GCC 4.7 绝对不是这种情况。

【参考方案7】：

测试肯定和someValue == 0不一样。浮点数的整个想法是它们存储一个指数和一个有效数。因此，它们表示具有一定数量的二进制有效数字精度的值（在 IEEE 双精度的情况下为 53）。可表示的值在 0 附近比在 1 附近更密集。

要使用更熟悉的十进制系统，假设您使用指数存储一个“到 4 位有效数字”的十进制值。那么下一个大于1的可表示值是1.001 * 10^0，而epsilon是1.000 * 10^-3。但是1.000 * 10^-4 也是可以表示的，假设指数可以存储-4。你可以相信我的话，IEEE double 可以存储的指数小于epsilon 的指数。

您无法仅从这段代码中判断将epsilon 专门用作绑定是否有意义，您需要查看上下文。可能epsilon 是对产生someValue 的计算中的错误的合理估计，也可能不是。

【讨论】：

好点，但即使是这种情况，更好的做法是将错误限制在一个合理命名的变量中并在比较中使用它。就目前而言，它与魔法常数没有什么不同。

也许我的问题应该更清楚：我没有质疑 epsilon 是否足够大的“阈值”来涵盖计算错误，而是这个比较是否等于 someValue == 0.0。

【参考方案8】：

有些数字存在于 0 和 epsilon 之间，因为 epsilon 是 1 和可以在 1 以上表示的下一个最高数字之间的差，而不是 0 和可以在 0 以上表示的下一个最高数字之间的差（如果它是，该代码将做的很少）：-

#include <limits>

int main ()

  struct Doubles
  
      double one;
      double epsilon;
      double half_epsilon;
   values;

  values.one = 1.0;
  values.epsilon = std::numeric_limits<double>::epsilon();
  values.half_epsilon = values.epsilon / 2.0;

使用调试器，在 main 结束时停止程序并查看结果，您会发现 epsilon / 2 与 epsilon、0 和 1 不同。

因此，此函数采用 +/- epsilon 之间的值并将它们设为零。

【参考方案9】：

由于尾​​数和指数部分，您不能将此应用于 0。 由于指数，您可以存储比 epsilon 小的数字， 但是当您尝试执行 (1.0 - "very small number") 之类的操作时，您将得到 1.0。 Epsilon 不是值的指标，而是值精度的指标，以尾数表示。 它显示了我们可以存储多少个正确的后继十进制数字。

【参考方案10】：

我认为这取决于您计算机的precision。 看看这个table：可以看到如果你的epsilon用double表示，但是你的精度更高，比较不等于

someValue == 0.0

好问题！

【参考方案11】：

假设系统无法区分 1.000000000000000000000 和 1.000000000000000000001。即 1.0 和 1.0 + 1e-20。您认为在-1e-20 和+1e-20 之间还有一些值可以表示吗？

【讨论】：

除了零，我认为-1e-20和+1e-20之间没有值。但仅仅因为我认为这并不成立。

@SebastianKrysmanski：这不是真的，在 0 和 epsilon 之间有很多浮点值。因为它是浮点点，而不是定点。

不同于零的最小可表示值受分配用于表示指数的位数限制。因此，如果 double 具有 11 位指数，则最小的数字将是 1e-1023。