The C99 standard documents how the <math.h> functions should behave in
certain corner cases. These are doucemnted in Appendix F, for those
with access to standard, and supposedly based on IEC 60559 which is
more widely known as IEEE 754. Our implementation of pow(3) gets some
of them wrong, most notably that:
- 1**y == 1, even if y is NaN
- (-1)**+-Inf == 1
Interestly enough the cephes testsuite that's in regress/lib/libm has
checks for these that are wrong. Noted this because the
quad-precision implementation did this right causing some tests to
fail ;).
The diff below fixes the corner cases for the single-, double- and
80-bit extended precision and adjusts the tests.
ok?
Index: lib/libm/src/e_pow.c
===================================================================
RCS file: /home/cvs/src/lib/libm/src/e_pow.c,v
retrieving revision 1.11
diff -u -p -r1.11 e_pow.c
--- lib/libm/src/e_pow.c 3 Jul 2013 04:46:36 -0000 1.11
+++ lib/libm/src/e_pow.c 1 Aug 2013 12:46:02 -0000
@@ -109,6 +109,9 @@ pow(double x, double y)
/* y==zero: x**0 = 1 */
if((iy|ly)==0) return one;
+ /* x==1: 1**y = 1, even if y is NaN */
+ if (hx==0x3ff00000 && lx == 0) return one;
+
/* +-NaN return x+y */
if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
@@ -138,7 +141,7 @@ pow(double x, double y)
if(ly==0) {
if (iy==0x7ff00000) { /* y is +-inf */
if(((ix-0x3ff00000)|lx)==0)
- return y - y; /* inf**+-1 is NaN */
+ return one; /* (-1)**+-inf is NaN */
else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
return (hy>=0)? y: zero;
else /* (|x|<1)**-,+inf = inf,0 */
Index: lib/libm/src/e_powf.c
===================================================================
RCS file: /home/cvs/src/lib/libm/src/e_powf.c,v
retrieving revision 1.6
diff -u -p -r1.6 e_powf.c
--- lib/libm/src/e_powf.c 27 Oct 2009 23:59:29 -0000 1.6
+++ lib/libm/src/e_powf.c 1 Aug 2013 12:46:02 -0000
@@ -64,6 +64,9 @@ powf(float x, float y)
/* y==zero: x**0 = 1 */
if(iy==0) return one;
+ /* x==1: 1**y = 1, even if y is NaN */
+ if (hx==0x3f800000) return one;
+
/* +-NaN return x+y */
if(ix > 0x7f800000 ||
iy > 0x7f800000)
@@ -87,7 +90,7 @@ powf(float x, float y)
/* special value of y */
if (iy==0x7f800000) { /* y is +-inf */
if (ix==0x3f800000)
- return y - y; /* inf**+-1 is NaN */
+ return one; /* (-1)**+-inf is NaN */
else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */
return (hy>=0)? y: zero;
else /* (|x|<1)**-,+inf = inf,0 */
Index: lib/libm/src/ld80/e_powl.c
===================================================================
RCS file: /home/cvs/src/lib/libm/src/ld80/e_powl.c,v
retrieving revision 1.3
diff -u -p -r1.3 e_powl.c
--- lib/libm/src/ld80/e_powl.c 20 Jul 2011 21:02:51 -0000 1.3
+++ lib/libm/src/ld80/e_powl.c 1 Aug 2013 12:54:29 -0000
@@ -211,6 +211,9 @@ long e;
if( y == 0.0L )
return( 1.0L );
+if( x == 1.0L )
+ return( 1.0L );
+
if( isnan(x) )
return( x );
if( isnan(y) )
@@ -219,10 +222,7 @@ if( isnan(y) )
if( y == 1.0L )
return( x );
-if( !isfinite(y) && (x == -1.0L || x == 1.0L) )
- return y - y; /* +-1**inf is NaN */
-
-if( x == 1.0L )
+if( !isfinite(y) && x == -1.0L )
return( 1.0L );
if( y >= LDBL_MAX )
Index: regress/lib/libm/cephes/testvect.c
===================================================================
RCS file: /home/cvs/src/regress/lib/libm/cephes/testvect.c,v
retrieving revision 1.1
diff -u -p -r1.1 testvect.c
--- regress/lib/libm/cephes/testvect.c 30 May 2011 20:23:35 -0000 1.1
+++ regress/lib/libm/cephes/testvect.c 1 Aug 2013 12:46:02 -0000
@@ -332,12 +332,12 @@ static struct twoarguments test2[] =
{"pow", pow, &MINF, &MTHREE, &MZERO, 0},
{"pow", pow, &MINF, &MTWO, &ZERO, 0},
{"pow", pow, &NAN, &ONE, &NAN, 0},
- {"pow", pow, &ONE, &NAN, &NAN, 0},
+ {"pow", pow, &ONE, &NAN, &ONE, 0},
{"pow", pow, &NAN, &NAN, &NAN, 0},
- {"pow", pow, &ONE, &INF, &NAN, 0},
- {"pow", pow, &MONE, &INF, &NAN, 0},
- {"pow", pow, &ONE, &MINF, &NAN, 0},
- {"pow", pow, &MONE, &MINF, &NAN, 0},
+ {"pow", pow, &ONE, &INF, &ONE, 0},
+ {"pow", pow, &MONE, &INF, &ONE, 0},
+ {"pow", pow, &ONE, &MINF, &ONE, 0},
+ {"pow", pow, &MONE, &MINF, &ONE, 0},
{"pow", pow, &MTWO, &HALF, &NAN, 0},
{"pow", pow, &ZERO, &MTHREE, &INF, 0},
{"pow", pow, &MZERO, &MTHREE, &MINF, 0},
Index: regress/lib/libm/cephes/testvectl.c
===================================================================
RCS file: /home/cvs/src/regress/lib/libm/cephes/testvectl.c,v
retrieving revision 1.2
diff -u -p -r1.2 testvectl.c
--- regress/lib/libm/cephes/testvectl.c 8 Jul 2011 16:49:05 -0000 1.2
+++ regress/lib/libm/cephes/testvectl.c 1 Aug 2013 12:46:02 -0000
@@ -325,12 +325,12 @@ static struct twoarguments test2[] =
{"powl", powl, &MINFL, &MTHREEL, &NEGZEROL, 0},
{"powl", powl, &MINFL, &MTWOL, &ZEROL, 0},
{"powl", powl, &NANL, &ONEL, &NANL, 0},
- {"powl", powl, &ONEL, &NANL, &NANL, 0},
+ {"powl", powl, &ONEL, &NANL, &ONEL, 0},
{"powl", powl, &NANL, &NANL, &NANL, 0},
- {"powl", powl, &ONEL, &INFINITYL, &NANL, 0},
- {"powl", powl, &MONEL, &INFINITYL, &NANL, 0},
- {"powl", powl, &ONEL, &MINFL, &NANL, 0},
- {"powl", powl, &MONEL, &MINFL, &NANL, 0},
+ {"powl", powl, &ONEL, &INFINITYL, &ONEL, 0},
+ {"powl", powl, &MONEL, &INFINITYL, &ONEL, 0},
+ {"powl", powl, &ONEL, &MINFL, &ONEL, 0},
+ {"powl", powl, &MONEL, &MINFL, &ONEL, 0},
{"powl", powl, &MTWOL, &HALFL, &NANL, 0},
{"powl", powl, &ZEROL, &MTHREEL, &INFINITYL, 0},
{"powl", powl, &NEGZEROL, &MTHREEL, &MINFL, 0},
Index: regress/lib/libm/cephes/testvectll.c
===================================================================
RCS file: /home/cvs/src/regress/lib/libm/cephes/testvectll.c,v
retrieving revision 1.2
diff -u -p -r1.2 testvectll.c
--- regress/lib/libm/cephes/testvectll.c 8 Jul 2011 16:49:05 -0000
1.2
+++ regress/lib/libm/cephes/testvectll.c 1 Aug 2013 12:46:02 -0000
@@ -328,12 +328,12 @@ static struct twoarguments test2[] =
{"powl", powl, &MINFL, &MTHREEL, &NEGZEROL, 0},
{"powl", powl, &MINFL, &MTWOL, &ZEROL, 0},
{"powl", powl, &NANL, &ONEL, &NANL, 0},
- {"powl", powl, &ONEL, &NANL, &NANL, 0},
+ {"powl", powl, &ONEL, &NANL, &ONEL, 0},
{"powl", powl, &NANL, &NANL, &NANL, 0},
- {"powl", powl, &ONEL, &INFINITYL, &NANL, 0},
- {"powl", powl, &MONEL, &INFINITYL, &NANL, 0},
- {"powl", powl, &ONEL, &MINFL, &NANL, 0},
- {"powl", powl, &MONEL, &MINFL, &NANL, 0},
+ {"powl", powl, &ONEL, &INFINITYL, &ONEL, 0},
+ {"powl", powl, &MONEL, &INFINITYL, &ONEL, 0},
+ {"powl", powl, &ONEL, &MINFL, &ONEL, 0},
+ {"powl", powl, &MONEL, &MINFL, &ONEL, 0},
{"powl", powl, &MTWOL, &HALFL, &NANL, 0},
{"powl", powl, &ZEROL, &MTHREEL, &INFINITYL, 0},
{"powl", powl, &NEGZEROL, &MTHREEL, &MINFL, 0},
