С плавающей точкой C против C ++

1 min


Этот пост оказался не таким, как я ожидал. Вначале это был пост о причудах арифметики с плавающей запятой, и он превратился в пост о C против C ++.

Все началось с примера Жана-Майкла Мюллера, который я нашел в книге Джона Густафсона Конец ошибки, Задача состоит в том, чтобы оценить следующую функцию в 15, 16, 17 и 9999.

 begin {align *} e (x) & =  frac { exp (x) - 1} {x} \ q (x) & =  left | x -  sqrt {x ^ 2 + 1}  right | -  frac {1} {x +  sqrt {x ^ 2 + 1}} \ h (x) & = e (q (x) ^ 2) \  end {align *}

Вот е(0) определяется непрерывностью, чтобы быть 1.

Этот пример был чертовски разработан, чтобы показать, что может пойти не так с арифметикой с плавающей запятой, даже с высокой точностью.

Если вы непосредственно реализуете функции выше в C, вы получите 0 в качестве выходных данных, независимо от того, используете ли вы одинарную, двойную или даже четырехкратную точность, как показано в следующем коде. Тем не менее, правильный ответ в каждом случае равен 1. (По крайней мере, я каждый раз получал 0. Я не удивлюсь, если разные компиляторы дадут разные результаты. Подробнее об этом в ближайшее время.)

#include 
#include 
#include 

#define T __float128

T e(T x) {
  return x == 0. ? 1. : (exp(x) - 1.)/x;
}

T q(T x) {
  return fabs(x - sqrt(x*x + 1.)) - 1./(x + sqrt(x*x + 1.));
}

T h(T x) {
  return e( q(x)*q(x) );
}

int main() {

  int x() = {15, 16, 17, 9999};
  int i;
  for (i = 0; i < 4; i++) {
    printf("%fn", h((T) x(i)));    
  }
}

Вот T определяется для __float128 для четверной точности. Вы можете определить T быть float или же double для одинарной или двойной точности соответственно.

Вот программа на C ++, которая использует шаблоны вместо препроцессора C.

#include 
#include 
#include 

using namespace std;

// Example by Jean-Michel Muller
// Popularized by William Kahan

template  e(T x) {
  return x == 0. ? 1. : (exp(x) - 1.)/x;
}

template  q(T x) {
  return fabs(x - sqrt(x*x + 1.)) - 1./(x + sqrt(x*x + 1.));
}

template  h(T x) {
  return e( q(x)*q(x) );
}

int main() {

  int x() = {15, 16, 17, 9999};

  for (int i = 0; i < 4; i++) {
      cout << h(     float(x(i)) ) << endl;
      cout << h(    double(x(i)) ) << endl;
      cout << h(__float128(x(i)) ) << endl;
  }
}

Этот код по сути такой же, и все же он возвращает 1 в каждом случае! Я собираю оба с gcc 4.9.2. Для программы на C ++ я привожу аргумент -lstdc++ сказать это файл C ++.

Маленькая алгебра показывает, что функция Q(Икс) вернет 0 в точной арифметике, но не в арифметике с плавающей точкой. Возвращает чрезвычайно малое число, и, по-видимому, по крайней мере в версии C числитель (exp(x) - 1.)/x оценивает в 0. Возможно, C ++ использует другую математическую библиотеку, которая не оценивает в 0.

Я пытался заменить exp(x) - 1 с участием expm1(x), Это не помогло коду C, и это сломало код C ++. То есть обе программы вернули все нули. (Еще expm1 Вот.)

Между прочим, bc -l дает правильный результат, независимо от того, как мало вы установите scale,

define ee(x) { if (x == 0) return 1 else return (e(x) - 1)/x }
define abs(x) { if (x > 0) return x else return -x }
define q(x) { return abs(x - sqrt(x^2 + 1)) - 1/(x + sqrt(x^2 + 1)) }
define h(x) { return ee( q(x)^2 ) }

Похожие сообщения


0 Comments

Ваш e-mail не будет опубликован. Обязательные поля помечены *