diff options
author | Thomas Fitzsimmons <fitzsim@redhat.com> | 2002-06-27 20:41:49 +0000 |
---|---|---|
committer | Thomas Fitzsimmons <fitzsim@redhat.com> | 2002-06-27 20:41:49 +0000 |
commit | 54be629f41ad4c5e46a736cf988535a794c664bd (patch) | |
tree | f285adadcb54176986d0002d17d1677e2439a8b5 /newlib/libm | |
parent | c36e6dd754453b8f57767b19c58d2f832bac8bb0 (diff) | |
download | newlib-54be629f41ad4c5e46a736cf988535a794c664bd.zip newlib-54be629f41ad4c5e46a736cf988535a794c664bd.tar.gz newlib-54be629f41ad4c5e46a736cf988535a794c664bd.tar.bz2 |
* libm/mathfp/s_pow.c (pow): Fix checks on variable k. Add
exponent_is_even_int variable. Handle case where x is
negative, and y is an odd integer.
* libm/mathfp/sf_pow.c (powf): Likewise.
Diffstat (limited to 'newlib/libm')
-rw-r--r-- | newlib/libm/mathfp/s_pow.c | 71 | ||||
-rw-r--r-- | newlib/libm/mathfp/sf_pow.c | 65 |
2 files changed, 85 insertions, 51 deletions
diff --git a/newlib/libm/mathfp/s_pow.c b/newlib/libm/mathfp/s_pow.c index 7c0a38a..b4b6c89 100644 --- a/newlib/libm/mathfp/s_pow.c +++ b/newlib/libm/mathfp/s_pow.c @@ -52,57 +52,67 @@ PORTABILITY double pow (double x, double y) { - double d, t, r = 1.0; - int n, k, sign = 0; + double d, k, t, r = 1.0; + int n, sign, exponent_is_even_int = 0; __uint32_t px; GET_HIGH_WORD (px, x); k = modf (y, &d); - if (k == 0.0) + + if (k == 0.0) { + /* Exponent y is an integer. */ if (modf (ldexp (y, -1), &t)) - sign = 0; + { + /* y is odd. */ + exponent_is_even_int = 0; + } else - sign = 1; + { + /* y is even. */ + exponent_is_even_int = 1; + } } - if (x == 0.0 && y <= 0.0) - errno = EDOM; - + if (x == 0.0 && y <= 0.0) + { + errno = EDOM; + } else if ((t = y * log (fabs (x))) >= BIGX) { errno = ERANGE; if (px & 0x80000000) { - if (!k) + /* x is negative. */ + if (k) { + /* y is not an integer. */ errno = EDOM; x = 0.0; } - else if (sign) - x = -z_infinity.d; + else if (exponent_is_even_int) + x = z_infinity.d; else - x = z_infinity.d; + x = -z_infinity.d; } - - else - x = z_infinity.d; - } - + else + { + x = z_infinity.d; + } + } else if (t < SMALLX) { errno = ERANGE; x = 0.0; } - else { - if ( k && fabs(d) <= 32767 ) + if ( !k && fabs(d) <= 32767 ) { n = (int) d; - if (sign = (n < 0)) + if ((sign = (n < 0))) n = -n; while ( n > 0 ) @@ -118,13 +128,14 @@ double pow (double x, double y) return r; } - else { if ( px & 0x80000000 ) { - if ( !k ) + /* x is negative. */ + if ( k ) { + /* y is not an integer. */ errno = EDOM; return 0.0; } @@ -132,13 +143,19 @@ double pow (double x, double y) x = exp (t); - if ( sign ) - { - px ^= 0x80000000; - SET_HIGH_WORD (x, px); + if (!exponent_is_even_int) + { + if (px & 0x80000000) + { + /* y is an odd integer, and x is negative, + so the result is negative. */ + GET_HIGH_WORD (px, x); + px |= 0x80000000; + SET_HIGH_WORD (x, px); + } } } - } + } return x; } diff --git a/newlib/libm/mathfp/sf_pow.c b/newlib/libm/mathfp/sf_pow.c index 2b3bed3..7f40186 100644 --- a/newlib/libm/mathfp/sf_pow.c +++ b/newlib/libm/mathfp/sf_pow.c @@ -7,56 +7,66 @@ float powf (float x, float y) { float d, t, r = 1.0; - int n, k, sign = 0; + int n, k, sign = 0, exponent_is_even_int = 0; __int32_t px; GET_FLOAT_WORD (px, x); k = modff (y, &d); + if (k == 0.0) { + /* Exponent y is an integer. */ if (modff (ldexpf (y, -1), &t)) - sign = 0; + { + /* y is odd. */ + exponent_is_even_int = 0; + } else - sign = 1; + { + /* y is even. */ + exponent_is_even_int = 1; + } } - if (x == 0.0 && y <= 0.0) - errno = EDOM; - + if (x == 0.0 && y <= 0.0) + { + errno = EDOM; + } else if ((t = y * log (fabsf (x))) >= BIGX) { errno = ERANGE; if (px & 0x80000000) { - if (!k) + /* x is negative. */ + if (k) { + /* y is not an integer. */ errno = EDOM; x = 0.0; } - else if (sign) - x = -z_infinity_f.f; + else if (exponent_is_even_int) + x = z_infinity_f.f; else - x = z_infinity_f.f; + x = -z_infinity_f.f; } - - else - x = z_infinity_f.f; - } - + else + { + x = z_infinity_f.f; + } + } else if (t < SMALLX) { errno = ERANGE; x = 0.0; } - else { - if ( k && fabsf (d) <= 32767 ) + if ( !k && fabsf (d) <= 32767 ) { n = (int) d; - if (sign = (n < 0)) + if ((sign = (n < 0))) n = -n; while ( n > 0 ) @@ -72,13 +82,14 @@ float powf (float x, float y) return r; } - else { if ( px & 0x80000000 ) { - if ( !k ) + /* x is negative. */ + if (k) { + /* y is not an integer. */ errno = EDOM; return 0.0; } @@ -86,13 +97,19 @@ float powf (float x, float y) x = exp (t); - if ( sign ) + if (!exponent_is_even_int) { - px ^= 0x80000000; - SET_FLOAT_WORD (x, px); + if (px & 0x80000000) + { + /* y is an odd integer, and x is negative, + so the result is negative. */ + GET_FLOAT_WORD (px, x); + px |= 0x80000000; + SET_FLOAT_WORD (x, px); + } } } - } + } return x; } |