fluencelabs/musl
1
0
Fork 0
mirror of https://github.com/fluencelabs/musl synced 2025-06-30 23:21:56 +00:00
Code Issues Projects Releases Wiki Activity

clean up pow.c and powf.c

Browse Source 
fix comments about special cases
This commit is contained in:
nsz
2012-03-20 20:04:53 +01:00
parent 615bbd365f
commit f1347a3a45
2 changed files with 43 additions and 57 deletions
Download Patch File Download Diff File Expand all files Collapse all files

65
src/math/pow.c
View File

@ -21,24 +21,28 @@
* *
* Special cases: * Special cases:
* 1. (anything) ** 0 is 1 * 1. (anything) ** 0 is 1
* 2. (anything) ** 1 is itself * 2. 1 ** (anything) is 1
* 3. (anything except 1) ** NAN is NAN, 1 ** NAN is 1 * 3. (anything except 1) ** NAN is NAN
* 4. NAN ** (anything except 0) is NAN * 4. NAN ** (anything except 0) is NAN
* 5. +-(|x| > 1) ** +INF is +INF * 5. +-(|x| > 1) ** +INF is +INF
* 6. +-(|x| > 1) ** -INF is +0 * 6. +-(|x| > 1) ** -INF is +0
* 7. +-(|x| < 1) ** +INF is +0 * 7. +-(|x| < 1) ** +INF is +0
* 8. +-(|x| < 1) ** -INF is +INF * 8. +-(|x| < 1) ** -INF is +INF
* 9. +-1 ** +-INF is 1 * 9. -1 ** +-INF is 1
* 10. +0 ** (+anything except 0, NAN) is +0 * 10. +0 ** (+anything except 0, NAN) is +0
* 11. -0 ** (+anything except 0, NAN, odd integer) is +0 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0
* 12. +0 ** (-anything except 0, NAN) is +INF * 12. +0 ** (-anything except 0, NAN) is +INF, raise divbyzero
* 13. -0 ** (-anything except 0, NAN, odd integer) is +INF * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF, raise divbyzero
* 14. -0 ** (odd integer) = -( +0 ** (odd integer) ) * 14. -0 ** (+odd integer) is -0
* 15. +INF ** (+anything except 0,NAN) is +INF * 15. -0 ** (-odd integer) is -INF, raise divbyzero
* 16. +INF ** (-anything except 0,NAN) is +0 * 16. +INF ** (+anything except 0,NAN) is +INF
* 17. -INF ** (anything) = -0 ** (-anything) * 17. +INF ** (-anything except 0,NAN) is +0
* 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) * 18. -INF ** (+odd integer) is -INF
* 19. (-anything except 0 and inf) ** (non-integer) is NAN * 19. -INF ** (anything) = -0 ** (-anything), (anything except odd integer)
* 20. (anything) ** 1 is (anything)
* 21. (anything) ** -1 is 1/(anything)
* 22. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
* 23. (-anything except 0 and inf) ** (non-integer) is NAN
* *
* Accuracy: * Accuracy:
* pow(x,y) returns x**y nearly rounded. In particular * pow(x,y) returns x**y nearly rounded. In particular
@ -98,18 +102,16 @@ double pow(double x, double y)
ix = hx & 0x7fffffff; ix = hx & 0x7fffffff;
iy = hy & 0x7fffffff; iy = hy & 0x7fffffff;
/* y == 0.0: x**0 = 1 */ /* x**0 = 1, even if x is NaN */
if ((iy|ly) == 0) if ((iy|ly) == 0)
return 1.0; return 1.0;
/* 1**y = 1, even if y is NaN */
/* x == 1: 1**y = 1, even if y is NaN */
if (hx == 0x3ff00000 && lx == 0) if (hx == 0x3ff00000 && lx == 0)
return 1.0; return 1.0;
/* NaN if either arg is NaN */
/* y != 0.0: result is NaN if either arg is NaN */
if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0) || if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0) ||
iy > 0x7ff00000 || (iy == 0x7ff00000 && ly != 0)) iy > 0x7ff00000 || (iy == 0x7ff00000 && ly != 0))
return (x+0.0) + (y+0.0); return x + y;
/* determine if y is an odd int when x < 0 /* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer * yisint = 0 ... y is not an integer
@ -142,13 +144,10 @@ double pow(double x, double y)
else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */ else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
return hy >= 0 ? y : 0.0; return hy >= 0 ? y : 0.0;
else /* (|x|<1)**+-inf = 0,inf */ else /* (|x|<1)**+-inf = 0,inf */
return hy < 0 ? -y : 0.0; return hy >= 0 ? 0.0 : -y;
}
if (iy == 0x3ff00000) { /* y is +-1 */
if (hy < 0)
return 1.0/x;
return x;
} }
if (iy == 0x3ff00000) /* y is +-1 */
return hy >= 0 ? x : 1.0/x;
if (hy == 0x40000000) /* y is 2 */ if (hy == 0x40000000) /* y is 2 */
return x*x; return x*x;
if (hy == 0x3fe00000) { /* y is 0.5 */ if (hy == 0x3fe00000) { /* y is 0.5 */
@ -174,19 +173,13 @@ double pow(double x, double y)
} }
} }
/* CYGNUS LOCAL + fdlibm-5.3 fix: This used to be s = 1.0; /* sign of result */
n = (hx>>31)+1; if (hx < 0) {
but ANSI C says a right shift of a signed negative quantity is if (yisint == 0) /* (x<0)**(non-int) is NaN */
implementation defined. */ return (x-x)/(x-x);
n = ((uint32_t)hx>>31) - 1; if (yisint == 1) /* (x<0)**(odd int) */
s = -1.0;
/* (x<0)**(non-int) is NaN */ }
if ((n|yisint) == 0)
return (x-x)/(x-x);
s = 1.0; /* s (sign of result -ve**odd) = -1 else = 1 */
if ((n|(yisint-1)) == 0)
s = -1.0;/* (-ve)**(odd int) */
/* |y| is huge */ /* |y| is huge */
if (iy > 0x41e00000) { /* if |y| > 2**31 */ if (iy > 0x41e00000) { /* if |y| > 2**31 */

35
src/math/powf.c
View File

@ -57,17 +57,15 @@ float powf(float x, float y)
ix = hx & 0x7fffffff; ix = hx & 0x7fffffff;
iy = hy & 0x7fffffff; iy = hy & 0x7fffffff;
/* y == 0: x**0 = 1 */ /* x**0 = 1, even if x is NaN */
if (iy == 0) if (iy == 0)
return 1.0f; return 1.0f;
/* 1**y = 1, even if y is NaN */
/* x == 1: 1**y = 1, even if y is NaN */
if (hx == 0x3f800000) if (hx == 0x3f800000)
return 1.0f; return 1.0f;
/* NaN if either arg is NaN */
/* y != 0: result is NaN if either arg is NaN */
if (ix > 0x7f800000 || iy > 0x7f800000) if (ix > 0x7f800000 || iy > 0x7f800000)
return (x+0.0f) + (y+0.0f); return x + y;
/* determine if y is an odd int when x < 0 /* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer * yisint = 0 ... y is not an integer
@ -93,13 +91,10 @@ float powf(float x, float y)
else if (ix > 0x3f800000) /* (|x|>1)**+-inf = inf,0 */ else if (ix > 0x3f800000) /* (|x|>1)**+-inf = inf,0 */
return hy >= 0 ? y : 0.0f; return hy >= 0 ? y : 0.0f;
else /* (|x|<1)**+-inf = 0,inf */ else /* (|x|<1)**+-inf = 0,inf */
return hy < 0 ? -y : 0.0f; return hy >= 0 ? 0.0f: -y;
}
if (iy == 0x3f800000) { /* y is +-1 */
if (hy < 0)
return 1.0f/x;
return x;
} }
if (iy == 0x3f800000) /* y is +-1 */
return hy >= 0 ? x : 1.0f/x;
if (hy == 0x40000000) /* y is 2 */ if (hy == 0x40000000) /* y is 2 */
return x*x; return x*x;
if (hy == 0x3f000000) { /* y is 0.5 */ if (hy == 0x3f000000) { /* y is 0.5 */
@ -122,15 +117,13 @@ float powf(float x, float y)
return z; return z;
} }
n = ((uint32_t)hx>>31) - 1; sn = 1.0f; /* sign of result */
if (hx < 0) {
/* (x<0)**(non-int) is NaN */ if (yisint == 0) /* (x<0)**(non-int) is NaN */
if ((n|yisint) == 0) return (x-x)/(x-x);
return (x-x)/(x-x); if (yisint == 1) /* (x<0)**(odd int) */
sn = -1.0f;
sn = 1.0f; /* s (sign of result -ve**odd) = -1 else = 1 */ }
if ((n|(yisint-1)) == 0) /* (-ve)**(odd int) */
sn = -1.0f;
/* |y| is huge */ /* |y| is huge */
if (iy > 0x4d000000) { /* if |y| > 2**27 */ if (iy > 0x4d000000) { /* if |y| > 2**27 */