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. ***
Problem with exp confirmed on x86_64 with current sources. Patch for exp proposed: http://sourceware.org/ml/libc-alpha/2012-02/msg00748.html Given that patch applied: * I cannot confirm the problem with sin or cos on x86_64 (though tests should be added to the testsuite). * pow (1.6, 1.6) does not segfault, but the result in round-upward mode is substantially inaccurate; pow will need a similar fix (and test in the testsuite). If other functions have problems in current sources, it would be best to have a separate bug for each function with problems (identifying clearly the platform on which the problem occurs).
Created attachment 6256 [details] Test a math function in the 4 rounding modes. (In reply to comment #8) > * I cannot confirm the problem with sin or cos on x86_64 (though tests should > be added to the testsuite). I still get the bug on the argument 100 under Debian (glibc 2.13). > * pow (1.6, 1.6) does not segfault, but the result in round-upward mode is > substantially inaccurate; I confirm, but pow(1.01,1.1) crashes: N: pow(1.01,1.1000000000000001) = 1.0110054835779234 Z: pow(1.01,1.1000000000000001) = 1.0110054835779232 D: pow(1.01,1.1000000000000001) = 1.0110054835779232 zsh: segmentation fault (core dumped) ./tfct-4rm 1.01 1.1 > pow will need a similar fix (and test in the testsuite). Yes, like the other functions. > If other functions have problems in current sources, [...] I would say that each function probably has the same problem. I did the tests with gcc -std=c99 tfct-4rm.c -o tfct-4rm -lm -DFCT=exp gcc -std=c99 tfct-4rm.c -o tfct-4rm -lm -DFCT=sin gcc -std=c99 tfct-4rm.c -o tfct-4rm -lm -DFCT=cos gcc -std=c99 tfct-4rm.c -o tfct-4rm -lm -DFCT=pow -DTWOARGS using the attached code.
(In reply to comment #9) > I would say that each function probably has the same problem. Actually log and atan seem to behave correctly (and even honor the rounding mode). But tan, sinh and cosh are buggy: $ ./tfct-4rm 1.6 N: tan(1.6000000000000001) = -34.232532735557314 Z: tan(1.6000000000000001) = -34.232532735557307 D: tan(1.6000000000000001) = -34.232532735557307 U: tan(1.6000000000000001) = -1.9000553624671392 $ ./tfct-4rm 100 N: sinh(100) = 1.3440585709080678e+43 Z: sinh(100) = 1.3440585709080676e+43 D: sinh(100) = 1.3440585709080676e+43 U: sinh(100) = -5.7576991416613074e+42 $ ./tfct-4rm 100 N: cosh(100) = 1.3440585709080678e+43 Z: cosh(100) = 1.3440585709080676e+43 D: cosh(100) = 1.3440585709080676e+43 U: cosh(100) = -5.7576991416613074e+42 tanh seems accurate, but it doesn't honor the rounding mode: $ ./tfct-4rm 0.8 N: tanh(0.80000000000000004) = 0.66403677026784913 Z: tanh(0.80000000000000004) = 0.6640367702678488 D: tanh(0.80000000000000004) = 0.66403677026784902 U: tanh(0.80000000000000004) = 0.66403677026784891
exp fixed by: commit 28afd92dbdb4fef4358051aad5cb944a9527a4b5 Author: Joseph Myers <joseph@codesourcery.com> Date: Fri Mar 2 15:12:53 2012 +0000 Fix exp in non-default rounding modes (bug 3976). This may fix some cases of sinh and cosh (where they use exp internally) though tests still need adding for those functions and other functions still need to be fixed.
sinh and cosh appear to have been fixed by the exp fix. I confirm what you report for tan and pow. For sin and cos, I don't see the problem for 100 with current sources, but do for 7: N: cos(7) = 0.7539022543433046 Z: cos(7) = 0.7539022543433046 D: cos(7) = 0.7539022543433046 U: cos(7) = -3.5593280928702544e+244 N: sin(7) = 0.65698659871878906 Z: sin(7) = 0.65698659871878906 D: sin(7) = 0.65698659871878906 U: sin(7) = 6.5993918533387462e+246
Examples for pow involving arguments that are exactly representable in all floating-point types (always nice for the testsuite): N: pow(1.0625,1.125) = 1.0705822930287614 Z: pow(1.0625,1.125) = 1.0705822930287612 D: pow(1.0625,1.125) = 1.0705822930287612 U: pow(1.0625,1.125) = 0.032633984947502387 N: pow(1.5,1.03125) = 1.5191270987147432 Z: pow(1.5,1.03125) = 1.0003045181943615 D: pow(1.5,1.03125) = 1.0003045181943615 U: pow(1.5,1.03125) = 1.5191270987147434
sin, cos, tan fixed by: commit 804360ed837e3347c9cd9738f25345d2587a1242 Author: Joseph Myers <joseph@codesourcery.com> Date: Fri Mar 2 20:51:39 2012 +0000 Fix sin, cos, tan in non-default rounding modes (bug 3976).
cosh, sinh tests (some errors in round-upward mode up to 27ulp) added by: commit ca811b2256d2e48c7288219e9e11dcbab3000f19 Author: Joseph Myers <joseph@codesourcery.com> Date: Mon Mar 5 12:20:24 2012 +0000 Test cosh, sinh in non-default rounding modes (bug 3976). pow fixed by: commit b7cd39e8f8c5cf2844f20eb03f545d19c4c25987 Author: Joseph Myers <joseph@codesourcery.com> Date: Mon Mar 5 12:22:46 2012 +0000 Fix pow in non-default rounding modes (bug 3976). So all the reported issues are fixed here. If cases are found of functions that still (in any rounding mode, for any of float, double, long double) produce wild results or crash with glibc after my patches - or of functions that should always produce correctly rounded results such as rint, nextafter, fma, sqrt, that fail to do so (in any rounding mode, for any of float, double, long double), please file them as separate bugs for each function found to have a problem.
(In reply to comment #10) > (In reply to comment #9) > > I would say that each function probably has the same problem. > > Actually log and atan seem to behave correctly (and even honor the rounding > mode). I forgot that since Ziv's strategy is used to provide correct rounding, my tests where incomplete (I tested only on a few random values). To understand the problem, here's what Ziv's strategy does: 1. Compute an approximation with slightly more precision than the target precision, typically by representing it as an exact sum of two floating-point values a+b, where |b| is much smaller than |a|. An error bound e1 is determined and a rounding test is performed. Most of the time, the correct rounding can be guaranteed by this test, so that the computed result can be returned. But there's a low probability of failure (problem known as the "Table Maker's Dilemma"), in which case a second step is necessary. 2. Ditto with even more precision and a lower error bound e2. In case of failure, a third step is necessary. 3. Computation with even more precision. IIRC, a 768-bit precision is used here in IBM's library (that's the mp*.c files), which assumes that this is sufficient. So, for functions that are believed to behave correctly in the directed rounding modes, and for which the active rounding mode is not changed temporarily in the function, one should consider the following remarks: One would need to test cases that fall in (2) and in (3). The probabilities of failure depend on the algorithm and implied error bounds e1 and e2, and possibly on the magnitude of the input (or some other criteria). To get cases in (2) from random tests, one may need to test several dozens of thousands of inputs (analyzing the code can give a better idea of what is needed -- it doesn't seem to be much documented). Ideally, tools like gcov to get code coverage would be useful. From my work on the Table Maker's Dilemma, I also have a list of the hardest-to-round cases ("worst cases") for some functions in some domains, which should probably fall in (3). AFAIK, in IBM's library, (3) uses integer arithmetic, so that it shouldn't depend on the active rounding mode (though I'm not completely sure because FP may also be used, e.g. at the end). Note that testing cases for (3) only is not sufficient, because with such cases, the results from (2) are discarded. It is also possible that in directed rounding modes, (1) and/or (2) always fail, meaning that a slower step would be performed. This wouldn't be incorrect, but rather inefficient. This should be checked too...
*** Bug 15010 has been marked as a duplicate of this bug. ***
*** Bug 260998 has been marked as a duplicate of this bug. *** Seen from the domain http://volichat.com Page where seen: http://volichat.com/adult-chat-rooms Marked for reference. Resolved as fixed @bugzilla.