c++ sqrt guaranteed precision, upper/lower bound

Question

I have to check an inequality containing square roots. To avoid incorrect results due to floating point inaccuracy and rounding, I use std::nextafter() to get an upper/lower bound:

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

a) Is y*y >= x guaranteed using GCC compiler?

b) Will this work for other operations like + - * / or even std::cos() and std::acos()?

c) Are there better ways to get upper/lower bounds?

Update: I read this is not guaranteed by the C++ Standard, but should work according to IEEE-754. Will this work with the GCC compiler?

Answer 1

In general, floating point operations will result in some ULP error. IEEE 754 requires that results for most operations be correct to within 0.5 ULP, but errors can accumulate, which means a result may not be within one ULP of the the exact result. There are limits to precision as well, so depending on the number of digits there are in resulting values, you also may not be working with values of the same magnitudes. Transcendental functions are also somewhat notorious for introducing error into calculations.

However, if you're using GNU glibc, sqrt will be correct to within 0.5 ULP (rounded), so you're specific example would work (neglecting NaN, +/-0, +/-Inf). Although, it's probably better to define some epsilon as your error tolerance and use that as your bound. For exmaple,

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

  return (a > b - eps);
}

Depending on the level of precision you need in calculations, you also may want to use long double instead.

So, to answer your questions...

a) Is y*y >= x guaranteed using GCC compiler?

Assuming you use GNU glibc or SSE2 intrinsics, yes.

b) Will this work for other operations like + - * / or even std::cos() and std::acos()?

Assuming you use GNU glibc and one operation, yes. Although some transcendentals are not guaranteed correctly rounded.

c) Are there better ways to get upper/lower bounds?

You need to know what your error tolerance in calculations is, and use that as an epsilon (which may be larger than one ULP).

Answer 2

For GCC this page suggests that it will work if you use the GCC builtin sqrt function __builtin_sqrt.

Additionally this behavior will be dependent on how you compile your code and the machine that it is run on

1. If the processor supports SSE2 then you should compile your code with the flags -mfpmath=sse -msse2 to ensure that all floating point operations are done using the SSE registers.

2. If the processor doesn't support SSE2 then you should use the long double type for the floating point values and compile with the flag -ffloat-store to force GCC to not use registers to store floating point values (you'll have a performance penalty for doing this)