Is there a bug in controlled rounding using `exp`?

user_13036124 0
Viewed 117

I'm observing incorrect (IMO) rounding behaviour on some platforms as follows:

Calculate the value of log(2) under rounding modes to FE_DOWNWARD and FE_UPWARD (see <fenv.h>). In all cases I've seen, the rounded downward result is less than the rounded upward result, which is what I would expect (result is inexact).

Now using the same rounding modes, call exp with each result. Since exp is monotonically increasing (slope is ~2 in the area of interest), I would expect the these results to be even further apart. But on some platforms tested, the two results from exp are equivalent.

Here's a simple test program:

#include <stdio.h>
#include <math.h>
#include <fenv.h>

int main(void) {
  fesetround(FE_DOWNWARD);
  double xlo = log(2);
  double xxlo = exp(xlo);
  printf("xlo=%.16e xxlo=%.16e\n",xlo,xxlo);

  fesetround(FE_UPWARD);
  double xhi = log(2);
  double xxhi = exp(xhi);
  printf("xhi=%.16e xxhi=%.16e\n",xhi,xxhi); 

  printf("Delta x=%.16e\n",xhi-xlo);
  if (xxlo == xxhi) 
    printf("FAIL\n");
  else 
    printf("Delta xx=%.16e\n",xxhi-xxlo);

  return 0;
}

One (good) result on Xeon E5520 @ 2.27GHz:

xlo=6.9314718055994528e-01 xxlo=1.9999999999999997e+00
xhi=6.9314718055994540e-01 xxhi=2.0000000000000005e+00
Delta x=1.1102230246251566e-16
Delta xx=6.6613381477509393e-16

Other processors that have passed the test: AMD Ryzen 9 3950X, Intel i7-5557U, Intel i7-3770, Intel m7-6Y75.

But a couple have reportedly failed (xxlo==xxhi): Intel E5-2650 and AMD EPYC 7401, coincidentally both server class machines. I don't have direct access to those but running the test at https://repl.it/languages/c also fails:

clang version 7.0.0-3~ubuntu0.18.04.1 (tags/RELEASE_700/final)
 clang-7 -pthread -lm -o main main.c
 ./main
xlo=6.93147180559945286e-01 xxlo=2.00000000000000000e+00
xhi=6.93147180559945398e-01 xxhi=2.00000000000000000e+00
Delta x=1.11022302462515655e-16
FAIL

I'm not so much interested in the actual precision of the results, but more in the relative values. So am I missing something, or is this just wrong? And if it's wrong, is it hardware rounding error or issues with the software math libraries, or even the compiler/settings? I'm also interested in collecting any results from running the simple test to see if any additional conclusions can be drawn from that data.

ieee-754 floating-point rounding c
3 Answers

The spread (round-down to round-up) gets magnified by the derivative of the function. So, if your original bound is 1ULP, I would expect exp(x) at 2.0 to have a spread of 3. That's (original spread) * (derivative) + 1. The +1 is because you are doing another round-down and round-up.

If the derivative is 2.5, I would expect the output spread to be between 3 and 4.

Please provide the hex for the numbers; this should be definitive by avoiding the bin<->dec conversions (which seem to be flawed on some of the machines, too).

None of what I have said invalidates your conclusion that the servers have inconsistent hardware, hence flawed IEEE implementation.

Oh, one more thing; check to see if the optimizer worked too hard -- by eliminating the exp after seeing the log. Or maybe using an FMA, which sometimes gives you better rounding errors.

user_17668310

@Rick James

I agree with your diagnosis. Here's the output from my original program with additional print statements providing the hex data. I'm not much of a C programmer, so if this isn't what you requested, be more explicit.

xlo=6.9314718055994528e-01 xxlo=2.0000000000000000e+00
xlo=0x1.62e42fefa39efp-1 xxlo=0x1p+1
xhi=6.9314718055994540e-01 xxhi=2.0000000000000000e+00
xhi=0x1.62e42fefa39fp-1 xxhi=0x1p+1
Delta x=1.1102230246251566e-16
Delta x=0x1p-53
FAIL
user_130361240

The OP asks for

My objective is to produce a pair of floating point values guaranteed to contain the answer.

All C implementations are lazy enough to just (barely) satisfy the C standard, as mentioned in the comments, and this is enough to satisfy the vast majority of users, and with a very good performance given the difficulties to comply to all platforms, toolchains, etc.

If you require upper and lower bounds to some mathematical univariate function, then there are very good interval arithmetic libraries such as GNU MPRIA, and most recently the Arb Library, which provide rigourous enclosing bounds for the evaluation of many mathematical functions with very good independency of platform (mainly because they rely on the GNU Multiple Precision Arithmetic Library GMP).

user_172281820
Related