Fix pow (0.0, -Inf) (bug 14241)

[Andreas Jaeger]

On 06/14/2012 10:57 PM, Joseph S. Myers wrote:
> Bug 14241 is a regression in pow where pow (0.0, -Inf) returns the
> wrong result on architectures using the dbl-64/e_pow.c
> implementation.  The problem is in the use of ABS(x) (inside the case
> where x is a zero); the ABS macro in mydefs.h is
>
> #define ABS(x)   (((x)>0)?(x):-(x))
>
> which gives -0.0 when x is 0.0.  (This also affects the similar code
> in the case for finite negative powers with absolute value up to 1e20
> and not odd integers, but this doesn't show as testsuite failures
> because in those cases the correct return value ends up coming from
> __kernel_standard.)
>
> I propose this patch fixing this by using 0.0 instead of ABS(x) in the
> relevant places (not changing the definition of the ABS macro, given
> that I don't know the reasons behind the present choice and given all
> the other uses of that macro).  Tested x86 and x86_64.  Note that the
> use of DIVIDE_BY_ZERO_EXCEPTION_OK is in accordance with C11 where
> this exception is permitted but not required in this case (a
> correction to C99 following N1515; ideally the exception would not be
> present for these cases of infinite exponent, but now is not the time
> to ensure and test consistency in this regard across the various pow
> implementations).
>
> 2012-06-14  Joseph Myers  <joseph@codesourcery.com>
>
> 	[BZ #14241]
> 	* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Use 0.0 instead
> 	of ABS(x) in calculating zero to negative powers other than odd
> 	integers.
> 	* math/libm-test.inc (pow_test): Add more tests.
>
> diff --git a/math/libm-test.inc b/math/libm-test.inc
> index d32dde3..8e4d02e 100644
> --- a/math/libm-test.inc
> +++ b/math/libm-test.inc
> @@ -6290,6 +6290,8 @@ pow_test (void)
>     TEST_ff_f (pow, minus_zero, 11.1L, 0);
>     TEST_ff_f (pow, 0, plus_infty, 0);
>     TEST_ff_f (pow, minus_zero, plus_infty, 0);
> +  TEST_ff_f (pow, 0, minus_infty, plus_infty, DIVIDE_BY_ZERO_EXCEPTION_OK);
> +  TEST_ff_f (pow, minus_zero, minus_infty, plus_infty, DIVIDE_BY_ZERO_EXCEPTION_OK);
>
>   #ifndef TEST_INLINE
>     /* pow (x, +inf) == +inf for |x| > 1.  */
> diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c
> index 6c41af9..3fd5e65 100644
> --- a/sysdeps/ieee754/dbl-64/e_pow.c
> +++ b/sysdeps/ieee754/dbl-64/e_pow.c
> @@ -112,12 +112,12 @@ __ieee754_pow(double x, double y) {
>       if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0)
>   	|| (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000)
>         return y;
> -    if (ABS(y) > 1.0e20) return (y>0)?0:1.0/ABS(x);
> +    if (ABS(y) > 1.0e20) return (y>0)?0:1.0/0.0;
>       k = checkint(y);
>       if (k == -1)
>         return y < 0 ? 1.0/x : x;
>       else
> -      return y < 0 ? 1.0/ABS(x) : 0.0;                               /* return 0 */
> +      return y < 0 ? 1.0/0.0 : 0.0;                               /* return 0 */
>     }
>
>     qx = u.i[HIGH_HALF]&0x7fffffff;  /*   no sign   */
>

Your previous change was:

commit 2460d3aa21f04cdf28497683bd3e29183189f779
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Wed Mar 21 12:14:57 2012 +0000

     Fix pow of zero and infinity to large powers.

diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c
index f936a72..26ffaad 100644
--- a/sysdeps/ieee754/dbl-64/e_pow.c
+++ b/sysdeps/ieee754/dbl-64/e_pow.c
@@ -111,7 +111,7 @@ __ieee754_pow(double x, double y) {
      if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] 
!= 0)
         || (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000)
        return y;
-    if (ABS(y) > 1.0e20) return (y>0)?0:INF.x;
+    if (ABS(y) > 1.0e20) return (y>0)?0:1.0/ABS(x);
      k = checkint(y);
      if (k == -1)
        return y < 0 ? 1.0/x : x;

Can't we just revert this - just using INF.x in both places? Why do we 
need the division here?

Andreas
-- 
  Andreas Jaeger aj@{suse.com,opensuse.org} Twitter/Identica: jaegerandi
   SUSE LINUX Products GmbH, Maxfeldstr. 5, 90409 Nürnberg, Germany
    GF: Jeff Hawn,Jennifer Guild,Felix Imendörffer,HRB16746 (AG Nürnberg)
     GPG fingerprint = 93A3 365E CE47 B889 DF7F  FED1 389A 563C C272 A126