NetBSD Problem Report #49240
From he@smistad.uninett.no Sun Sep 28 16:24:23 2014
Return-Path: <he@smistad.uninett.no>
Received: from mail.netbsd.org (mail.netbsd.org [149.20.53.66])
(using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits))
(Client CN "mail.netbsd.org", Issuer "Postmaster NetBSD.org" (verified OK))
by mollari.NetBSD.org (Postfix) with ESMTPS id 31E01A6554
for <gnats-bugs@gnats.NetBSD.org>; Sun, 28 Sep 2014 16:24:23 +0000 (UTC)
Message-Id: <20140928162416.DC18C3D0B8@smistad.uninett.no>
Date: Sun, 28 Sep 2014 18:24:16 +0200 (CEST)
From: he@NetBSD.org
Reply-To: he@NetBSD.org
To: gnats-bugs@gnats.NetBSD.org
Subject: Some corner cases for pow() fail to work correctly
X-Send-Pr-Version: 3.95
>Number: 49240
>Category: lib
>Synopsis: Some corner cases for pow() fail to work correctly
>Confidential: no
>Severity: serious
>Priority: medium
>Responsible: lib-bug-people
>State: open
>Class: sw-bug
>Submitter-Id: net
>Arrival-Date: Sun Sep 28 16:25:01 +0000 2014
>Last-Modified: Sat Oct 04 14:25:00 +0000 2014
>Originator: Havard Eidnes
>Release: NetBSD 6.1_STABLE
>Organization:
>Environment:
System: NetBSD smistad.uninett.no 6.1_STABLE NetBSD 6.1_STABLE (MAANEN) #0: Mon Sep 1 14:42:18 CEST 2014 root@smistad.uninett.no:/usr/obj/sys/arch/i386/compile/MAANEN i386
Architecture: i386
Machine: i386
>Description:
While working on pkgsrc-wip/v8's test suite, I found one of
them which failed which excercise the pow() function, by
first doing some simple sanity checking and then trying out
various corner cases.
The v8 test bails out on the first test failure. I've
translated the tests into C, which runs through all the tests
counting up the successes and failures, displaying the result
of each test, with details for the failing tests.
Mostly, the failing tests on NetBSD appear also to contradict
the comments about the "Special cases" near the top of our
lib/libm/src/e_pow.c. The test failures and the special case
confirming the test is OK and our code is not are:
Test failure: NaN == pow(1, Infinity)
a == nan
b == 1.000000
Test failure: NaN == pow(1, -Infinity)
a == nan
b == 1.000000
* 9. +-1 ** +-INF is NAN
Test failure: Infinity == pow(+0.0, -1.1)
a == inf
b == -inf
Test failure: Infinity == pow(+0.0, -2)
a == inf
b == -inf
* 12. +0 ** (-anything except 0, NAN) is +INF
Test failure: Infinity == pow(-0.0, -3.1)
a == inf
b == -inf
Test failure: Infinity == pow(-0.0, -2)
a == inf
b == -inf
Test failure: +Infinity == pow(-0.0, -0.5)
a == inf
b == -inf
Test failure: +Infinity == pow(-0.0, -0.6)
a == inf
b == -inf
* 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
>How-To-Repeat:
Compile and run this program, and watch it fail 8 of our tests
on netbsd-6 and netbsd-7.
(The meat of the code is (c) Google with a 3-clause BSD-like
license, but this is a translation into C by myself.)
------------------------------ snip
#include <string.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define Infinity INFINITY
#define NaN NAN
#define assertEquals(a, b) \
do { \
if ((a) != (b)) { \
char s1[64], s2[64]; \
sprintf(s1, "%f", (double)(a)); \
sprintf(s2, "%f", (double)(b)); \
if ((strcmp(s1, "nan") == 0) || \
(strcmp(s2, "nan") == 0)) { \
if (strcmp(s1, s2) != 0) { \
printf("Test failure: %s == %s\n", #a, #b); \
printf(" a == %s\n", s1); \
printf(" b == %s\n", s2); \
failures++; \
} else { \
printf("Test OK : %s == %s\n", #a, #b); \
ok++; \
} \
} else { \
printf("Test failure: %s == %s\n", #a, #b); \
printf(" a == %f\n", (double)(a)); \
printf(" b == %f\n", (double)(b)); \
failures++; \
} \
} else { \
printf("Test OK : %s == %s\n", #a, #b); \
ok++; \
} \
} while(0)
int
main(int argc, char **argv)
{
int failures = 0;
int ok = 0;
// These two need to work for the tests to work
// Verify the "nan" dance in the assertEquals macro...
assertEquals(NaN, NaN);
assertEquals(Infinity, Infinity);
// Spec tests
assertEquals(NaN, pow(2, NaN));
assertEquals(NaN, pow(+0.0, NaN));
assertEquals(NaN, pow(-0.0, NaN));
assertEquals(NaN, pow(Infinity, NaN));
assertEquals(NaN, pow(-Infinity, NaN));
assertEquals(1, pow(NaN, +0.0));
assertEquals(1, pow(NaN, -0.0));
assertEquals(NaN, pow(NaN, NaN));
assertEquals(NaN, pow(NaN, 2.2));
assertEquals(NaN, pow(NaN, 1));
assertEquals(NaN, pow(NaN, -1));
assertEquals(NaN, pow(NaN, -2.2));
assertEquals(NaN, pow(NaN, Infinity));
assertEquals(NaN, pow(NaN, -Infinity));
assertEquals(Infinity, pow(1.1, Infinity));
assertEquals(Infinity, pow(-1.1, Infinity));
assertEquals(Infinity, pow(2, Infinity));
assertEquals(Infinity, pow(-2, Infinity));
// Because +0 == -0, we need to compare 1/{+,-}0 to {+,-}Infinity
assertEquals(+Infinity, 1/pow(1.1, -Infinity));
assertEquals(+Infinity, 1/pow(-1.1, -Infinity));
assertEquals(+Infinity, 1/pow(2, -Infinity));
assertEquals(+Infinity, 1/pow(-2, -Infinity));
assertEquals(NaN, pow(1, Infinity));
assertEquals(NaN, pow(1, -Infinity));
assertEquals(NaN, pow(-1, Infinity));
assertEquals(NaN, pow(-1, -Infinity));
assertEquals(+0.0, pow(0.1, Infinity));
assertEquals(+0.0, pow(-0.1, Infinity));
assertEquals(+0.0, pow(0.999, Infinity));
assertEquals(+0.0, pow(-0.999, Infinity));
assertEquals(Infinity, pow(0.1, -Infinity));
assertEquals(Infinity, pow(-0.1, -Infinity));
assertEquals(Infinity, pow(0.999, -Infinity));
assertEquals(Infinity, pow(-0.999, -Infinity));
assertEquals(Infinity, pow(Infinity, 0.1));
assertEquals(Infinity, pow(Infinity, 2));
assertEquals(+Infinity, 1/pow(Infinity, -0.1));
assertEquals(+Infinity, 1/pow(Infinity, -2));
assertEquals(-Infinity, pow(-Infinity, 3));
assertEquals(-Infinity, pow(-Infinity, 13));
assertEquals(Infinity, pow(-Infinity, 3.1));
assertEquals(Infinity, pow(-Infinity, 2));
assertEquals(-Infinity, 1/pow(-Infinity, -3));
assertEquals(-Infinity, 1/pow(-Infinity, -13));
assertEquals(+Infinity, 1/pow(-Infinity, -3.1));
assertEquals(+Infinity, 1/pow(-Infinity, -2));
assertEquals(+Infinity, 1/pow(+0.0, 1.1));
assertEquals(+Infinity, 1/pow(+0.0, 2));
assertEquals(Infinity, pow(+0.0, -1.1));
assertEquals(Infinity, pow(+0.0, -2));
assertEquals(-Infinity, 1/pow(-0.0, 3));
assertEquals(-Infinity, 1/pow(-0.0, 13));
assertEquals(+Infinity, 1/pow(-0.0, 3.1));
assertEquals(+Infinity, 1/pow(-0.0, 2));
assertEquals(-Infinity, pow(-0.0, -3));
assertEquals(-Infinity, pow(-0.0, -13));
assertEquals(Infinity, pow(-0.0, -3.1));
assertEquals(Infinity, pow(-0.0, -2));
assertEquals(NaN, pow(-0.00001, 1.1));
assertEquals(NaN, pow(-0.00001, -1.1));
assertEquals(NaN, pow(-1.1, 1.1));
assertEquals(NaN, pow(-1.1, -1.1));
assertEquals(NaN, pow(-2, 1.1));
assertEquals(NaN, pow(-2, -1.1));
assertEquals(NaN, pow(-1000, 1.1));
assertEquals(NaN, pow(-1000, -1.1));
assertEquals(+Infinity, 1/pow(-0.0, 0.5));
assertEquals(+Infinity, 1/pow(-0.0, 0.6));
assertEquals(-Infinity, 1/pow(-0.0, 1));
assertEquals(-Infinity, 1/pow(-0.0, 10000000001.0));
assertEquals(+Infinity, pow(-0.0, -0.5));
assertEquals(+Infinity, pow(-0.0, -0.6));
assertEquals(-Infinity, pow(-0.0, -1));
assertEquals(-Infinity, pow(-0.0, -10000000001.0));
assertEquals(4, pow(16, 0.5));
assertEquals(NaN, pow(-16, 0.5));
assertEquals(0.25, pow(16, -0.5));
assertEquals(NaN, pow(-16, -0.5));
// Test detecting and converting integer value as double.
assertEquals(8, pow(2, sqrt(9)));
// Tests from Mozilla 15.8.2.13.
assertEquals(Infinity, pow(-Infinity, Infinity));
assertEquals(0, pow(-Infinity, -Infinity));
assertEquals(1, pow(0, 0));
assertEquals(0, pow(0, Infinity));
assertEquals(NaN, pow(NaN, 0.5));
assertEquals(NaN, pow(NaN, -0.5));
// Tests from Sputnik S8.5_A13_T1.
assertEquals(
(1*((pow(2,53))-1)*(pow(2,-1074))),
4.4501477170144023e-308);
assertEquals(
(1*(pow(2,52))*(pow(2,-1074))),
2.2250738585072014e-308);
assertEquals(
(-1*(pow(2,52))*(pow(2,-1074))),
-2.2250738585072014e-308);
printf("\n");
printf("Failed tests: %d\n", failures);
printf("Successful tests: %d\n", ok);
exit(failures != 0);
}
------------------------------ snip
>Fix:
Sorry, no fix proposed.
The pcc maintainer tested this with pcc on netbsd-6, and
interestingly these two tests succeed with pcc but fail with
gcc (with the same libm!):
Test OK : NaN == pow(1, Infinity)
Test OK : NaN == pow(1, -Infinity)
Also of interest is that there's been a regression between
netbsd-5 and netbsd-6 on the number of test failures: these
two succeed on netbsd-5:
Test OK : Infinity == pow(-0.0, -2)
Test OK : Infinity == pow(+0.0, -2)
but fail on netbsd-6 or netbsd-7:
Test failure: Infinity == pow(+0.0, -2)
a == inf
b == -inf
Test failure: Infinity == pow(-0.0, -2)
a == inf
b == -inf
This is seen both for macppc and amd64.
>Audit-Trail:
From: Havard Eidnes <he@NetBSD.org>
To: gnats-bugs@NetBSD.org
Cc:
Subject: Re: lib/49240: Some corner cases for pow() fail to work correctly
Date: Tue, 30 Sep 2014 14:20:55 +0200 (CEST)
Hi,
I've dug some more, and it looks like the following four tests:
assertEquals(NaN, pow(1, Infinity));
assertEquals(NaN, pow(1, -Infinity));
assertEquals(NaN, pow(-1, Infinity));
assertEquals(NaN, pow(-1, -Infinity));
are in fact erroneous, at least reportedly if one is to follow
C99, and some (FreeBSD commit messages) claim the revised IEEE
754 of 2008. These are now all defined to give +1.0 as result.
I'll follow that up with Google's V8 project.
It seems our libm is based on fdlibm, but we have no entry in
doc/3RDPARTY documentind that, so "which version" is unknown to
me so far.
Interestingly, the test when run on FreeBSD 10/amd64 failed only
the four above tests, with results which were according to C99,
apparently.
Some pointers:
FreeBSD diff for the C99 update:
http://svnweb.freebsd.org/base/head/lib/msun/src/e_pow.c?r1=176276&r2=268597
A pointer to fdlibm 5.3 and an IEEE test suite from UCB (I've not
looked at that yet -- will do time permitting) can be found at
http://www.validlab.com/software/
There appears to be another bug in fdlibm 5.3's pow function,
talked about here:
https://blogs.oracle.com/darcy/entry/finding_a_bug_in_fdlibm
but unfortunately the URL pointer to the fix is dead. That
document also says that the 4 above tests should give +1.0 in C99.
A minor tweak and copying of our math-private.h makes the FreeBSD
version build outside of the NetBSD source tree, and by using the
result and correcting the test cases for the above 4, I now have
a "0 fails" testrun on NetBSD.
I'll work out a minimal diff for our e_pow.c.
Regards,
- Havard
From: Havard Eidnes <he@NetBSD.org>
To: gnats-bugs@netbsd.org
Cc:
Subject: Re: lib/49240: Some corner cases for pow() fail to work correctly
Date: Tue, 30 Sep 2014 17:23:27 +0200 (CEST)
> I've dug some more, and it looks like the following four tests:
>
> assertEquals(NaN, pow(1, Infinity));
> assertEquals(NaN, pow(1, -Infinity));
> assertEquals(NaN, pow(-1, Infinity));
> assertEquals(NaN, pow(-1, -Infinity));
>
> are in fact erroneous, at least reportedly if one is to follow
> C99, and some (FreeBSD commit messages) claim the revised IEEE
> 754 of 2008. These are now all defined to give +1.0 as result.
> I'll follow that up with Google's V8 project.
The feedback I received there is that ECMAscript is specific
about mandating the old IEEE 754 rules in these 4 test cases.
I'm currently testing a wrapper patch for this issue in v8, so
that C99 systems also can do this "correctly".
Regards,
- Havard
From: Havard Eidnes <he@NetBSD.org>
To: gnats-bugs@NetBSD.org
Cc:
Subject: Re: lib/49240: Some corner cases for pow() fail to work correctly
Date: Wed, 01 Oct 2014 12:26:54 +0200 (CEST)
> I'll work out a minimal diff for our e_pow.c.
This follows below. With that and the test correction for C99
behaviour, the tests now all pass when I test locally.
The changes rolled into this are:
* Change yy1 to y1; no need for yy1.
* Add a case for 1**y = 1, even if y is NaN (or +-Inf)
* Use (x+0.0)+(y+0.0) instead of x+y when mixing NaNs. There's
a long explanation of why this is done in FreeBSD's commit log
for e_pow.c at
http://svnweb.freebsd.org/base/head/lib/msun/src/e_pow.c?view=log
* Add a case for -1**+-Inf, since C99 and IEEE 754 revised 2008(?)
reportedly state that this should be 1
* Add a fix to get correct sign for some over- and underflow
cases. This seems to come from fdlib 5.3, FreeBSD revision
129956. (Ref. above SVN log.)
* Add a fix to get well-defined behaviour when right-shifting a
possibly negative value by adding a cast.
* Fix comment detailing the expected result for the special cases
--- /usr/src/lib/libm/src/e_pow.c 2013-03-27 19:12:07.000000000 +0100
+++ e_pow.c 2014-10-01 11:41:10.000000000 +0200
@@ -29,13 +29,13 @@
* Special cases:
* 1. (anything) ** 0 is 1
* 2. (anything) ** 1 is itself
- * 3. (anything) ** NAN is NAN
+ * 3. (anything) ** NAN is NAN except 1 ** NAN = 1
* 4. NAN ** (anything except 0) is NAN
* 5. +-(|x| > 1) ** +INF is +INF
* 6. +-(|x| > 1) ** -INF is +0
* 7. +-(|x| < 1) ** +INF is +0
* 8. +-(|x| < 1) ** -INF is +INF
- * 9. +-1 ** +-INF is NAN
+ * 9. +-1 ** +-INF is 1
* 10. +0 ** (+anything except 0, NAN) is +0
* 11. -0 ** (+anything except 0, NAN, odd integer) is +0
* 12. +0 ** (-anything except 0, NAN) is +INF
@@ -98,10 +98,10 @@
ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
double
-__ieee754_pow(double x, double y)
+pow(double x, double y)
{
double z,ax,z_h,z_l,p_h,p_l;
- double yy1,t1,t2,r,s,t,u,v,w;
+ double y1,t1,t2,r,s,t,u,v,w;
int32_t i,j,k,yisint,n;
int32_t hx,hy,ix,iy;
u_int32_t lx,ly;
@@ -113,10 +113,13 @@
/* y==zero: x**0 = 1 */
if((iy|ly)==0) return one;
- /* +-NaN return x+y */
+ /* x==1: 1**y = 1, even if y is NaN */
+ if (hx==0x3ff00000 && lx == 0) return one;
+
+ /* y!=zero: result is NaN if either arg is NaN */
if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
- return x+y;
+ return (x+0.0)+(y+0.0);
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer
@@ -142,7 +145,7 @@
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 1 */
else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
return (hy>=0)? y: zero;
else /* (|x|<1)**-,+inf = inf,0 */
@@ -174,7 +177,11 @@
}
}
+ /* CYGNUS LOCAL + fdlibm-5.3 fix: This used to be
n = (hx>>31)+1;
+ but ANSI C says a right shift of a signed negative quantity is
+ implementation defined. */
+ n = ((u_int32_t)hx>>31)-1;
/* (x<0)**(non-int) is NaN */
if((n|yisint)==0) return (x-x)/(x-x);
@@ -250,11 +257,11 @@
t2 = z_l-(((t1-t)-dp_h[k])-z_h);
}
- /* split up y into yy1+y2 and compute (yy1+y2)*(t1+t2) */
- yy1 = y;
- SET_LOW_WORD(yy1,0);
- p_l = (y-yy1)*t1+y*t2;
- p_h = yy1*t1;
+ /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
+ y1 = y;
+ SET_LOW_WORD(y1,0);
+ p_l = (y-y1)*t1+y*t2;
+ p_h = y1*t1;
z = p_l+p_h;
EXTRACT_WORDS(j,i,z);
if (j>=0x40900000) { /* z >= 1024 */
From: christos@zoulas.com (Christos Zoulas)
To: gnats-bugs@NetBSD.org, lib-bug-people@netbsd.org,
gnats-admin@netbsd.org, netbsd-bugs@netbsd.org, he@NetBSD.org
Cc:
Subject: Re: lib/49240: Some corner cases for pow() fail to work correctly
Date: Wed, 1 Oct 2014 08:39:16 -0400
On Oct 1, 10:30am, he@NetBSD.org (Havard Eidnes) wrote:
-- Subject: Re: lib/49240: Some corner cases for pow() fail to work correctly
| The changes rolled into this are:
|
| * Change yy1 to y1; no need for yy1.
Please change that back. Compilers warn:
/usr/include/math.h:double y1(double);
christos
From: Havard Eidnes <he@NetBSD.org>
To: gnats-bugs@NetBSD.org
Cc:
Subject: Re: lib/49240: Some corner cases for pow() fail to work correctly
Date: Sat, 04 Oct 2014 16:24:18 +0200 (CEST)
Hmm.
It appears that we build a "multi-standard fdlibm as libm", and that
we by default implement _POSIX_MODE (ref. src/lib/lim/Makefile).
So ... the test program, if it should test the IEEE compliance for the
corner cases, needs to tell libm to do the IEEE operations instead of
the POSIX operations via
_LIB_VERSION =3D _IEEE_;
(is that the officially sanctioned method, BTW?)
and the only failures flagged by the test program I posted earlier
which our original libm then fail with is
Test failure: 1.0 =3D=3D pow(-1, Infinity)
a =3D=3D 1.000000
b =3D=3D nan
Test failure: 1.0 =3D=3D pow(-1, -Infinity)
a =3D=3D 1.000000
b =3D=3D nan
which correspond to this change:
* Add a case for -1**+-Inf, since C99 and IEEE 754 revised 2008(?)
reportedly state that this should be 1
Further, it looks like the reason for the "yy1" variable name is that
a use of "y1" appears to mask a global "y1", at least on netbsd-6.
Test cases to justify the other changes should be found.
Regards,
- H=E5vard
(Contact us)
$NetBSD: query-full-pr,v 1.39 2013/11/01 18:47:49 spz Exp $
$NetBSD: gnats_config.sh,v 1.8 2006/05/07 09:23:38 tsutsui Exp $
Copyright © 1994-2014
The NetBSD Foundation, Inc. ALL RIGHTS RESERVED.
Home