尝试编写代码来查找机器 epsilon答案

【问题标题】：Trying to write a code for finding the machine epsilon尝试编写代码来查找机器 epsilon
【发布时间】：2011-03-29 13:39:57
【问题描述】：

我正在尝试找出 C 中各种浮点格式的精度级别（即 float、double 和 long double）。这是我目前正在使用的代码：

#include <stdio.h>
#define N 100000

int main(void)
{
   float max = 1.0, min = 0.0, test;
   int i;                              /* Counter for the conditional loop */

   for (i = 0; i < N; i++) {
      test = (max + min) / 2.0;
      if( (1.0 + test) != 1.0)         /* If too high, set max to test and try again */
     max = test;
  if( (1.0 + test) == 1.0)     /* If too low, set min to test and try again */
         min = test;
   }
   printf("The epsilon machine is %.50lf\n", max);
   return 0;
}

这给出了大约 ~2^-64 的值，正如预期的那样。但是，当我将减速更改为双打或“长双打”时，我得到相同的答案，我应该得到一个较小的值，但我没有。有人有什么想法吗？

【问题讨论】：

浮点数没有 23 位尾数吗？为什么你会期望 2^-64？

标签： c floating-point double floating-accuracy

【解决方案1】：

这取决于您所说的“精度级别”。

浮点数具有“常规”（正常）值，但也有特殊的、次正常的数字。如果您想找出不同的限制，C 标准有预定义的常量：

#include <math.h>
#include <stdio.h>
#include <float.h>

int main(void)
{
    printf("%30s: %g\n", "FLT_EPSILON", FLT_EPSILON);
    printf("%30s: %g\n", "FLT_MIN", FLT_MIN);
    printf("%30s: %g\n", "nextafterf(0.0, 1.0)", nextafterf(0.0, 1.0));
    printf("%30s: %g\n", "nextafterf(1.0, 2.0)-1", (nextafterf(1.0, 2.0) - 1.0f));
    puts("");
    printf("%30s: %g\n", "DBL_EPSILON", DBL_EPSILON);
    printf("%30s: %g\n", "DBL_MIN", DBL_MIN);
    printf("%30s: %g\n", "nextafter(0.0, 1.0)", nextafter(0.0, 1.0));
    printf("%30s: %g\n", "nextafter(1.0, 2.0)-1", (nextafter(1.0, 2.0) - 1.0));
    puts("");
    printf("%30s: %Lg\n", "LDBL_EPSILON", LDBL_EPSILON);
    printf("%30s: %Lg\n", "LDBL_MIN", LDBL_MIN);
    printf("%30s: %Lg\n", "nextafterl(0.0, 1.0)", nextafterl(0.0, 1.0));
    printf("%30s: %Lg\n", "nextafterl(1.0, 2.0)-1", (nextafterl(1.0, 2.0) - 1.0));
    return 0;
}

上述程序为每种类型打印 4 个值：

1 与该类型中大于 1 的最小值之间的差 (TYPE_EPSILON)，
给定类型中的最小正标准化值 (TYPE_MIN)。这不包括subnormal numbers，
给定类型中的最小正值 (nextafter*(0...))。这包括次正规数，
大于 1 的最小数。这与 TYPE_EPSILON 相同，但计算方式不同。

根据您所说的“精确度”的含义，以上任何一项对您都有用或没有任何用处。

这是上面程序在我电脑上的输出：

               FLT_EPSILON: 1.19209e-07
                   FLT_MIN: 1.17549e-38
      nextafterf(0.0, 1.0): 1.4013e-45
    nextafterf(1.0, 2.0)-1: 1.19209e-07

               DBL_EPSILON: 2.22045e-16
                   DBL_MIN: 2.22507e-308
       nextafter(0.0, 1.0): 4.94066e-324
     nextafter(1.0, 2.0)-1: 2.22045e-16

              LDBL_EPSILON: 1.0842e-19
                  LDBL_MIN: 3.3621e-4932
      nextafterl(0.0, 1.0): 3.6452e-4951
    nextafterl(1.0, 2.0)-1: 1.0842e-19

【讨论】：

干杯，我很高兴C为此内置了东西。但是，我的任务是编写代码来找到大于可以用浮点数表示的最小值。那里的第一个数字：1.19209e-07 是我所期望的，但由于某种原因，我的代码没有给我这个数字。非常感谢
@Jack：好的。然后您应该确保计算中使用的所有浮点数都是float 值。所以我不会做1.0 + test != 1.0，而是做：float try = 1.0 + test; if (try != 1.0)等。
干杯，伙计们，现在开始给出明智的答案

【解决方案2】：

猜猜为什么你会得到相同的答案：

if( (1.0 + test) != 1.0)

这里的 1.0 是一个双精度常数，因此它将浮点数提升为双精度，并将加法作为双精度执行。您可能希望在此处声明一个临时浮点数来执行加法，或者将这些浮点数设为数字常量 (1.0fIIRC)。

您还可能陷入临时浮点数超精度问题，可能需要强制将中间值存储在内存中以降低到正确的精度。

这是重做范围搜索方法但以正确类型计算测试的快速方法。不过，我得到的答案有点过大。

#include <stdio.h>
#define N 100000
#define TYPE float

int main(void)
{
   TYPE max = 1.0, min = 0.0, test;
   int i;

   for (i = 0; i < N; i++)
   {
      TYPE one_plus_test;

      test = (max + min) / ((TYPE)2.0);
      one_plus_test = ((TYPE)1.0) + test;
      if (one_plus_test == ((TYPE)1.0))
      {
         min = test;
      }
      else
      {
         max = test;
      }
   }
   printf("The epsilon machine is %.50lf\n", max);
   return 0;
}

【讨论】：

如何将其“投射”为浮点数？我会试一试，看看会发生什么
是的，我试过了，但它仍然给我和 2^-64 的 epsilon 值
好的，尝试将其存储到 volatile float 变量中，然后从中读取。即：volatile float tmp = 1.0 + test; if (tmp == 1.0) ...

【解决方案3】：

我不确定您的算法应该如何工作。这个（C++）给出了正确答案：

#include <iostream>

template<typename T>
int epsilon() {
    int pow = 0;
    T eps = 1;
    while (eps + 1 != 1) {
        eps /= 2;
        --pow;
    }
    return pow + 1;
}

int main() {
    std::cout << "Epsilon for float: 2^" << epsilon<float>() << '\n';
    std::cout << "Epsilon for double: 2^" << epsilon<double>() << '\n';
}

这会计算最小值，使得当加到 1 时仍可与 1 区分开来。

输出：

Epsilon for float: 2^-23
Epsilon for double: 2^-52

【讨论】：

恐怕我不懂任何 C++
如果你知道 C 应该不难理解。模板只允许我写 T 并用 float 或 double 代替它。而且打印效果不同，但不用担心。

【解决方案4】：

IEEE 754 浮点格式的特性是，当重新解释为相同宽度的二进制补码整数时，它们在正值上单调递增，在负值上单调递减（参见 32 位浮点数的二进制表示）。它们还具有 0

typedef union {
  long long i64;
  double d64;
} dbl_64;

double machine_eps (double value)
{
    dbl_64 s;
    s.d64 = value;
    s.i64++;
    return s.d64 - value;
}

来自https://en.wikipedia.org/wiki/Machine_epsilon

【讨论】：

【解决方案5】：

我想补充一点，您可以使用 long double 从浮点计算中获得最高精度。

要将其应用于@Rup 的解决方案，只需将TYPE 更改为long double 并将printf 语句更改为：

printf("The epsilon machine is %.50Lf\n", max);

这是我机器上使用 float 的 Epsilon：

0.00000005960465188081798260100185871124267578125000

并使用long double：

0.00000000000000000005421010862427522170625011179761

差异很大。

【讨论】：

【解决方案6】：

此类代码的一个问题是编译器会将浮点变量加载到微处理器的浮点寄存器中。如果您的微处理器只有双精度浮点寄存器，float 和 double 的精度将相同。

您需要找到一种方法来强制编译器在每两次计算之间将浮点值存储回内存（存储到正确类型的变量中）。这样，它必须丢弃寄存器的额外精度。但是今天的编译器在优化你的代码方面很聪明。所以这可能很难实现。

【讨论】：

为什么不在调试模式或类似模式下编译，不执行优化？