riscv
/
musl
mirror of https://git.musl-libc.org/git/musl
1
0
Fork 0
Code Issues Projects Releases Wiki Activity
Browse Source

Merge remote branch 'nsz/master'

rs-1.0
Rich Felker 14 years ago
parent
ad47d45e9d 91c28f61f4
commit
58bf74850f
10 changed files with 127 additions and 198 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 2
      src/math/__signbit.c
  2. 2
      src/math/__signbitf.c
  3. 31
      src/math/cbrtl.c
  4. 23
      src/math/nearbyint.c
  5. 14
      src/math/nearbyintf.c
  6. 11
      src/math/nearbyintl.c
  7. 63
      src/math/pow.c
  8. 31
      src/math/powf.c
  9. 129
      src/math/powl.c
  10. 1
      src/math/remquo.c

2
src/math/__signbit.c
View File

@ -1,6 +1,6 @@
#include "libm.h" #include "libm.h"
// FIXME: macro // FIXME: macro in math.h
int __signbit(double x) int __signbit(double x)
{ {
union { union {

2
src/math/__signbitf.c
View File

@ -1,6 +1,6 @@
#include "libm.h" #include "libm.h"
// FIXME // FIXME: macro in math.h
int __signbitf(float x) int __signbitf(float x)
{ {
union { union {

31
src/math/cbrtl.c
View File

@ -23,9 +23,9 @@ long double cbrtl(long double x)
return cbrt(x); return cbrt(x);
} }
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#define BIAS (LDBL_MAX_EXP - 1) #define BIAS (LDBL_MAX_EXP - 1)
static const unsigned static const unsigned B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */
B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */
long double cbrtl(long double x) long double cbrtl(long double x)
{ {
@ -48,25 +48,10 @@ long double cbrtl(long double x)
if (k == BIAS + LDBL_MAX_EXP) if (k == BIAS + LDBL_MAX_EXP)
return x + x; return x + x;
// FIXME: extended precision is default on linux..
#undef __i386__
#ifdef __i386__
fp_prec_t oprec;
oprec = fpgetprec();
if (oprec != FP_PE)
fpsetprec(FP_PE);
#endif
if (k == 0) { if (k == 0) {
/* If x = +-0, then cbrt(x) = +-0. */ /* If x = +-0, then cbrt(x) = +-0. */
if ((u.bits.manh | u.bits.manl) == 0) { if ((u.bits.manh | u.bits.manl) == 0)
#ifdef __i386__ return x;
if (oprec != FP_PE)
fpsetprec(oprec);
#endif
return (x);
}
/* Adjust subnormal numbers. */ /* Adjust subnormal numbers. */
u.e *= 0x1.0p514; u.e *= 0x1.0p514;
k = u.bits.exp; k = u.bits.exp;
@ -118,7 +103,7 @@ long double cbrtl(long double x)
* Round it away from zero to 32 bits (32 so that t*t is exact, and * Round it away from zero to 32 bits (32 so that t*t is exact, and
* away from zero for technical reasons). * away from zero for technical reasons).
*/ */
t = dt + (0x1.0p32L + 0x1.0p-32L) - 0x1.0p32; t = dt + (0x1.0p32L + 0x1.0p-31L) - 0x1.0p32;
#elif LDBL_MANT_DIG == 113 #elif LDBL_MANT_DIG == 113
/* /*
* Round dt away from zero to 47 bits. Since we don't trust the 47, * Round dt away from zero to 47 bits. Since we don't trust the 47,
@ -129,8 +114,6 @@ long double cbrtl(long double x)
* dt is much more accurate than needed. * dt is much more accurate than needed.
*/ */
t = dt + 0x2.0p-46 + 0x1.0p60L - 0x1.0p60; t = dt + 0x2.0p-46 + 0x1.0p60L - 0x1.0p60;
#else
#error "Unsupported long double format"
#endif #endif
/* /*
@ -144,10 +127,6 @@ long double cbrtl(long double x)
t = t+t*r; /* error <= 0.5 + 0.5/3 + epsilon */ t = t+t*r; /* error <= 0.5 + 0.5/3 + epsilon */
t *= v.e; t *= v.e;
#ifdef __i386__
if (oprec != FP_PE)
fpsetprec(oprec);
#endif
return t; return t;
} }
#endif #endif

23
src/math/nearbyint.c
View File

@ -1,20 +1,19 @@
#include <fenv.h> #include <fenv.h>
#include <math.h> #include <math.h>
/* /* nearbyint is the same as rint, but it must not raise the inexact exception */
rint may raise inexact (and it should not alter the fenv otherwise)
nearbyint must not raise inexact
(according to ieee754r section 7.9 both functions should raise invalid double nearbyint(double x)
when the input is signaling nan, but c99 does not define snan so saving {
and restoring the entire fenv should be fine) #ifdef FE_INEXACT
*/ int e;
double nearbyint(double x) { e = fetestexcept(FE_INEXACT);
fenv_t e; #endif
fegetenv(&e);
x = rint(x); x = rint(x);
fesetenv(&e); #ifdef FE_INEXACT
if (!e)
feclearexcept(FE_INEXACT);
#endif
return x; return x;
} }

14
src/math/nearbyintf.c
View File

@ -1,11 +1,17 @@
#include <fenv.h> #include <fenv.h>
#include <math.h> #include <math.h>
float nearbyintf(float x) { float nearbyintf(float x)
fenv_t e; {
#ifdef FE_INEXACT
int e;
fegetenv(&e); e = fetestexcept(FE_INEXACT);
#endif
x = rintf(x); x = rintf(x);
fesetenv(&e); #ifdef FE_INEXACT
if (!e)
feclearexcept(FE_INEXACT);
#endif
return x; return x;
} }

11
src/math/nearbyintl.c
View File

@ -10,11 +10,16 @@ long double nearbyintl(long double x)
#include <fenv.h> #include <fenv.h>
long double nearbyintl(long double x) long double nearbyintl(long double x)
{ {
fenv_t e; #ifdef FE_INEXACT
int e;
fegetenv(&e); e = fetestexcept(FE_INEXACT);
#endif
x = rintl(x); x = rintl(x);
fesetenv(&e); #ifdef FE_INEXACT
if (!e)
feclearexcept(FE_INEXACT);
#endif
return x; return x;
} }
#endif #endif

63
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. */
n = ((uint32_t)hx>>31) - 1;
/* (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) */
s = 1.0; /* s (sign of result -ve**odd) = -1 else = 1 */ s = -1.0;
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 */

31
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; /* s (sign of result -ve**odd) = -1 else = 1 */
if ((n|(yisint-1)) == 0) /* (-ve)**(odd int) */
sn = -1.0f; sn = -1.0f;
}
/* |y| is huge */ /* |y| is huge */
if (iy > 0x4d000000) { /* if |y| > 2**27 */ if (iy > 0x4d000000) { /* if |y| > 2**27 */

129
src/math/powl.c
View File

@ -78,8 +78,6 @@ long double powl(long double x, long double y)
/* Table size */ /* Table size */
#define NXT 32 #define NXT 32
/* log2(Table size) */
#define LNXT 5
/* log(1+x) = x - .5x^2 + x^3 * P(z)/Q(z) /* log(1+x) = x - .5x^2 + x^3 * P(z)/Q(z)
* on the domain 2^(-1/32) - 1 <= x <= 2^(1/32) - 1 * on the domain 2^(-1/32) - 1 <= x <= 2^(1/32) - 1
@ -203,38 +201,35 @@ long double powl(long double x, long double y)
volatile long double z=0; volatile long double z=0;
long double w=0, W=0, Wa=0, Wb=0, ya=0, yb=0, u=0; long double w=0, W=0, Wa=0, Wb=0, ya=0, yb=0, u=0;
if (y == 0.0) /* make sure no invalid exception is raised by nan comparision */
if (isnan(x)) {
if (!isnan(y) && y == 0.0)
return 1.0; return 1.0;
if (isnan(x))
return x; return x;
if (isnan(y)) }
if (isnan(y)) {
if (x == 1.0)
return 1.0;
return y; return y;
}
if (x == 1.0)
return 1.0; /* 1**y = 1, even if y is nan */
if (x == -1.0 && !isfinite(y))
return 1.0; /* -1**inf = 1 */
if (y == 0.0)
return 1.0; /* x**0 = 1, even if x is nan */
if (y == 1.0) if (y == 1.0)
return x; return x;
// FIXME: this is wrong, see pow special cases in c99 F.9.4.4
if (!isfinite(y) && (x == -1.0 || x == 1.0) )
return y - y; /* +-1**inf is NaN */
if (x == 1.0)
return 1.0;
if (y >= LDBL_MAX) { if (y >= LDBL_MAX) {
if (x > 1.0) if (x > 1.0 || x < -1.0)
return INFINITY;
if (x > 0.0 && x < 1.0)
return 0.0;
if (x < -1.0)
return INFINITY; return INFINITY;
if (x > -1.0 && x < 0.0) if (x != 0.0)
return 0.0; return 0.0;
} }
if (y <= -LDBL_MAX) { if (y <= -LDBL_MAX) {
if (x > 1.0) if (x > 1.0 || x < -1.0)
return 0.0; return 0.0;
if (x > 0.0 && x < 1.0) if (x != 0.0)
return INFINITY;
if (x < -1.0)
return 0.0;
if (x > -1.0 && x < 0.0)
return INFINITY; return INFINITY;
} }
if (x >= LDBL_MAX) { if (x >= LDBL_MAX) {
@ -244,6 +239,7 @@ long double powl(long double x, long double y)
} }
w = floorl(y); w = floorl(y);
/* Set iyflg to 1 if y is an integer. */ /* Set iyflg to 1 if y is an integer. */
iyflg = 0; iyflg = 0;
if (w == y) if (w == y)
@ -271,42 +267,32 @@ long double powl(long double x, long double y)
return 0.0; return 0.0;
} }
} }
nflg = 0; /* (x<0)**(odd int) */
nflg = 0; /* flag = 1 if x<0 raised to integer power */
if (x <= 0.0) { if (x <= 0.0) {
if (x == 0.0) { if (x == 0.0) {
if (y < 0.0) { if (y < 0.0) {
if (signbit(x) && yoddint) if (signbit(x) && yoddint)
return -INFINITY; /* (-0.0)**(-odd int) = -inf, divbyzero */
return INFINITY; return -1.0/0.0;
/* (+-0.0)**(negative) = inf, divbyzero */
return 1.0/0.0;
} }
if (y > 0.0) {
if (signbit(x) && yoddint) if (signbit(x) && yoddint)
return -0.0; return -0.0;
return 0.0; return 0.0;
} }
if (y == 0.0)
return 1.0; /* 0**0 */
return 0.0; /* 0**y */
}
if (iyflg == 0) if (iyflg == 0)
return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */ return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */
nflg = 1; /* (x<0)**(integer) */
if (yoddint)
nflg = 1; /* negate result */
x = -x;
} }
/* (+integer)**(integer) */
/* Integer power of an integer. */ if (iyflg && floorl(x) == x && fabsl(y) < 32768.0) {
if (iyflg) {
i = w;
w = floorl(x);
if (w == x && fabsl(y) < 32768.0) {
w = powil(x, (int)y); w = powil(x, (int)y);
return w; return nflg ? -w : w;
} }
}
if (nflg)
x = fabsl(x);
/* separate significand from exponent */ /* separate significand from exponent */
x = frexpl(x, &i); x = frexpl(x, &i);
@ -354,9 +340,7 @@ long double powl(long double x, long double y)
z += x; z += x;
/* Compute exponent term of the base 2 logarithm. */ /* Compute exponent term of the base 2 logarithm. */
w = -i; w = -i / NXT;
// TODO: use w * 0x1p-5;
w = scalbnl(w, -LNXT); /* divide by NXT */
w += e; w += e;
/* Now base 2 log of x is w + z. */ /* Now base 2 log of x is w + z. */
@ -381,7 +365,7 @@ long double powl(long double x, long double y)
H = Fb + Gb; H = Fb + Gb;
Ha = reducl(H); Ha = reducl(H);
w = scalbnl( Ga+Ha, LNXT ); w = (Ga + Ha) * NXT;
/* Test the power of 2 for overflow */ /* Test the power of 2 for overflow */
if (w > MEXP) if (w > MEXP)
@ -418,18 +402,8 @@ long double powl(long double x, long double y)
z = z + w; z = z + w;
z = scalbnl(z, i); /* multiply by integer power of 2 */ z = scalbnl(z, i); /* multiply by integer power of 2 */
if (nflg) { if (nflg)
/* For negative x, z = -z;
* find out if the integer exponent
* is odd or even.
*/
w = 0.5*y;
w = floorl(w);
w = 2.0*w;
if (w != y)
z = -z; /* odd exponent */
}
return z; return z;
} }
@ -439,15 +413,14 @@ static long double reducl(long double x)
{ {
long double t; long double t;
t = scalbnl(x, LNXT); t = x * NXT;
t = floorl(t); t = floorl(t);
t = scalbnl(t, -LNXT); t = t / NXT;
return t; return t;
} }
/* powil.c /*
* * Positive real raised to integer power, long double precision
* Real raised to integer power, long double precision
* *
* *
* SYNOPSIS: * SYNOPSIS:
@ -460,7 +433,7 @@ static long double reducl(long double x)
* *
* DESCRIPTION: * DESCRIPTION:
* *
* Returns argument x raised to the nth power. * Returns argument x>0 raised to the nth power.
* The routine efficiently decomposes n as a sum of powers of * The routine efficiently decomposes n as a sum of powers of
* two. The desired power is a product of two-to-the-kth * two. The desired power is a product of two-to-the-kth
* powers of x. Thus to compute the 32767 power of x requires * powers of x. Thus to compute the 32767 power of x requires
@ -482,24 +455,10 @@ static long double powil(long double x, int nn)
{ {
long double ww, y; long double ww, y;
long double s; long double s;
int n, e, sign, asign, lx; int n, e, sign, lx;
if (x == 0.0) {
if (nn == 0) if (nn == 0)
return 1.0; return 1.0;
else if (nn < 0)
return LDBL_MAX;
return 0.0;
}
if (nn == 0)
return 1.0;
if (x < 0.0) {
asign = -1;
x = -x;
} else
asign = 0;
if (nn < 0) { if (nn < 0) {
sign = -1; sign = -1;
@ -539,10 +498,8 @@ static long double powil(long double x, int nn)
/* First bit of the power */ /* First bit of the power */
if (n & 1) if (n & 1)
y = x; y = x;
else { else
y = 1.0; y = 1.0;
asign = 0;
}
ww = x; ww = x;
n >>= 1; n >>= 1;
@ -553,8 +510,6 @@ static long double powil(long double x, int nn)
n >>= 1; n >>= 1;
} }
if (asign)
y = -y; /* odd power of negative number */
if (sign < 0) if (sign < 0)
y = 1.0/y; y = 1.0/y;
return y; return y;

1
src/math/remquo.c
View File

@ -35,7 +35,6 @@ double remquo(double x, double y, int *quo)
hy &= 0x7fffffff; /* |y| */ hy &= 0x7fffffff; /* |y| */
/* purge off exception values */ /* purge off exception values */
// FIXME: signed shift
if ((hy|ly) == 0 || hx >= 0x7ff00000 || /* y = 0, or x not finite */ if ((hy|ly) == 0 || hx >= 0x7ff00000 || /* y = 0, or x not finite */
(hy|((ly|-ly)>>31)) > 0x7ff00000) /* or y is NaN */ (hy|((ly|-ly)>>31)) > 0x7ff00000) /* or y is NaN */
return (x*y)/(x*y); return (x*y)/(x*y);

Write Preview
Loading…
Cancel
Save