aboutsummaryrefslogtreecommitdiff
path: root/newlib/libm
diff options
context:
space:
mode:
authorThomas Fitzsimmons <fitzsim@redhat.com>2002-06-27 20:41:49 +0000
committerThomas Fitzsimmons <fitzsim@redhat.com>2002-06-27 20:41:49 +0000
commit54be629f41ad4c5e46a736cf988535a794c664bd (patch)
treef285adadcb54176986d0002d17d1677e2439a8b5 /newlib/libm
parentc36e6dd754453b8f57767b19c58d2f832bac8bb0 (diff)
downloadnewlib-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.c71
-rw-r--r--newlib/libm/mathfp/sf_pow.c65
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;
}