[BZ #258]

Update.
2004-07-19  Jakub Jelinek  <jakub@redhat.com>

	[BZ #258]
	* math/libm-test.inc (max_value, min_value): New variables.
	(initialize): Initialize them.
	(pow_test): Add a couple of new tests.
	* sysdeps/i386/fpu/e_powf.S (__ieee754_powf): Don't generate invalid
	exception if |y| >= 1U<<31.
	* sysdeps/i386/fpu/e_pow.S (__ieee754_pow): Don't generate invalid
	exception if |y| >= 1L<<63.
	* sysdeps/i386/fpu/e_powl.S (__ieee754_powl): Likewise.
	If y*log2(x) overflows to +-inf, return still +inf/+0 instead of NaN.
	* sysdeps/x86_64/fpu/e_powl.S (__ieee754_powl): Likewise.

diff --git a/ChangeLog b/ChangeLog
index 6083c205e7ff9cb4c110ec5bbce68c0cd01c2df8..5f4fc034a5752671cee31c42f48f89bbf7d5d71d 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,17 @@
+2004-07-19  Jakub Jelinek  <jakub@redhat.com>
+
+       [BZ #258]
+       * math/libm-test.inc (max_value, min_value): New variables.
+       (initialize): Initialize them.
+       (pow_test): Add a couple of new tests.
+       * sysdeps/i386/fpu/e_powf.S (__ieee754_powf): Don't generate invalid
+       exception if |y| >= 1U<<31.
+       * sysdeps/i386/fpu/e_pow.S (__ieee754_pow): Don't generate invalid
+       exception if |y| >= 1L<<63.
+       * sysdeps/i386/fpu/e_powl.S (__ieee754_powl): Likewise.
+       If y*log2(x) overflows to +-inf, return still +inf/+0 instead of NaN.
+       * sysdeps/x86_64/fpu/e_powl.S (__ieee754_powl): Likewise.
+
 2004-07-18  Ulrich Drepper  <drepper@redhat.com>
 
        * nscd/pwdcache.c (cache_addpw): Optimize case of unsuccessful

diff --git a/math/libm-test.inc b/math/libm-test.inc
index bad3948f1dea125529389c1f0697e80af67def40..729fce936993073b69bf0926f01630a3b362115c 100644 (file)
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -169,7 +169,7 @@ static int output_points;   /* Should the single function results printed?  */
 static int ignore_max_ulp;     /* Should we ignore max_ulp?  */
 
 static FLOAT minus_zero, plus_zero;
-static FLOAT plus_infty, minus_infty, nan_value;
+static FLOAT plus_infty, minus_infty, nan_value, max_value, min_value;
 
 static FLOAT max_error, real_max_error, imag_max_error;
 
@@ -3593,6 +3593,28 @@ pow_test (void)
   TEST_ff_f (pow, -1, plus_infty, 1);
   TEST_ff_f (pow, 1, minus_infty, 1);
   TEST_ff_f (pow, -1, minus_infty, 1);
+  TEST_ff_f (pow, 1, 1, 1);
+  TEST_ff_f (pow, 1, -1, 1);
+  TEST_ff_f (pow, 1, 1.25, 1);
+  TEST_ff_f (pow, 1, -1.25, 1);
+  TEST_ff_f (pow, 1, 0x1p62L, 1);
+  TEST_ff_f (pow, 1, 0x1p63L, 1);
+  TEST_ff_f (pow, 1, 0x1p64L, 1);
+  TEST_ff_f (pow, 1, 0x1p72L, 1);
+
+  /* pow (x, +-0) == 1.  */
+  TEST_ff_f (pow, plus_infty, 0, 1);
+  TEST_ff_f (pow, plus_infty, minus_zero, 1);
+  TEST_ff_f (pow, minus_infty, 0, 1);
+  TEST_ff_f (pow, minus_infty, minus_zero, 1);
+  TEST_ff_f (pow, 32.75L, 0, 1);
+  TEST_ff_f (pow, 32.75L, minus_zero, 1);
+  TEST_ff_f (pow, -32.75L, 0, 1);
+  TEST_ff_f (pow, -32.75L, minus_zero, 1);
+  TEST_ff_f (pow, 0x1p72L, 0, 1);
+  TEST_ff_f (pow, 0x1p72L, minus_zero, 1);
+  TEST_ff_f (pow, 0x1p-72L, 0, 1);
+  TEST_ff_f (pow, 0x1p-72L, minus_zero, 1);
 
   TEST_ff_f (pow, -0.1L, 1.1L, nan_value, INVALID_EXCEPTION);
   TEST_ff_f (pow, -0.1L, -1.1L, nan_value, INVALID_EXCEPTION);
@@ -3609,6 +3631,10 @@ pow_test (void)
   TEST_ff_f (pow, minus_zero, -2, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
   TEST_ff_f (pow, minus_zero, -11.1L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
 
+  TEST_ff_f (pow, 0x1p72L, 0x1p72L, plus_infty);
+  TEST_ff_f (pow, 10, -0x1p72L, 0);
+  TEST_ff_f (pow, max_value, max_value, plus_infty);
+  TEST_ff_f (pow, 10, -max_value, 0);
 
   TEST_ff_f (pow, 0, 1, 0);
   TEST_ff_f (pow, 0, 11, 0);
@@ -3623,6 +3649,8 @@ pow_test (void)
 
   TEST_ff_f (pow, minus_zero, 2, 0);
   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);
 
 #ifndef TEST_INLINE
   /* pow (x, +inf) == +inf for |x| > 1.  */
@@ -3667,6 +3695,11 @@ pow_test (void)
   /* pow (-0, y) == +0 for y > 0 and not an odd integer.  */
   TEST_ff_f (pow, minus_zero, 4, 0.0);
 
+  TEST_ff_f (pow, 16, 0.25L, 2);
+  TEST_ff_f (pow, 0x1p64L, 0.125L, 256);
+  TEST_ff_f (pow, 2, 4, 16);
+  TEST_ff_f (pow, 256, 8, 0x1p64L);
+
   TEST_ff_f (pow, 0.75L, 1.25L, 0.697953644326574699205914060237425566L);
 
 #if defined TEST_DOUBLE || defined TEST_LDOUBLE
@@ -4312,12 +4345,18 @@ initialize (void)
                       HUGE_VALL, HUGE_VAL, HUGE_VALF);
   minus_infty = CHOOSE (-HUGE_VALL, -HUGE_VAL, -HUGE_VALF,
                        -HUGE_VALL, -HUGE_VAL, -HUGE_VALF);
+  max_value = CHOOSE (LDBL_MAX, DBL_MAX, FLT_MAX,
+                     LDBL_MAX, DBL_MAX, FLT_MAX);
+  min_value = CHOOSE (LDBL_MIN, DBL_MIN, FLT_MIN,
+                     LDBL_MIN, DBL_MIN, FLT_MIN);
 
   (void) &plus_zero;
   (void) &nan_value;
   (void) &minus_zero;
   (void) &plus_infty;
   (void) &minus_infty;
+  (void) &max_value;
+  (void) &min_value;
 
   /* Clear all exceptions.  From now on we must not get random exceptions.  */
   feclearexcept (FE_ALL_EXCEPT);

diff --git a/sysdeps/i386/fpu/e_pow.S b/sysdeps/i386/fpu/e_pow.S
index 997dd30e9dba4799350d6a14ca124eef9a35c264..d100c0ac571dedca3f425ee94c716e81abb077b2 100644 (file)
--- a/sysdeps/i386/fpu/e_pow.S
+++ b/sysdeps/i386/fpu/e_pow.S
@@ -1,5 +1,6 @@
 /* ix87 specific implementation of pow function.
-   Copyright (C) 1996, 1997, 1998, 1999, 2001 Free Software Foundation, Inc.
+   Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004
+   Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
 
@@ -48,6 +49,9 @@ one:  .double 1.0
        ASM_TYPE_DIRECTIVE(limit,@object)
 limit: .double 0.29
        ASM_SIZE_DIRECTIVE(limit)
+       ASM_TYPE_DIRECTIVE(p63,@object)
+p63:   .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
+       ASM_SIZE_DIRECTIVE(p63)
 
 #ifdef PIC
 #define MO(op) op##@GOTOFF(%ecx)
@@ -96,6 +100,14 @@ ENTRY(__ieee754_pow)
 
        fxch                    // y : x
 
+       /* fistpll raises invalid exception for |y| >= 1L<<63.  */
+       fld     %st             // y : y : x
+       fabs                    // |y| : y : x
+       fcompl  MO(p63)         // y : x
+       fnstsw
+       sahf
+       jnc     2f
+
        /* First see whether `y' is a natural number.  In this case we
           can use a more precise algorithm.  */
        fld     %st             // y : y : x

diff --git a/sysdeps/i386/fpu/e_powl.S b/sysdeps/i386/fpu/e_powl.S
index 5d825711288f932af5ec877d20544cc8cfe96655..080764b84a6296f7489df1548e063b7cfbb3285d 100644 (file)
--- a/sysdeps/i386/fpu/e_powl.S
+++ b/sysdeps/i386/fpu/e_powl.S
@@ -1,5 +1,6 @@
 /* ix87 specific implementation of pow function.
-   Copyright (C) 1996, 1997, 1998, 1999, 2001 Free Software Foundation, Inc.
+   Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004
+   Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
 
@@ -48,6 +49,9 @@ one:  .double 1.0
        ASM_TYPE_DIRECTIVE(limit,@object)
 limit: .double 0.29
        ASM_SIZE_DIRECTIVE(limit)
+       ASM_TYPE_DIRECTIVE(p63,@object)
+p63:   .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
+       ASM_SIZE_DIRECTIVE(p63)
 
 #ifdef PIC
 #define MO(op) op##@GOTOFF(%ecx)
@@ -96,6 +100,14 @@ ENTRY(__ieee754_powl)
 
        fxch                    // y : x
 
+       /* fistpll raises invalid exception for |y| >= 1L<<63.  */
+       fld     %st             // y : y : x
+       fabs                    // |y| : y : x
+       fcompl  MO(p63)         // y : x
+       fnstsw
+       sahf
+       jnc     2f
+
        /* First see whether `y' is a natural number.  In this case we
           can use a more precise algorithm.  */
        fld     %st             // y : y : x
@@ -161,6 +173,11 @@ ENTRY(__ieee754_powl)
 
 7:     fyl2x                   // log2(x) : y
 8:     fmul    %st(1)          // y*log2(x) : y
+       fxam
+       fnstsw
+       andb    $0x45, %ah
+       cmpb    $0x05, %ah      // is y*log2(x) == ±inf ?
+       je      28f
        fst     %st(1)          // y*log2(x) : y*log2(x)
        frndint                 // int(y*log2(x)) : y*log2(x)
        fsubr   %st, %st(1)     // int(y*log2(x)) : fract(y*log2(x))
@@ -172,6 +189,12 @@ ENTRY(__ieee754_powl)
        fstp    %st(1)          // 2^fract(y*log2(x))*2^int(y*log2(x))
        ret
 
+28:    fstp    %st(1)          // y*log2(x)
+       fldl    MO(one)         // 1 : y*log2(x)
+       fscale                  // 2^(y*log2(x)) : y*log2(x)
+       addl    $8, %esp
+       fstp    %st(1)          // 2^(y*log2(x))
+       ret
 
        // pow(x,±0) = 1
        .align ALIGNARG(4)

diff --git a/sysdeps/x86_64/fpu/e_powl.S b/sysdeps/x86_64/fpu/e_powl.S
index 8c690e16ccfa143245fe77a81d8910342ff0ff98..85f4deb3c7e0772db16c8dd5f5a0a061c6906aaf 100644 (file)
--- a/sysdeps/x86_64/fpu/e_powl.S
+++ b/sysdeps/x86_64/fpu/e_powl.S
@@ -1,5 +1,5 @@
 /* ix87 specific implementation of pow function.
-   Copyright (C) 1996, 1997, 1998, 1999, 2001 Free Software Foundation, Inc.
+   Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
 
@@ -48,6 +48,10 @@ one: .double 1.0
        ASM_TYPE_DIRECTIVE(limit,@object)
 limit: .double 0.29
        ASM_SIZE_DIRECTIVE(limit)
+       ASM_TYPE_DIRECTIVE(p63,@object)
+p63:
+       .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
+       ASM_SIZE_DIRECTIVE(p63)
 
 #ifdef PIC
 #define MO(op) op##(%rip)
@@ -87,6 +91,14 @@ ENTRY(__ieee754_powl)
 
        fxch                    // y : x
 
+       /* fistpll raises invalid exception for |y| >= 1L<<63.  */
+       fldl    MO(p63)         // 1L<<63 : y : x
+       fld     %st(1)          // y : 1L<<63 : y : x
+       fabs                    // |y| : 1L<<63 : y : x
+       fcomip  %st(1), %st     // 1L<<63 : y : x
+       fstp    %st(0)          // y : x
+       jnc     2f
+
        /* First see whether `y' is a natural number.  In this case we
           can use a more precise algorithm.  */
        fld     %st             // y : y : x
@@ -148,6 +160,11 @@ ENTRY(__ieee754_powl)
 
 7:     fyl2x                   // log2(x) : y
 8:     fmul    %st(1)          // y*log2(x) : y
+       fxam
+       fnstsw
+       andb    $0x45, %ah
+       cmpb    $0x05, %ah      // is y*log2(x) == ±inf ?
+       je      28f
        fst     %st(1)          // y*log2(x) : y*log2(x)
        frndint                 // int(y*log2(x)) : y*log2(x)
        fsubr   %st, %st(1)     // int(y*log2(x)) : fract(y*log2(x))
@@ -158,6 +175,11 @@ ENTRY(__ieee754_powl)
        fstp    %st(1)          // 2^fract(y*log2(x))*2^int(y*log2(x))
        ret
 
+28:    fstp    %st(1)          // y*log2(x)
+       fldl    MO(one)         // 1 : y*log2(x)
+       fscale                  // 2^(y*log2(x)) : y*log2(x)
+       fstp    %st(1)          // 2^(y*log2(x))
+       ret
 
        // pow(x,±0) = 1
        .align ALIGNARG(4)