View Bug Activity | Format For Printing
It seems that changing the rounding modes make some functions (like exp, or sin, or ...) buggy on many archs. The testing code is: =============================================================================== #include <stdio.h> #include <stdlib.h> #include <math.h> #include <fenv.h> static int rnd[4] = { FE_TONEAREST, FE_TOWARDZERO, FE_DOWNWARD, FE_UPWARD }; static char rc[4] = "NZDU"; int main (int argc, char *argv[]) { int i; for (i = 1; i < argc; i++) { int r; double x; char *end; x = strtod (argv[i], &end); if (*end != '\0') exit (EXIT_FAILURE); for (r = 0; r < 4; r++) { double y; fesetround (rnd[r]); y = exp (x); printf ("%c: exp(%.17g) = %.17g\n", rc[r], x, y); } } return 0; } =============================================================================== I've checked on my amd64 that it indeed works in 32bits mode (and it also work on an i386 machine) but it does not in 64bits mode: [madcoder mad] gcc -m32 -lm -o a a.c ; ./a 1 2 N: exp(1) = 2.7182818284590451 Z: exp(1) = 2.7182818284590451 D: exp(1) = 2.7182818284590451 U: exp(1) = 2.7182818284590455 N: exp(2) = 7.3890560989306504 Z: exp(2) = 7.3890560989306495 D: exp(2) = 7.3890560989306495 U: exp(2) = 7.3890560989306504 [madcoder mad] gcc -m64 -lm -o a a.c ; ./a 1 2 N: exp(1) = 2.7182818284590451 Z: exp(1) = 2.7182818284590451 D: exp(1) = 2.7182818284590451 U: exp(1) = 0.04788398250919005 N: exp(2) = 7.3890560989306504 Z: exp(2) = 4.0037745305985499 D: exp(2) = 4.0037745305985499 U: exp(2) = 7.3890560989306504 I tested it on other archs, here is the summary: - i386, m68k, ia64 have been tested OK. - arm: only FE_ROUNDTONEAREST exists (which is ok as per C standard) - amd64, mips, ppc, hppa, s390, sparc produce wrong results, in their own unique way. I've not been able to test alpha yet though.
(In reply to comment #0) > I've not been able to test alpha yet though. alpha gives curious results as all values are exactly the same, which is quite reasonnable, but may hide a bug too. It gave (for 1 2 3 4 5): N: exp(1) = 2.7182818284590451 Z: exp(1) = 2.7182818284590451 D: exp(1) = 2.7182818284590451 U: exp(1) = 2.7182818284590451 N: exp(2) = 7.3890560989306504 Z: exp(2) = 7.3890560989306504 D: exp(2) = 7.3890560989306504 U: exp(2) = 7.3890560989306504 N: exp(3) = 20.085536923187668 Z: exp(3) = 20.085536923187668 D: exp(3) = 20.085536923187668 U: exp(3) = 20.085536923187668 N: exp(4) = 54.598150033144236 Z: exp(4) = 54.598150033144236 D: exp(4) = 54.598150033144236 U: exp(4) = 54.598150033144236 N: exp(5) = 148.4131591025766 Z: exp(5) = 148.4131591025766 D: exp(5) = 148.4131591025766 U: exp(5) = 148.4131591025766
You can see other results in the bug I originally reported on the Debian BTS: http://bugs.debian.org/cgi-bin/bugreport.cgi?bug=153022 In short, sin can even give out-of-range results, and pow can segfault.
At least on my Core2 Duo it's also not working with still other values (32bit version is ok): gcc -m64 -lm -o a a.c; ./a 1 2 N: exp(1) = 2.7182818284590451 Z: exp(1) = 2.7182818284590451 D: exp(1) = 2.7182818284590451 U: exp(1) = 7.1387612927397726 N: exp(2) = 7.3890560989306504 Z: exp(2) = 4.0037745305985499 D: exp(2) = 4.0037745305985499 U: exp(2) = 7.3890560989306504 As far as I dug into it, the 32bit and 64bit versions use other code. 64bit comes from the IBM Accurate Mathematical Library. For exp sysdeps/ieee754/dbl-64/e_exp.c produces the wrong results.
The functions have defined behavior only in the default rounding mode (round-to-even), anything else is undefined behavior and completely programmer's fault for calling the functions in those rounding modes.
(In reply to comment #4) > The functions have defined behavior only in the default rounding mode > (round-to-even), anything else is undefined behavior and completely programmer's > fault for calling the functions in those rounding modes. No, it isn't. There's no undefined behavior there. The C standard says (F.9): "Whether the functions honor the rounding direction mode is implementation-defined." So, this just means that an implementation doesn't need to return a result rounded to the correct direction (this is implementation-defined). Still it should return an approximate value whatever the rounding direction mode is, and not behave erratically.
Subject: Re: libm rounding modes do not work correctly for many archs On Mon, 27 Oct 2008, vincent+libc at vinc17 dot org wrote: > So, this just means that an implementation doesn't need to return a result > rounded to the correct direction (this is implementation-defined). Still it > should return an approximate value whatever the rounding direction mode is, and > not behave erratically. Furthermore, some functions whose implementations only work in round-to-nearest mode do save the mode then set round-to-nearest using feholdexcept and fesetround (e.g. sysdeps/ieee754/dbl-64/e_exp2.c). I think this is the best thing to do in such cases. (Regarding the issue of rounding modes in signal handlers, I think the best thing is for ABI documents to make explicit that library functions may temporarily save and restore rounding modes; that's what is being done in the Power Architecture ABI document being worked on, to document what the de facto ABI is right now.)
*** Bug 6869 has been marked as a duplicate of this bug. ***