3976
2007-02-06 13:14:10 +0000
libm rounding modes do not work correctly for many archs
2016-02-26 15:07:58 +0000
1
1
1
Unclassified
glibc
math
unspecified
All
All
RESOLVED
FIXED
P2
normal
---
1
madcoder
unassigned
bagnara
c_keil
debian-glibc
glibc-bugs
JoshuaHopp
vincent-srcware
oldest_to_newest
14882
0
madcoder
2007-02-06 13:14:10 +0000
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.
14921
1
madcoder
2007-02-06 22:13:19 +0000
(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
21833
2
vincent-srcware
2008-02-25 13:11:22 +0000
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.
25502
3
c_keil
2008-10-27 08:17:04 +0000
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.
25503
4
jakub
2008-10-27 08:25:06 +0000
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.
25504
5
vincent-srcware
2008-10-27 09:33:27 +0000
(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.
25506
6
joseph
2008-10-27 12:42:01 +0000
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.)
42130
7
vincent-srcware
2010-03-21 22:24:35 +0000
*** Bug 6869 has been marked as a duplicate of this bug. ***
53632
8
jsm28
2012-02-29 01:14:47 +0000
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).
53634
9
6256
vincent-srcware
2012-02-29 02:35:24 +0000
Created attachment 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.
53635
10
vincent-srcware
2012-02-29 02:54:10 +0000
(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
53730
11
jsm28
2012-03-02 15:16:39 +0000
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.
53732
12
jsm28
2012-03-02 16:04:12 +0000
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
53734
13
jsm28
2012-03-02 16:08:45 +0000
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
53739
14
jsm28
2012-03-02 20:54:41 +0000
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).
53775
15
jsm28
2012-03-05 12:27:36 +0000
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.
53776
16
vincent-srcware
2012-03-05 13:32:24 +0000
(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...
59417
17
aj
2013-01-10 19:22:58 +0000
*** Bug 15010 has been marked as a duplicate of this bug. ***
66767
18
jackie.rosen
2014-02-16 19:37:05 +0000
*** 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.
6256
2012-02-29 02:35:24 +0000
2012-02-29 02:35:24 +0000
Test a math function in the 4 rounding modes.
tfct-4rm.c
text/x-csrc
834
vincent-srcware
I2lmbmRlZiBGQ1QKIyBlcnJvciAiRkNUIG11c3QgYmUgdGhlIHRlc3RlZCBmdW5jdGlvbiIKI2Vu
ZGlmCgojaW5jbHVkZSA8c3RkaW8uaD4KI2luY2x1ZGUgPHN0ZGxpYi5oPgojaW5jbHVkZSA8bWF0
aC5oPgojaW5jbHVkZSA8ZmVudi5oPgoKI2RlZmluZSBTVFJJTkdJRlkoUykgI1MKI2RlZmluZSBN
QUtFX1NUUihTKSBTVFJJTkdJRlkoUykKCnN0YXRpYyBpbnQgcm5kWzRdID0geyBGRV9UT05FQVJF
U1QsIEZFX1RPV0FSRFpFUk8sIEZFX0RPV05XQVJELCBGRV9VUFdBUkQgfTsKc3RhdGljIGNoYXIg
cmNbNF0gPSAiTlpEVSI7CgppbnQgbWFpbiAoaW50IGFyZ2MsIGNoYXIgKmFyZ3ZbXSkKewogIGlu
dCBpOwoKICBmb3IgKGkgPSAxOyBpIDwgYXJnYzsgaSsrKQogICAgewogICAgICBpbnQgcjsKICAg
ICAgZG91YmxlIHg7CiAgICAgIGNoYXIgKmVuZDsKI2lmIFRXT0FSR1MKICAgICAgZG91YmxlIHk7
CiNlbmRpZgoKICAgICAgeCA9IHN0cnRvZCAoYXJndltpXSwgJmVuZCk7CiAgICAgIGlmICgqZW5k
ICE9ICdcMCcpCiAgICAgICAgZXhpdCAoRVhJVF9GQUlMVVJFKTsKCiNpZiBUV09BUkdTCiAgICAg
IGlmICgrK2kgPj0gYXJnYykKICAgICAgICBleGl0IChFWElUX0ZBSUxVUkUpOwogICAgICB5ID0g
c3RydG9kIChhcmd2W2ldLCAmZW5kKTsKICAgICAgaWYgKCplbmQgIT0gJ1wwJykKICAgICAgICBl
eGl0IChFWElUX0ZBSUxVUkUpOwojZW5kaWYKCiAgICAgIGZvciAociA9IDA7IHIgPCA0OyByKysp
CiAgICAgICAgewogICAgICAgICAgZG91YmxlIHY7CgogICAgICAgICAgZmVzZXRyb3VuZCAocm5k
W3JdKTsKI2lmIFRXT0FSR1MKICAgICAgICAgIHYgPSBGQ1QgKHgsIHkpOwogICAgICAgICAgcHJp
bnRmICgiJWM6ICIgTUFLRV9TVFIoRkNUKSAiKCUuMTdnLCUuMTdnKSA9ICUuMTdnXG4iLAogICAg
ICAgICAgICAgICAgICByY1tyXSwgeCwgeSwgdik7CiNlbHNlCiAgICAgICAgICB2ID0gRkNUICh4
KTsKICAgICAgICAgIHByaW50ZiAoIiVjOiAiIE1BS0VfU1RSKEZDVCkgIiglLjE3ZykgPSAlLjE3
Z1xuIiwgcmNbcl0sIHgsIHYpOwojZW5kaWYKICAgICAgICB9CiAgICB9CgogIHJldHVybiAwOwp9
CgovKiBFeGFtcGxlcyBvbiBhIERlYmlhbiBHTlUvTGludXggeDg2XzY0IG1hY2hpbmUgd2l0aCBn
bGliYyAyLjEzOgoKZ2NjIC1zdGQ9Yzk5IHRmY3QtNHJtLmMgLW8gdGZjdC00cm0gLWxtIC1ERkNU
PWV4cAouL3RmY3QtNHJtIDE1Ck46IGV4cCgxNSkgPSAzMjY5MDE3LjM3MjQ3MjExMDcKWjogZXhw
KDE1KSA9IDIwOTg2OTcuNjUzMTk0MDUzNgpEOiBleHAoMTUpID0gMjA5ODY5Ny42NTMxOTQwNTM2
ClU6IGV4cCgxNSkgPSAzMjY5MDE3LjM3MjQ3MjExMDcKCmdjYyAtc3RkPWM5OSB0ZmN0LTRybS5j
IC1vIHRmY3QtNHJtIC1sbSAtREZDVD1zaW4KLi90ZmN0LTRybSAxMDAKTjogc2luKDEwMCkgPSAt
MC41MDYzNjU2NDExMDk3NTg3OQpaOiBzaW4oMTAwKSA9IC02LjUxMzAwNTgyNDY2Nzg0NjYKRDog
c2luKDEwMCkgPSAtNi41MTMwMDU4MjQ2Njc4NDY2ClU6IHNpbigxMDApID0gLTAuNTA2MzY1NjQx
MTA5NzU4OQoKZ2NjIC1zdGQ9Yzk5IHRmY3QtNHJtLmMgLW8gdGZjdC00cm0gLWxtIC1ERkNUPWNv
cwpOOiBjb3MoMTAwKSA9IDAuODYyMzE4ODcyMjg3NjgzODkKWjogY29zKDEwMCkgPSAwLjMxMDg4
NTE1NDczODU3NTk4CkQ6IGNvcygxMDApID0gMC4zMTA4ODUxNTQ3Mzg1NzU5OApVOiBjb3MoMTAw
KSA9IDAuODYyMzE4ODcyMjg3Njg0CgpnY2MgLXN0ZD1jOTkgdGZjdC00cm0uYyAtbyB0ZmN0LTRy
bSAtbG0gLURGQ1Q9cG93IC1EVFdPQVJHUwouL3RmY3QtNHJtIDEuMDEgMS4xCk46IHBvdygxLjAx
LDEuMTAwMDAwMDAwMDAwMDAwMSkgPSAxLjAxMTAwNTQ4MzU3NzkyMzQKWjogcG93KDEuMDEsMS4x
MDAwMDAwMDAwMDAwMDAxKSA9IDEuMDExMDA1NDgzNTc3OTIzMgpEOiBwb3coMS4wMSwxLjEwMDAw
MDAwMDAwMDAwMDEpID0gMS4wMTEwMDU0ODM1Nzc5MjMyCnpzaDogc2VnbWVudGF0aW9uIGZhdWx0
IChjb3JlIGR1bXBlZCkgIC4vdGZjdC00cm0gMS4wMSAxLjEKCiovCg==