[Bug math/3976] New: libm rounding modes do not work correctly for many archs

[Bug math/3976] New: libm rounding modes do not work correctly for many archs

[madcoder at debian dot org]

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.

-- 
           Summary: libm rounding modes do not work correctly for many archs
           Product: glibc
           Version: unspecified
            Status: NEW
          Severity: normal
          Priority: P2
         Component: math
        AssignedTo: aj at suse dot de
        ReportedBy: madcoder at debian dot org
                CC: glibc-bugs at sources dot redhat dot com


http://sourceware.org/bugzilla/show_bug.cgi?id=3976

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/3976] libm rounding modes do not work correctly for many archs

[madcoder at debian dot org]

------- Additional Comments From madcoder at debian dot org  2007-02-06 22:13 -------
(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


-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=3976

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/3976] libm rounding modes do not work correctly for many archs

[debian-glibc at lists dot debian dot org]

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |debian-glibc at lists dot
                   |                            |debian dot org


http://sourceware.org/bugzilla/show_bug.cgi?id=3976

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/3976] libm rounding modes do not work correctly for many archs

[vincent+libc at vinc17 dot org]

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |vincent+libc at vinc17 dot
                   |                            |org


http://sourceware.org/bugzilla/show_bug.cgi?id=3976

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/3976] libm rounding modes do not work correctly for many archs

[vincent+libc at vinc17 dot org]

------- Additional Comments From vincent+libc at vinc17 dot org  2008-02-25 13:11 -------
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.

-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=3976

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/3976] libm rounding modes do not work correctly for many archs

[c_keil at yahoo dot de]

------- Additional Comments From c_keil at yahoo dot de  2008-10-27 08:17 -------
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.

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |c_keil at yahoo dot de


http://sourceware.org/bugzilla/show_bug.cgi?id=3976

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/3976] libm rounding modes do not work correctly for many archs

[jakub at redhat dot com]

------- Additional Comments From jakub at redhat dot com  2008-10-27 08:25 -------
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.

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |RESOLVED
         Resolution|                            |INVALID


http://sourceware.org/bugzilla/show_bug.cgi?id=3976

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/3976] libm rounding modes do not work correctly for many archs

[vincent+libc at vinc17 dot org]

------- Additional Comments From vincent+libc at vinc17 dot org  2008-10-27 09:33 -------
(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.

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|RESOLVED                    |REOPENED
         Resolution|INVALID                     |


http://sourceware.org/bugzilla/show_bug.cgi?id=3976

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/3976] libm rounding modes do not work correctly for many archs

[joseph at codesourcery dot com]

------- Additional Comments From joseph at codesourcery dot com  2008-10-27 12:42 -------
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.)



-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=3976

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/3976] libm rounding modes do not work correctly for many archs

[vincent+libc at vinc17 dot org]

------- Additional Comments From vincent+libc at vinc17 dot org  2010-03-21 22:24 -------
*** Bug 6869 has been marked as a duplicate of this bug. ***

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |bagnara at cs dot unipr dot
                   |                            |it


http://sourceware.org/bugzilla/show_bug.cgi?id=3976

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.