试图编写查找机器epsilon的代码

猜猜你为什么得到同样的答案:

if( (1.0 + test) != 1.0)

这里1.0是一个双常数,所以它将把浮点提升为双常数,并将加法作为双常数执行。您可能想在这里声明一个临时的浮点来执行加法,或者使这些浮点数字常量( 1.0f IIRC)。

您可能还陷入了临时浮动问题中的过度精度问题,可能需要强制它将中间值存储在内存中,以降低到正确的精度。

这里有一个快速的步骤,可以重新使用你的范围搜索方法,但是要用正确的类型计算测试。不过,我得到的答案有点太大了。

#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;
}

这取决于你所说的“精确水平”。

浮点数有“常规”(正常)值,但也有特殊的、次正常的值。如果要找出不同的限制,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的最小值之间的差异( 类型 _EPSILON )
• 给定类型中的最小正标准化值( 类型 _MIN )这不包括 subnormal numbers ,
• 给定类型中的最小正值( nextafter * (0 … ) )这包括次正规数,
• 最小值大于1。这和 类型 εε 但计算方式不同。

根据你所说的“精确”的意思,上述任何一项或任何一项对你都是有用的。

以下是我的计算机上上述程序的输出:

               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++)给出正确的答案:

#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

IEEE754浮点格式具有这样的特性:当重新解释为相同宽度的两个补整数时,它们单调地增加到正值,单调地减少到负值(参见32位浮点的二进制表示)。它们还具有0<f(x)<、和f(x+1)f(x)f(x)f(x1)(其中f(x)是上述x的整数重新解释)。在允许使用类型punning并总是使用IEEE754-1985的语言中,我们可以利用它在恒定时间内计算机器epsilon。例如,在c中:

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

我想补充一点,您可以通过使用 long double .

要将其应用于@rup的解决方案,只需更改 TYPE 到 长双倍 以及 printf 声明:

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

这是我机器上的epsilon float :

0.00000005960465188081798260100185871124267578125000

并使用 长双倍 :

0.00000000000000000005421010862427522170625011179761

两者之间的差异非常显著。

这种代码的一个问题是编译器将把浮点变量加载到微处理器的浮点寄存器中。如果你的微处理器只有双精度浮点寄存器,那么它的精度将与 float 和 double .

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