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

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

[jsm28 at gcc dot gnu.org]

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

Joseph Myers <jsm28 at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
         AssignedTo|aj at suse dot de           |unassigned at sourceware
                   |                            |dot org

--- Comment #8 from Joseph Myers <jsm28 at gcc dot gnu.org> 2012-02-29 01:14:47 UTC ---
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).

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[vincent-srcware at vinc17 dot net]

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

--- Comment #9 from Vincent Lefèvre <vincent-srcware at vinc17 dot net> 2012-02-29 02:35:24 UTC ---
Created attachment 6256
  --> http://sourceware.org/bugzilla/attachment.cgi?id=6256
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.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[vincent-srcware at vinc17 dot net]

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

--- Comment #10 from Vincent Lefèvre <vincent-srcware at vinc17 dot net> 2012-02-29 02:54:10 UTC ---
(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

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[jsm28 at gcc dot gnu.org]

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

--- Comment #11 from Joseph Myers <jsm28 at gcc dot gnu.org> 2012-03-02 15:16:39 UTC ---
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.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[jsm28 at gcc dot gnu.org]

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

--- Comment #12 from Joseph Myers <jsm28 at gcc dot gnu.org> 2012-03-02 16:04:12 UTC ---
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

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[jsm28 at gcc dot gnu.org]

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

--- Comment #13 from Joseph Myers <jsm28 at gcc dot gnu.org> 2012-03-02 16:08:45 UTC ---
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

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[jsm28 at gcc dot gnu.org]

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

--- Comment #14 from Joseph Myers <jsm28 at gcc dot gnu.org> 2012-03-02 20:54:41 UTC ---
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).

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[jsm28 at gcc dot gnu.org]

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

Joseph Myers <jsm28 at gcc dot gnu.org> changed:

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

--- Comment #15 from Joseph Myers <jsm28 at gcc dot gnu.org> 2012-03-05 12:27:36 UTC ---
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.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[vincent-srcware at vinc17 dot net]

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

--- Comment #16 from Vincent Lefèvre <vincent-srcware at vinc17 dot net> 2012-03-05 13:32:24 UTC ---
(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...

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[aj at suse dot de]

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

Andreas Jaeger <aj at suse dot de> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |JoshuaHopp at web dot de

--- Comment #17 from Andreas Jaeger <aj at suse dot de> 2013-01-10 19:22:58 UTC ---
*** Bug 15010 has been marked as a duplicate of this bug. ***

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[jackie.rosen at hushmail dot com]

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

Jackie Rosen <jackie.rosen at hushmail dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |jackie.rosen at hushmail dot com

--- Comment #18 from Jackie Rosen <jackie.rosen at hushmail dot com> ---
*** 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.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

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

[schwab at sourceware dot org]

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

Andreas Schwab <schwab at sourceware dot org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|jackie.rosen at hushmail dot com   |

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[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.

[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  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

[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

[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

[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

[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

[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

[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.