Summary:  libm rounding modes do not work correctly for many archs  

Product:  glibc  Reporter:  Pierre Habouzit <madcoder> 
Component:  math  Assignee:  Not yet assigned to anyone <unassigned> 
Status:  RESOLVED FIXED  
Severity:  normal  CC:  bagnara, c_keil, debianglibc, glibcbugs, JoshuaHopp, vincentsrcware 
Priority:  P2  Flags:  fweimer:
security

Version:  unspecified  
Target Milestone:    
Host:  Target:  
Build:  Last reconfirmed:  
Attachments:  Test a math function in the 4 rounding modes. 
Description
Pierre Habouzit
20070206 13:14:10 UTC
(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/cgibin/bugreport.cgi?bug=153022 In short, sin can even give outofrange 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/dbl64/e_exp.c produces the wrong results. The functions have defined behavior only in the default rounding mode (roundtoeven), 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 > (roundtoeven), 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 implementationdefined." So, this just means that an implementation doesn't need to return a result rounded to the correct direction (this is implementationdefined). 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 implementationdefined). 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
roundtonearest mode do save the mode then set roundtonearest using
feholdexcept and fesetround (e.g. sysdeps/ieee754/dbl64/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/libcalpha/201202/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 roundupward 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 roundupward 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) ./tfct4rm 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 tfct4rm.c o tfct4rm lm DFCT=exp gcc std=c99 tfct4rm.c o tfct4rm lm DFCT=sin gcc std=c99 tfct4rm.c o tfct4rm lm DFCT=cos gcc std=c99 tfct4rm.c o tfct4rm 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: $ ./tfct4rm 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 $ ./tfct4rm 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 $ ./tfct4rm 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: $ ./tfct4rm 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 nondefault 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 floatingpoint 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 nondefault rounding modes (bug 3976). cosh, sinh tests (some errors in roundupward 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 nondefault 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 nondefault 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 floatingpoint 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 768bit 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 hardesttoround 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/adultchatrooms Marked for reference. Resolved as fixed @bugzilla. 