Re: exp2(3) bug?
You might want to read up on floating point arithmetic. (rounding and
representation)
On 06/02/14 05:13, Daniel Dickman wrote:
I hit this problem while working with the numpy 1.8.1 regress suite
which has some tests that are currently failing.
Here is a reduced test case of the logaddexp2 python function which
ends up calling exp2. Is this a bug in the openbsd exp2
implementation?
---8---
#include stdio.h
#include stdlib.h
#include math.h
int main(void) {
double x;
double y;
// x = log2(5)
x = 2.32192809489;
// y = 2**(log2(5))
y = exp2(x);
printf(expected: 5.0\n);
printf(actual: %f\n, y);
---8---
on a linux/x86_64 machine:
# gcc -lm test.c ./a.out
expected: 5.0
actual: 5.00
on an openbsd/i386 machine:
# gcc -lm test.c ./a.out
expected: 5.0
actual: 4.994404
Re: exp2(3) bug?
Date: Mon, 02 Jun 2014 09:34:20 +0200
From: Benjamin Baier program...@netzbasis.de
You might want to read up on floating point arithmetic. (rounding and
representation)
Well, the difference between 4.994404 and 5.0 is a bit large to blame
rounding and binary representation. And other OpenBSD platforms
(amd64, sparc64, powerpc) return the expected result. So I'd say that
there is a bug in the i386-specific implementation of exp2(3).
On 06/02/14 05:13, Daniel Dickman wrote:
I hit this problem while working with the numpy 1.8.1 regress suite
which has some tests that are currently failing.
Here is a reduced test case of the logaddexp2 python function which
ends up calling exp2. Is this a bug in the openbsd exp2
implementation?
---8---
#include stdio.h
#include stdlib.h
#include math.h
int main(void) {
double x;
double y;
// x = log2(5)
x = 2.32192809489;
// y = 2**(log2(5))
y = exp2(x);
printf(expected: 5.0\n);
printf(actual: %f\n, y);
---8---
on a linux/x86_64 machine:
# gcc -lm test.c ./a.out
expected: 5.0
actual: 5.00
on an openbsd/i386 machine:
# gcc -lm test.c ./a.out
expected: 5.0
actual: 4.994404
Re: exp2(3) bug?
Ok i might have been quick to judge, because i don't really trust
floating point arithmetic.
btw. exp2f works correct on i386.
On 06/02/14 10:17, Mark Kettenis wrote:
Date: Mon, 02 Jun 2014 09:34:20 +0200
From: Benjamin Baier program...@netzbasis.de
You might want to read up on floating point arithmetic. (rounding and
representation)
Well, the difference between 4.994404 and 5.0 is a bit large to blame
rounding and binary representation. And other OpenBSD platforms
(amd64, sparc64, powerpc) return the expected result. So I'd say that
there is a bug in the i386-specific implementation of exp2(3).
On 06/02/14 05:13, Daniel Dickman wrote:
I hit this problem while working with the numpy 1.8.1 regress suite
which has some tests that are currently failing.
Here is a reduced test case of the logaddexp2 python function which
ends up calling exp2. Is this a bug in the openbsd exp2
implementation?
---8---
#include stdio.h
#include stdlib.h
#include math.h
int main(void) {
double x;
double y;
// x = log2(5)
x = 2.32192809489;
// y = 2**(log2(5))
y = exp2(x);
printf(expected: 5.0\n);
printf(actual: %f\n, y);
---8---
on a linux/x86_64 machine:
# gcc -lm test.c ./a.out
expected: 5.0
actual: 5.00
on an openbsd/i386 machine:
# gcc -lm test.c ./a.out
expected: 5.0
actual: 4.994404
Re: exp2(3) bug?
Date: Mon, 2 Jun 2014 10:17:53 +0200 (CEST)
From: Mark Kettenis mark.kette...@xs4all.nl
Date: Mon, 02 Jun 2014 09:34:20 +0200
From: Benjamin Baier program...@netzbasis.de
You might want to read up on floating point arithmetic. (rounding and
representation)
Well, the difference between 4.994404 and 5.0 is a bit large to blame
rounding and binary representation. And other OpenBSD platforms
(amd64, sparc64, powerpc) return the expected result. So I'd say that
there is a bug in the i386-specific implementation of exp2(3).
And here is a fix. There actually isn't any i386-specific code, but
i386 is special and needs STRICT_ALIGN() to work properly for double
as well as float. FreeBSD made the same change a while ago:
http://svnweb.FreeBSD.org/base/head/lib/msun/src/math_private.h?revision=240827view=markup
Haven't run the regression tests yet with this change.
Index: src/math_private.h
===
RCS file: /cvs/src/lib/libm/src/math_private.h,v
retrieving revision 1.16
diff -u -p -r1.16 math_private.h
--- src/math_private.h 12 Nov 2013 20:35:09 - 1.16
+++ src/math_private.h 2 Jun 2014 09:30:13 -
@@ -349,7 +349,7 @@ do {
\
#defineSTRICT_ASSIGN(type, lval, rval) do {\
volatile type __lval; \
\
- if (sizeof(type) = sizeof(double)) \
+ if (sizeof(type) = sizeof(long double))\
(lval) = (rval);\
else { \
__lval = (rval);\
Re: exp2(3) bug?
Hi Benjamin,
On Mon, Jun 2, 2014 at 3:34 AM, Benjamin Baier program...@netzbasis.de wrote:
You might want to read up on floating point arithmetic. (rounding and
representation)
The error in floating point calculations is typically measured in
ulps. The error shown in this example is far too big for the
existing calculation to be correct. I'm linking to the wikipedia
article for some background reading for you but there's a lot more
info on the web on this topic:
http://en.wikipedia.org/wiki/Unit_in_the_last_place
On 06/02/14 05:13, Daniel Dickman wrote:
I hit this problem while working with the numpy 1.8.1 regress suite
which has some tests that are currently failing.
Here is a reduced test case of the logaddexp2 python function which
ends up calling exp2. Is this a bug in the openbsd exp2
implementation?
---8---
#include stdio.h
#include stdlib.h
#include math.h
int main(void) {
double x;
double y;
// x = log2(5)
x = 2.32192809489;
// y = 2**(log2(5))
y = exp2(x);
printf(expected: 5.0\n);
printf(actual: %f\n, y);
---8---
on a linux/x86_64 machine:
# gcc -lm test.c ./a.out
expected: 5.0
actual: 5.00
on an openbsd/i386 machine:
# gcc -lm test.c ./a.out
expected: 5.0
actual: 4.994404
Re: exp2(3) bug?
On Mon, Jun 2, 2014 at 5:33 AM, Mark Kettenis mark.kette...@xs4all.nl wrote:
Date: Mon, 2 Jun 2014 10:17:53 +0200 (CEST)
From: Mark Kettenis mark.kette...@xs4all.nl
Date: Mon, 02 Jun 2014 09:34:20 +0200
From: Benjamin Baier program...@netzbasis.de
You might want to read up on floating point arithmetic. (rounding and
representation)
Well, the difference between 4.994404 and 5.0 is a bit large to blame
rounding and binary representation. And other OpenBSD platforms
(amd64, sparc64, powerpc) return the expected result. So I'd say that
there is a bug in the i386-specific implementation of exp2(3).
And here is a fix. There actually isn't any i386-specific code, but
i386 is special and needs STRICT_ALIGN() to work properly for double
as well as float. FreeBSD made the same change a while ago:
Fixes my test case and the numpy test. So ok daniel@ for what it's worth...
http://svnweb.FreeBSD.org/base/head/lib/msun/src/math_private.h?revision=240827view=markup
Haven't run the regression tests yet with this change.
Index: src/math_private.h
===
RCS file: /cvs/src/lib/libm/src/math_private.h,v
retrieving revision 1.16
diff -u -p -r1.16 math_private.h
--- src/math_private.h 12 Nov 2013 20:35:09 - 1.16
+++ src/math_private.h 2 Jun 2014 09:30:13 -
@@ -349,7 +349,7 @@ do {
\
#defineSTRICT_ASSIGN(type, lval, rval) do {\
volatile type __lval; \
\
- if (sizeof(type) = sizeof(double)) \
+ if (sizeof(type) = sizeof(long double))\
(lval) = (rval);\
else { \
__lval = (rval);\
Re: exp2(3) bug?
On 6/2/14, Mark Kettenis mark.kette...@xs4all.nl wrote:
Date: Mon, 2 Jun 2014 10:17:53 +0200 (CEST)
From: Mark Kettenis mark.kette...@xs4all.nl
Date: Mon, 02 Jun 2014 09:34:20 +0200
From: Benjamin Baier program...@netzbasis.de
You might want to read up on floating point arithmetic. (rounding and
representation)
Well, the difference between 4.994404 and 5.0 is a bit large to blame
rounding and binary representation. And other OpenBSD platforms
(amd64, sparc64, powerpc) return the expected result. So I'd say that
there is a bug in the i386-specific implementation of exp2(3).
And here is a fix. There actually isn't any i386-specific code, but
i386 is special and needs STRICT_ALIGN() to work properly for double
as well as float. FreeBSD made the same change a while ago:
http://svnweb.FreeBSD.org/base/head/lib/msun/src/math_private.h?revision=240827view=markup
You are quick. Thanks for figuring this out.
Ugh ... it was intentionally left broken as an optimization. That's crazy.
OK martynas@.
Haven't run the regression tests yet with this change.
Index: src/math_private.h
===
RCS file: /cvs/src/lib/libm/src/math_private.h,v
retrieving revision 1.16
diff -u -p -r1.16 math_private.h
--- src/math_private.h12 Nov 2013 20:35:09 - 1.16
+++ src/math_private.h2 Jun 2014 09:30:13 -
@@ -349,7 +349,7 @@ do {
\
#define STRICT_ASSIGN(type, lval, rval) do {\
volatile type __lval; \
\
- if (sizeof(type) = sizeof(double)) \
+ if (sizeof(type) = sizeof(long double))\
(lval) = (rval);\
else { \
__lval = (rval);\
