Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
A fix has been committed, but there's still a problem on loongson with libm updated: [...] This is definitely a problem in gdtoa. The values are always computed the same, but do not get printed correctly, as shown by this modified test program. Both lines should print 0.606531 for the long doubles, but don't always; this smells like an uninitialized value somewhere. Well, it turned out this was caused by the program stack being aligned to an 8 byte boundary, instead of a 16 byte boundary, causing long double variables on the stack to be misaligned 50% of the time. This has been fixed now. Miod
Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
A fix has been committed, but there's still a problem on loongson with
libm updated:
$ ls -l /usr/lib/libm.so.*
-r--r--r-- 1 root bin 926033 Feb 12 12:17 /usr/lib/libm.so.9.0
$ cc -o expl expl.c -O2 -pipe -lm
$ for in in 1 2 3 4 5 6 ; do ./expl ; done
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
theta == 1, pr == 0.606531, pr2 == 0.606531
theta == 1, pr == 0.606531, pr2 == 0.606531
theta == 1, pr == 0.606531, pr2 == 0.606531
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
theta == 1, pr == 0.606531, pr2 == 0.606531
This is definitely a problem in gdtoa. The values are always computed
the same, but do not get printed correctly, as shown by this modified
test program. Both lines should print 0.606531 for the long doubles, but
don't always; this smells like an uninitialized value somewhere.
#include stdio.h
#include stdlib.h
#include math.h
int
main(void)
{
double theta = 1;
long double lambda, pr, pr2;
static const uint64_t zpr[2] = {
#ifdef __MIPSEL__
0xa000ULL, 0x3ffe368b2fc6f960ULL
#else
0x3ffe368b2fc6f960ULL, 0xa000ULL
#endif
}, zpr2[2] = {
#ifdef __MIPSEL__
0x9fe7aceb46aa619cULL, 0x3ffe368b2fc6f960ULL
#else
0x3ffe368b2fc6f960ULL, 0x9fe7aceb46aa619cULL
#endif
};
lambda = (0.5*theta);
pr = exp(-lambda);
pr2 = expl(-lambda);
printf(theta == %g, pr == %Lg, pr2 == %Lg\n, theta, pr, pr2);
printf(theta == %g, pr == %Lg, pr2 == %Lg\n, theta, zpr[0], zpr[1],
zpr2[0], zpr2[1]);
exit(0);
}
Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
On Thu, 6 Feb 2014, Donovan Watteau wrote:
David Coppa wrote:
Take the following reduced test-case, adapted from what R's code
does:
---8---
#include stdio.h
#include stdlib.h
#include math.h
int main(void) {
double theta = 1;
long double lambda, pr, pr2;
lambda = (0.5*theta);
pr = exp(-lambda);
pr2 = expl(-lambda);
printf(theta == %g, pr == %Lg, pr2 == %Lg\n, theta, pr, pr2);
exit(0);
}
---8---
This produces the following output on Linux (x86_64):
theta == 1, pr == 0.606531, pr2 == 0.606531
While on OpenBSD -current amd64:
theta == 1, pr == 0.606531, pr2 == nan
FWIW, it looks even stranger on loongson:
$ cc -o expl expl.c -O2 -pipe -lm
$ ./expl
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
$ ./expl
theta == 1, pr == 0.606531, pr2 == 0.606531
$ ./expl
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
$ sysctl kern.version
kern.version=OpenBSD 5.5-beta (GENERIC) #106: Mon Feb 3 01:47:15 MST 2014
t...@loongson.openbsd.org:/usr/src/sys/arch/loongson/compile/GENERIC
A fix has been committed, but there's still a problem on loongson with
libm updated:
$ ls -l /usr/lib/libm.so.*
-r--r--r-- 1 root bin 926033 Feb 12 12:17 /usr/lib/libm.so.9.0
$ cc -o expl expl.c -O2 -pipe -lm
$ for in in 1 2 3 4 5 6 ; do ./expl ; done
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
theta == 1, pr == 0.606531, pr2 == 0.606531
theta == 1, pr == 0.606531, pr2 == 0.606531
theta == 1, pr == 0.606531, pr2 == 0.606531
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
theta == 1, pr == 0.606531, pr2 == 0.606531
Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
On 2/12/14, Donovan Watteau tso...@gmail.com wrote:
On Thu, 6 Feb 2014, Donovan Watteau wrote:
David Coppa wrote:
Take the following reduced test-case, adapted from what R's code
does:
---8---
#include stdio.h
#include stdlib.h
#include math.h
int main(void) {
double theta = 1;
long double lambda, pr, pr2;
lambda = (0.5*theta);
pr = exp(-lambda);
pr2 = expl(-lambda);
printf(theta == %g, pr == %Lg, pr2 == %Lg\n, theta, pr, pr2);
exit(0);
}
---8---
This produces the following output on Linux (x86_64):
theta == 1, pr == 0.606531, pr2 == 0.606531
While on OpenBSD -current amd64:
theta == 1, pr == 0.606531, pr2 == nan
FWIW, it looks even stranger on loongson:
$ cc -o expl expl.c -O2 -pipe -lm
$ ./expl
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
$ ./expl
theta == 1, pr == 0.606531, pr2 == 0.606531
$ ./expl
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
$ sysctl kern.version
kern.version=OpenBSD 5.5-beta (GENERIC) #106: Mon Feb 3 01:47:15 MST
2014
t...@loongson.openbsd.org:/usr/src/sys/arch/loongson/compile/GENERIC
A fix has been committed, but there's still a problem on loongson with
libm updated:
$ ls -l /usr/lib/libm.so.*
-r--r--r-- 1 root bin 926033 Feb 12 12:17 /usr/lib/libm.so.9.0
$ cc -o expl expl.c -O2 -pipe -lm
$ for in in 1 2 3 4 5 6 ; do ./expl ; done
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
theta == 1, pr == 0.606531, pr2 == 0.606531
theta == 1, pr == 0.606531, pr2 == 0.606531
theta == 1, pr == 0.606531, pr2 == 0.606531
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
theta == 1, pr == 0.606531, pr2 == 0.606531
This isn't related to exp/expl. Looks like a bug in either compiler
or libc/printf/gdtoa. I don't have the hardware so I couldn't tell
much more.
As usual, there are a few ways to fix it:
1. Debug it and provide a diff,
2. Donate hardware. Oh, BTW I have a request for loongson in
http://www.openbsd.org/want.html sitting for over a year...
Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
Date: Mon, 10 Feb 2014 22:35:02 -0700 (MST)
From: Martynas Venckus marty...@cvs.openbsd.org
Here's a diff that sticks a bit closer to the original code. It's
equivalent to your diff, and admittedly purely a matter of taste which
version to prefer.
I prefer my version better. It's not '93 anymore and compilers are
able to convert 0.0L and -1.0L precisely, otherwise we have a huge
problem. There's no need to obfuscate here by manually converting
to floating point representation.
Here's a diff which also includes same fix for ld128. OK?
Fair enough. Please commit.
Index: ld128/s_floorl.c
===
RCS file: /cvs/src/lib/libm/src/ld128/s_floorl.c,v
retrieving revision 1.1
diff -u -r1.1 s_floorl.c
--- ld128/s_floorl.c 6 Jul 2011 00:02:42 - 1.1
+++ ld128/s_floorl.c 11 Feb 2014 05:24:15 -
@@ -34,10 +34,11 @@
jj0 = ((i048)0x7fff)-0x3fff;
if(jj048) {
if(jj00) { /* raise inexact if x != 0 */
- if(huge+x0.0) {/* return 0*sign(x) if |x|1 */
- if(i0=0) {i0=i1=0;}
+ if(huge+x0.0) {
+ if(i0=0)
+ return 0.0L;
else if(((i00x7fffLL)|i1)!=0)
- { i0=0xbfffULL;i1=0;}
+ return -1.0L;
}
} else {
i = (0xULL)jj0;
Index: ld80/s_floorl.c
===
RCS file: /cvs/src/lib/libm/src/ld80/s_floorl.c,v
retrieving revision 1.2
diff -u -r1.2 s_floorl.c
--- ld80/s_floorl.c 25 Jul 2011 16:20:09 - 1.2
+++ ld80/s_floorl.c 11 Feb 2014 05:24:01 -
@@ -35,10 +35,11 @@
jj0 = (se0x7fff)-0x3fff;
if(jj031) {
if(jj00) { /* raise inexact if x != 0 */
- if(huge+x0.0) {/* return 0*sign(x) if |x|1 */
- if(sx==0) {se=0;i0=i1=0;}
+ if(huge+x0.0) {
+ if(sx==0)
+ return 0.0L;
else if(((se0x7fff)|i0|i1)!=0)
- { se=0xbfff;i0=i1=0;}
+ return -1.0L;
}
} else {
i = (0x7fff)jj0;
Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
On Fri, Feb 7, 2014 at 4:03 PM, David Coppa dco...@gmail.com wrote:
On Fri, Feb 7, 2014 at 3:49 PM, Mark Kettenis mark.kette...@xs4all.nl wrote:
Date: Thu, 6 Feb 2014 23:07:58 -0800
From: Martynas Venckus marty...@venck.us
Yup.Does this diff fix it for you?
Here's a diff that sticks a bit closer to the original code. It's
equivalent to your diff, and admittedly purely a matter of taste which
version to prefer.
Index: s_floorl.c
===
RCS file: /home/cvs/src/lib/libm/src/ld80/s_floorl.c,v
retrieving revision 1.2
diff -u -p -r1.2 s_floorl.c
--- s_floorl.c 25 Jul 2011 16:20:09 - 1.2
+++ s_floorl.c 7 Feb 2014 14:43:19 -
@@ -38,7 +38,7 @@ floorl(long double x)
if(huge+x0.0) {/* return 0*sign(x) if |x|1 */
if(sx==0) {se=0;i0=i1=0;}
else if(((se0x7fff)|i0|i1)!=0)
- { se=0xbfff;i0=i1=0;}
+ { se=0xbfff;i0=0x8000;i1=0;}
}
} else {
i = (0x7fff)jj0;
Just tested, and this one works too, of course...
I'd be happy if this or the other one could be committed.
Ping
Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
Here's a diff that sticks a bit closer to the original code. It's
equivalent to your diff, and admittedly purely a matter of taste which
version to prefer.
I prefer my version better. It's not '93 anymore and compilers are
able to convert 0.0L and -1.0L precisely, otherwise we have a huge
problem. There's no need to obfuscate here by manually converting
to floating point representation.
Here's a diff which also includes same fix for ld128. OK?
Index: ld128/s_floorl.c
===
RCS file: /cvs/src/lib/libm/src/ld128/s_floorl.c,v
retrieving revision 1.1
diff -u -r1.1 s_floorl.c
--- ld128/s_floorl.c6 Jul 2011 00:02:42 - 1.1
+++ ld128/s_floorl.c11 Feb 2014 05:24:15 -
@@ -34,10 +34,11 @@
jj0 = ((i048)0x7fff)-0x3fff;
if(jj048) {
if(jj00) { /* raise inexact if x != 0 */
- if(huge+x0.0) {/* return 0*sign(x) if |x|1 */
- if(i0=0) {i0=i1=0;}
+ if(huge+x0.0) {
+ if(i0=0)
+ return 0.0L;
else if(((i00x7fffLL)|i1)!=0)
- { i0=0xbfffULL;i1=0;}
+ return -1.0L;
}
} else {
i = (0xULL)jj0;
Index: ld80/s_floorl.c
===
RCS file: /cvs/src/lib/libm/src/ld80/s_floorl.c,v
retrieving revision 1.2
diff -u -r1.2 s_floorl.c
--- ld80/s_floorl.c 25 Jul 2011 16:20:09 - 1.2
+++ ld80/s_floorl.c 11 Feb 2014 05:24:01 -
@@ -35,10 +35,11 @@
jj0 = (se0x7fff)-0x3fff;
if(jj031) {
if(jj00) { /* raise inexact if x != 0 */
- if(huge+x0.0) {/* return 0*sign(x) if |x|1 */
- if(sx==0) {se=0;i0=i1=0;}
+ if(huge+x0.0) {
+ if(sx==0)
+ return 0.0L;
else if(((se0x7fff)|i0|i1)!=0)
- { se=0xbfff;i0=i1=0;}
+ return -1.0L;
}
} else {
i = (0x7fff)jj0;
Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
On Tue, Feb 11, 2014 at 6:35 AM, Martynas Venckus
marty...@cvs.openbsd.org wrote:
Here's a diff that sticks a bit closer to the original code. It's
equivalent to your diff, and admittedly purely a matter of taste which
version to prefer.
I prefer my version better. It's not '93 anymore and compilers are
able to convert 0.0L and -1.0L precisely, otherwise we have a huge
problem. There's no need to obfuscate here by manually converting
to floating point representation.
Here's a diff which also includes same fix for ld128. OK?
OK with me.
Ciao,
David
Index: ld128/s_floorl.c
===
RCS file: /cvs/src/lib/libm/src/ld128/s_floorl.c,v
retrieving revision 1.1
diff -u -r1.1 s_floorl.c
--- ld128/s_floorl.c6 Jul 2011 00:02:42 - 1.1
+++ ld128/s_floorl.c11 Feb 2014 05:24:15 -
@@ -34,10 +34,11 @@
jj0 = ((i048)0x7fff)-0x3fff;
if(jj048) {
if(jj00) { /* raise inexact if x != 0 */
- if(huge+x0.0) {/* return 0*sign(x) if |x|1 */
- if(i0=0) {i0=i1=0;}
+ if(huge+x0.0) {
+ if(i0=0)
+ return 0.0L;
else if(((i00x7fffLL)|i1)!=0)
- { i0=0xbfffULL;i1=0;}
+ return -1.0L;
}
} else {
i = (0xULL)jj0;
Index: ld80/s_floorl.c
===
RCS file: /cvs/src/lib/libm/src/ld80/s_floorl.c,v
retrieving revision 1.2
diff -u -r1.2 s_floorl.c
--- ld80/s_floorl.c 25 Jul 2011 16:20:09 - 1.2
+++ ld80/s_floorl.c 11 Feb 2014 05:24:01 -
@@ -35,10 +35,11 @@
jj0 = (se0x7fff)-0x3fff;
if(jj031) {
if(jj00) { /* raise inexact if x != 0 */
- if(huge+x0.0) {/* return 0*sign(x) if |x|1 */
- if(sx==0) {se=0;i0=i1=0;}
+ if(huge+x0.0) {
+ if(sx==0)
+ return 0.0L;
else if(((se0x7fff)|i0|i1)!=0)
- { se=0xbfff;i0=i1=0;}
+ return -1.0L;
}
} else {
i = (0x7fff)jj0;
Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
On Fri, Feb 7, 2014 at 8:07 AM, Martynas Venckus marty...@venck.us wrote:
Yup.Does this diff fix it for you?
Yeah! It works.
And R-3's testsuite is also happy now.
Thanks a lot! And thanks to Daniel too...
ok dcoppa@ to commit it, obviously
Ciao,
David
On 2/6/14, Daniel Dickman didick...@gmail.com wrote:
I think I recently ran into a similar issue but I suspect the root cause
might be the same. I think the floorl function is wrong for numbers slightly
larger than -1 to numbers slightly below 0. In this range floorl returns -0
instead of -1.
On Feb 5, 2014, at 3:57 AM, David Coppa dco...@openbsd.org wrote:
Hi!
I hit this problem while working on updating math/R from version
2.15.3 to the latest version (3.0.2).
It started happening since upstream switched from double functions
to C99 long double functions (expl, fabsl, ...), during the R-3
development cycle.
Take the following reduced test-case, adapted from what R's code
does:
---8---
#include stdio.h
#include stdlib.h
#include math.h
int main(void) {
double theta = 1;
long double lambda, pr, pr2;
lambda = (0.5*theta);
pr = exp(-lambda);
pr2 = expl(-lambda);
printf(theta == %g, pr == %Lg, pr2 == %Lg\n, theta, pr, pr2);
exit(0);
}
---8---
This produces the following output on Linux (x86_64):
theta == 1, pr == 0.606531, pr2 == 0.606531
While on OpenBSD -current amd64:
theta == 1, pr == 0.606531, pr2 == nan
And indeed R-3's testsuite fails with the error message
NaNs produced:
Warning in pchisq(1e-300, df = 0, ncp = lam) : NaNs produced
stopifnot(all.equal(p00, exp(-lam/2)),
+ all.equal(p.0, exp(-lam/2)))
Error: all.equal(p.0, exp(-lam/2)) is not TRUE
Execution halted
Is this a bug in our expl() ?
Ciao,
David
Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
Date: Thu, 6 Feb 2014 23:07:58 -0800
From: Martynas Venckus marty...@venck.us
Yup.Does this diff fix it for you?
Here's a diff that sticks a bit closer to the original code. It's
equivalent to your diff, and admittedly purely a matter of taste which
version to prefer.
Index: s_floorl.c
===
RCS file: /home/cvs/src/lib/libm/src/ld80/s_floorl.c,v
retrieving revision 1.2
diff -u -p -r1.2 s_floorl.c
--- s_floorl.c 25 Jul 2011 16:20:09 - 1.2
+++ s_floorl.c 7 Feb 2014 14:43:19 -
@@ -38,7 +38,7 @@ floorl(long double x)
if(huge+x0.0) {/* return 0*sign(x) if |x|1 */
if(sx==0) {se=0;i0=i1=0;}
else if(((se0x7fff)|i0|i1)!=0)
- { se=0xbfff;i0=i1=0;}
+ { se=0xbfff;i0=0x8000;i1=0;}
}
} else {
i = (0x7fff)jj0;
Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
On Fri, Feb 7, 2014 at 3:49 PM, Mark Kettenis mark.kette...@xs4all.nl wrote:
Date: Thu, 6 Feb 2014 23:07:58 -0800
From: Martynas Venckus marty...@venck.us
Yup.Does this diff fix it for you?
Here's a diff that sticks a bit closer to the original code. It's
equivalent to your diff, and admittedly purely a matter of taste which
version to prefer.
Index: s_floorl.c
===
RCS file: /home/cvs/src/lib/libm/src/ld80/s_floorl.c,v
retrieving revision 1.2
diff -u -p -r1.2 s_floorl.c
--- s_floorl.c 25 Jul 2011 16:20:09 - 1.2
+++ s_floorl.c 7 Feb 2014 14:43:19 -
@@ -38,7 +38,7 @@ floorl(long double x)
if(huge+x0.0) {/* return 0*sign(x) if |x|1 */
if(sx==0) {se=0;i0=i1=0;}
else if(((se0x7fff)|i0|i1)!=0)
- { se=0xbfff;i0=i1=0;}
+ { se=0xbfff;i0=0x8000;i1=0;}
}
} else {
i = (0x7fff)jj0;
Just tested, and this one works too, of course...
I'd be happy if this or the other one could be committed.
ciao,
David
Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
David Coppa wrote:
Take the following reduced test-case, adapted from what R's code
does:
---8---
#include stdio.h
#include stdlib.h
#include math.h
int main(void) {
double theta = 1;
long double lambda, pr, pr2;
lambda = (0.5*theta);
pr = exp(-lambda);
pr2 = expl(-lambda);
printf(theta == %g, pr == %Lg, pr2 == %Lg\n, theta, pr, pr2);
exit(0);
}
---8---
This produces the following output on Linux (x86_64):
theta == 1, pr == 0.606531, pr2 == 0.606531
While on OpenBSD -current amd64:
theta == 1, pr == 0.606531, pr2 == nan
FWIW, it looks even stranger on loongson:
$ cc -o expl expl.c -O2 -pipe -lm
$ ./expl
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
$ ./expl
theta == 1, pr == 0.606531, pr2 == 0.606531
$ ./expl
theta == 1, pr == -9.15569e-2474, pr2 == 6.10667e-4944
$ sysctl kern.version
kern.version=OpenBSD 5.5-beta (GENERIC) #106: Mon Feb 3 01:47:15 MST 2014
t...@loongson.openbsd.org:/usr/src/sys/arch/loongson/compile/GENERIC
Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
I think I recently ran into a similar issue but I suspect the root cause might
be the same. I think the floorl function is wrong for numbers slightly larger
than -1 to numbers slightly below 0. In this range floorl returns -0 instead of
-1.
On Feb 5, 2014, at 3:57 AM, David Coppa dco...@openbsd.org wrote:
Hi!
I hit this problem while working on updating math/R from version
2.15.3 to the latest version (3.0.2).
It started happening since upstream switched from double functions
to C99 long double functions (expl, fabsl, ...), during the R-3
development cycle.
Take the following reduced test-case, adapted from what R's code
does:
---8---
#include stdio.h
#include stdlib.h
#include math.h
int main(void) {
double theta = 1;
long double lambda, pr, pr2;
lambda = (0.5*theta);
pr = exp(-lambda);
pr2 = expl(-lambda);
printf(theta == %g, pr == %Lg, pr2 == %Lg\n, theta, pr, pr2);
exit(0);
}
---8---
This produces the following output on Linux (x86_64):
theta == 1, pr == 0.606531, pr2 == 0.606531
While on OpenBSD -current amd64:
theta == 1, pr == 0.606531, pr2 == nan
And indeed R-3's testsuite fails with the error message
NaNs produced:
Warning in pchisq(1e-300, df = 0, ncp = lam) : NaNs produced
stopifnot(all.equal(p00, exp(-lam/2)),
+ all.equal(p.0, exp(-lam/2)))
Error: all.equal(p.0, exp(-lam/2)) is not TRUE
Execution halted
Is this a bug in our expl() ?
Ciao,
David
Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
Yup.Does this diff fix it for you?
On 2/6/14, Daniel Dickman didick...@gmail.com wrote:
I think I recently ran into a similar issue but I suspect the root cause
might be the same. I think the floorl function is wrong for numbers slightly
larger than -1 to numbers slightly below 0. In this range floorl returns -0
instead of -1.
On Feb 5, 2014, at 3:57 AM, David Coppa dco...@openbsd.org wrote:
Hi!
I hit this problem while working on updating math/R from version
2.15.3 to the latest version (3.0.2).
It started happening since upstream switched from double functions
to C99 long double functions (expl, fabsl, ...), during the R-3
development cycle.
Take the following reduced test-case, adapted from what R's code
does:
---8---
#include stdio.h
#include stdlib.h
#include math.h
int main(void) {
double theta = 1;
long double lambda, pr, pr2;
lambda = (0.5*theta);
pr = exp(-lambda);
pr2 = expl(-lambda);
printf(theta == %g, pr == %Lg, pr2 == %Lg\n, theta, pr, pr2);
exit(0);
}
---8---
This produces the following output on Linux (x86_64):
theta == 1, pr == 0.606531, pr2 == 0.606531
While on OpenBSD -current amd64:
theta == 1, pr == 0.606531, pr2 == nan
And indeed R-3's testsuite fails with the error message
NaNs produced:
Warning in pchisq(1e-300, df = 0, ncp = lam) : NaNs produced
stopifnot(all.equal(p00, exp(-lam/2)),
+ all.equal(p.0, exp(-lam/2)))
Error: all.equal(p.0, exp(-lam/2)) is not TRUE
Execution halted
Is this a bug in our expl() ?
Ciao,
David
Index: s_floorl.c
===
RCS file: /cvs/src/lib/libm/src/ld80/s_floorl.c,v
retrieving revision 1.2
diff -u -r1.2 s_floorl.c
--- s_floorl.c 25 Jul 2011 16:20:09 - 1.2
+++ s_floorl.c 7 Feb 2014 07:01:59 -
@@ -35,10 +35,12 @@
jj0 = (se0x7fff)-0x3fff;
if(jj031) {
if(jj00) { /* raise inexact if x != 0 */
- if(huge+x0.0) {/* return 0*sign(x) if |x|1 */
- if(sx==0) {se=0;i0=i1=0;}
- else if(((se0x7fff)|i0|i1)!=0)
- { se=0xbfff;i0=i1=0;}
+ if(huge+x0.0) {
+ if(sx==0) {
+return 0.0L;
+} else if(((se0x7fff)|i0|i1)!=0) {
+ return -1.0L;
+}
}
} else {
i = (0x7fff)jj0;
exp() / expl() on Linux and OpenBSD (expl() bug?)
Hi!
I hit this problem while working on updating math/R from version
2.15.3 to the latest version (3.0.2).
It started happening since upstream switched from double functions
to C99 long double functions (expl, fabsl, ...), during the R-3
development cycle.
Take the following reduced test-case, adapted from what R's code
does:
---8---
#include stdio.h
#include stdlib.h
#include math.h
int main(void) {
double theta = 1;
long double lambda, pr, pr2;
lambda = (0.5*theta);
pr = exp(-lambda);
pr2 = expl(-lambda);
printf(theta == %g, pr == %Lg, pr2 == %Lg\n, theta, pr, pr2);
exit(0);
}
---8---
This produces the following output on Linux (x86_64):
theta == 1, pr == 0.606531, pr2 == 0.606531
While on OpenBSD -current amd64:
theta == 1, pr == 0.606531, pr2 == nan
And indeed R-3's testsuite fails with the error message
NaNs produced:
Warning in pchisq(1e-300, df = 0, ncp = lam) : NaNs produced
stopifnot(all.equal(p00, exp(-lam/2)),
+ all.equal(p.0, exp(-lam/2)))
Error: all.equal(p.0, exp(-lam/2)) is not TRUE
Execution halted
Is this a bug in our expl() ?
Ciao,
David
Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
Date: Wed, 5 Feb 2014 01:57:33 -0700
From: David Coppa dco...@openbsd.org
Hi!
I hit this problem while working on updating math/R from version
2.15.3 to the latest version (3.0.2).
It started happening since upstream switched from double functions
to C99 long double functions (expl, fabsl, ...), during the R-3
development cycle.
Take the following reduced test-case, adapted from what R's code
does:
---8---
#include stdio.h
#include stdlib.h
#include math.h
int main(void) {
double theta = 1;
long double lambda, pr, pr2;
lambda = (0.5*theta);
pr = exp(-lambda);
pr2 = expl(-lambda);
printf(theta == %g, pr == %Lg, pr2 == %Lg\n, theta, pr, pr2);
exit(0);
}
---8---
This produces the following output on Linux (x86_64):
theta == 1, pr == 0.606531, pr2 == 0.606531
While on OpenBSD -current amd64:
theta == 1, pr == 0.606531, pr2 == nan
And indeed R-3's testsuite fails with the error message
NaNs produced:
Warning in pchisq(1e-300, df = 0, ncp = lam) : NaNs produced
stopifnot(all.equal(p00, exp(-lam/2)),
+ all.equal(p.0, exp(-lam/2)))
Error: all.equal(p.0, exp(-lam/2)) is not TRUE
Execution halted
Is this a bug in our expl() ?
Yes.
Re: exp() / expl() on Linux and OpenBSD (expl() bug?)
Date: Wed, 5 Feb 2014 01:57:33 -0700
From: David Coppa dco...@openbsd.org
Hi!
I hit this problem while working on updating math/R from version
2.15.3 to the latest version (3.0.2).
It started happening since upstream switched from double functions
to C99 long double functions (expl, fabsl, ...), during the R-3
development cycle.
Take the following reduced test-case, adapted from what R's code
does:
---8---
#include stdio.h
#include stdlib.h
#include math.h
int main(void) {
double theta = 1;
long double lambda, pr, pr2;
lambda = (0.5*theta);
pr = exp(-lambda);
pr2 = expl(-lambda);
printf(theta == %g, pr == %Lg, pr2 == %Lg\n, theta, pr, pr2);
exit(0);
}
---8---
This produces the following output on Linux (x86_64):
theta == 1, pr == 0.606531, pr2 == 0.606531
While on OpenBSD -current amd64:
theta == 1, pr == 0.606531, pr2 == nan
And indeed R-3's testsuite fails with the error message
NaNs produced:
Warning in pchisq(1e-300, df = 0, ncp = lam) : NaNs produced
stopifnot(all.equal(p00, exp(-lam/2)),
+ all.equal(p.0, exp(-lam/2)))
Error: all.equal(p.0, exp(-lam/2)) is not TRUE
Execution halted
Is this a bug in our expl() ?
Oh, btw, the quad-precision code used on sparc64 gets this right. So
the bug is probably somewhere in src/lib/limb/src/ld80.
