Пытаюсь написать код для поиска машины эпсилон - PullRequest
Новая
       21

Пытаюсь написать код для поиска машины эпсилон

5 голосов
/ 13 августа 2010

Я пытаюсь выяснить уровень точности для различных форматов с плавающей запятой в 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, как и ожидалось. Однако, когда я изменяю замедления на удвоения или «длинные удвоения», я получаю один и тот же ответ, я должен получить меньшее значение, но не получаю. У кого-нибудь есть идеи?

Ответы [ 6 ]

9 голосов
/ 13 августа 2010

Это зависит от того, что вы подразумеваете под "уровнем точности".

Числа с плавающей запятой имеют "обычные" (нормальные) значения, но есть и специальные, ненормальные числа.Если вы хотите узнать другие пределы, стандарт 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).Сюда не входят субнормальные числа ,
  • минимальное положительное значение для данного типа (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
2 голосов
/ 22 мая 2016

Форматы IEEE 754 имеют свойство, заключающееся в том, что при повторной интерпретации как целое число с дополнением до двух одинаковой ширины они монотонно увеличиваются по положительным значениям и монотонно уменьшаются по отрицательным значениям (см. Двоичное представление 32-разрядных чисел с плавающей точкой).У них также есть свойство, что 0 <| f (x) |<∞ и | f (x + 1) - f (x) |≥ | f (x) - f (x − 1) |(где f (x) - вышеупомянутая целочисленная реинтерпретация x).В языках, которые допускают штампование типов и всегда используют IEEE 754-1985, мы можем использовать это для вычисления машинного эпсилона в постоянное время.Например, в C: </p> 

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

From https://en.wikipedia.org/wiki/Machine_epsilon

2 голосов
/ 13 августа 2010

Я не уверен, как ваш алгоритм должен работать.Этот (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.1007 *

2 голосов
/ 13 августа 2010

Угадайте, почему вы получаете один и тот же ответ: 

if( (1.0 + test) != 1.0)

Здесь 1.0 - двойная константа, поэтому она повышает ваш float до двойного и выполняет сложение как двойное.Возможно, вы захотите объявить временное значение с плавающей точкой здесь, чтобы выполнить сложение, или сделать эти числовые константы с плавающей точкой (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;
}
1 голос
/ 25 января 2016

Я хотел бы добавить, что вы можете получить самую высокую точность вычислений с плавающей запятой, используя long double.

Чтобы применить это к решению @ Rup, просто измените TYPE на long double и оператор printf на: 

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

Это Epsilon на моей машине, используя float: 

0.00000005960465188081798260100185871124267578125000

И с использованием long double: 

0.00000000000000000005421010862427522170625011179761

Разница довольно значительная.

0 голосов
/ 13 августа 2010

Проблема с таким кодом состоит в том, что компилятор загружает переменные с плавающей запятой в регистры микропроцессора с плавающей запятой. И если ваш микропроцессор имеет только регистры с плавающей точкой двойной точности, точность будет одинаковой для float и double.

Вам нужно будет найти способ заставить компилятор сохранять значение с плавающей запятой обратно в память между каждыми двумя вычислениями (в переменную правильного типа). Таким образом, он должен отбросить дополнительную точность регистров. Но сегодняшние компиляторы умны в оптимизации вашего кода. Так что это может быть трудно достичь.

...