C++ sqrt гарантированная точность, верхняя / нижняя граница

Question

Я должен проверить неравенство, содержащее квадратные корни. Чтобы избежать неверных результатов из-за неточностей с плавающей точкой и округления, я использую std::nextafter() чтобы получить верхнюю / нижнюю границу:

#include <cfloat> // DBL_MAX
#include <cmath> // std::nextafter, std::sqrt

double x = 42.0; //just an example number
double y = std::nextafter(std::sqrt(x), DBL_MAX);

а) y*y >= x гарантировано использование компилятора GCC?

б) Будет ли это работать для других операций, таких как + - * / или даже std::cos() а также std::acos()?

в) Есть ли лучшие способы получить верхние / нижние границы?

Обновление: я прочитал, что это не гарантируется стандартом C++, но должно работать в соответствии с IEEE-754. Будет ли это работать с компилятором GCC?

Answer 1

В общем случае операции с плавающей запятой приводят к некоторой ошибке ULP. IEEE 754 требует, чтобы результаты для большинства операций были правильными с точностью до 0,5 ULP, но ошибки могут накапливаться, что означает, что результат может быть не в пределах одного ULP от точного результата. Также есть пределы точности, поэтому в зависимости от количества цифр в результирующих значениях, вы также можете не работать со значениями одинаковых величин. Трансцендентные функции также несколько печально известны введением ошибки в вычисления.

Однако, если вы используете GNU glibc, sqrt будет корректным с точностью до 0,5 ULP (округлено), так что ваш конкретный пример будет работать (пренебрегая NaN, +/-0, +/-Inf). Хотя, вероятно, лучше определить некоторый эпсилон в качестве допустимого отклонения и использовать его в качестве границы. Например,

bool gt(double a, double b, double eps) {

  return (a > b - eps);
}

В зависимости от уровня точности, который вам нужен в расчетах, вы также можете использовать long double вместо.

Итак, чтобы ответить на ваши вопросы...

а) Гарантируется ли y*y >= x с помощью компилятора GCC?

Если вы используете встроенные функции GNU glibc или SSE2, да.

б) Будет ли это работать для других операций, таких как + - * / или даже std::cos() и std::acos()?

Предполагая, что вы используете GNU glibc и одну операцию, да. Хотя некоторые трансцендентальные не гарантированно правильно округлены.

в) Есть ли лучшие способы получить верхние / нижние границы?

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

Answer 2

Для GCC эта страница предполагает, что он будет работать, если вы используете встроенную функцию GCC sqrt __builtin_sqrt,

Кроме того, это поведение будет зависеть от того, как вы компилируете свой код и на какой машине он запущен.

1. Если процессор поддерживает SSE2, вы должны скомпилировать свой код с флагами -mfpmath=sse -msse2 чтобы убедиться, что все операции с плавающей запятой выполняются с использованием регистров SSE.

2. Если процессор не поддерживает SSE2, вы должны использовать long double введите значения с плавающей запятой и скомпилируйте с флагом -ffloat-store вынудить GCC не использовать регистры для хранения значений с плавающей запятой (за это вам понизят производительность)

Answer 3

Что касается

в) Есть ли лучшие способы получить верхние / нижние границы?

Другой способ - использовать другой режим округления, т.е. FE_UPWARD или же FE_DOWNWARD вместо по умолчанию FE_TONEAREST, См. /questions/34665357/izmenit-rezhim-okrugleniya-s-plavayuschej-zapyatoj/34665377#34665377 Это может быть медленнее, но лучше верхней / нижней границы.