diff options
Diffstat (limited to 'libm')
| -rw-r--r-- | libm/e_lgamma_r.c | 42 | ||||
| -rw-r--r-- | libm/e_pow.c | 63 | ||||
| -rw-r--r-- | libm/e_scalb.c | 13 | ||||
| -rw-r--r-- | libm/float_wrappers.c | 3 | ||||
| -rw-r--r-- | libm/ldouble_wrappers.c | 207 | ||||
| -rw-r--r-- | libm/s_ilogb.c | 47 | ||||
| -rw-r--r-- | libm/s_modf.c | 15 | ||||
| -rw-r--r-- | libm/s_rint.c | 54 |
8 files changed, 177 insertions, 267 deletions
diff --git a/libm/e_lgamma_r.c b/libm/e_lgamma_r.c index 6d5d9c514..dbb0a0d5d 100644 --- a/libm/e_lgamma_r.c +++ b/libm/e_lgamma_r.c @@ -360,31 +360,35 @@ strong_alias(__ieee754_lgamma, gamma) #endif -/* FIXME! Looks like someone just used __ieee754_gamma_r, - * believing it's a "true" gamma function, but it was not! - * Our tgamma is WRONG. - */ - /* double tgamma(double x) * Return the Gamma function of x. */ double tgamma(double x) { - double y; - int local_signgam; + int sign_of_gamma; + int32_t hx; + u_int32_t lx; - y = __ieee754_lgamma_r(x, &local_signgam); // was __ieee754_gamma_r - if (local_signgam < 0) - y = -y; -#ifndef _IEEE_LIBM - if (_LIB_VERSION == _IEEE_) - return y; - if (!isfinite(y) && isfinite(x)) { - if (floor(x) == x && x <= 0.0) - return __kernel_standard(x, x, 41); /* tgamma pole */ - return __kernel_standard(x, x, 40); /* tgamma overflow */ + /* We don't have a real gamma implementation now. We'll use lgamma + and the exp function. But due to the required boundary + conditions we must check some values separately. */ + + EXTRACT_WORDS(hx, lx, x); + + if (((hx & 0x7fffffff) | lx) == 0) { + /* Return value for x == 0 is Inf with divide by zero exception. */ + return 1.0 / x; } -#endif - return y; + if (hx < 0 && (u_int32_t)hx < 0xfff00000 && rint(x) == x) { + /* Return value for integer x < 0 is NaN with invalid exception. */ + return (x - x) / (x - x); + } + if ((u_int32_t)hx == 0xfff00000 && lx == 0) { + /* x == -Inf. According to ISO this is NaN. */ + return x - x; + } + + x = exp(lgamma_r(x, &sign_of_gamma)); + return sign_of_gamma >= 0 ? x : -x; } libm_hidden_def(tgamma) diff --git a/libm/e_pow.c b/libm/e_pow.c index 137f600c3..3be900389 100644 --- a/libm/e_pow.c +++ b/libm/e_pow.c @@ -21,25 +21,26 @@ * 3. Return x**y = 2**n*exp(y'*log2) * * Special cases: - * 1. (anything) ** 0 is 1 - * 2. (anything) ** 1 is itself - * 3. (anything) ** NAN is NAN - * 4. NAN ** (anything except 0) is NAN - * 5. +-(|x| > 1) ** +INF is +INF - * 6. +-(|x| > 1) ** -INF is +0 - * 7. +-(|x| < 1) ** +INF is +0 - * 8. +-(|x| < 1) ** -INF is +INF - * 9. +-1 ** +-INF is NAN - * 10. +0 ** (+anything except 0, NAN) is +0 - * 11. -0 ** (+anything except 0, NAN, odd integer) is +0 - * 12. +0 ** (-anything except 0, NAN) is +INF - * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF - * 14. -0 ** (odd integer) = -( +0 ** (odd integer) ) - * 15. +INF ** (+anything except 0,NAN) is +INF - * 16. +INF ** (-anything except 0,NAN) is +0 - * 17. -INF ** (anything) = -0 ** (-anything) - * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) - * 19. (-anything except 0 and inf) ** (non-integer) is NAN + * 1. +-1 ** anything is 1.0 + * 2. +-1 ** +-INF is 1.0 + * 3. (anything) ** 0 is 1 + * 4. (anything) ** 1 is itself + * 5. (anything) ** NAN is NAN + * 6. NAN ** (anything except 0) is NAN + * 7. +-(|x| > 1) ** +INF is +INF + * 8. +-(|x| > 1) ** -INF is +0 + * 9. +-(|x| < 1) ** +INF is +0 + * 10 +-(|x| < 1) ** -INF is +INF + * 11. +0 ** (+anything except 0, NAN) is +0 + * 12. -0 ** (+anything except 0, NAN, odd integer) is +0 + * 13. +0 ** (-anything except 0, NAN) is +INF + * 14. -0 ** (-anything except 0, NAN, odd integer) is +INF + * 15. -0 ** (odd integer) = -( +0 ** (odd integer) ) + * 16. +INF ** (+anything except 0,NAN) is +INF + * 17. +INF ** (-anything except 0,NAN) is +0 + * 18. -INF ** (anything) = -0 ** (-anything) + * 19. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) + * 20. (-anything except 0 and inf) ** (non-integer) is NAN * * Accuracy: * pow(x,y) returns x**y nearly rounded. In particular @@ -99,8 +100,14 @@ double attribute_hidden __ieee754_pow(double x, double y) u_int32_t lx,ly; EXTRACT_WORDS(hx,lx,x); + /* x==1: 1**y = 1 (even if y is NaN) */ + if (hx==0x3ff00000 && lx==0) { + return x; + } + ix = hx&0x7fffffff; + EXTRACT_WORDS(hy,ly,y); - ix = hx&0x7fffffff; iy = hy&0x7fffffff; + iy = hy&0x7fffffff; /* y==zero: x**0 = 1 */ if((iy|ly)==0) return one; @@ -132,13 +139,13 @@ double attribute_hidden __ieee754_pow(double x, double y) /* special value of y */ if(ly==0) { - if (iy==0x7ff00000) { /* y is +-inf */ - if(((ix-0x3ff00000)|lx)==0) - return y - y; /* inf**+-1 is NaN */ - else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */ - return (hy>=0)? y: zero; - else /* (|x|<1)**-,+inf = inf,0 */ - return (hy<0)?-y: zero; + if (iy==0x7ff00000) { /* y is +-inf */ + if (((ix-0x3ff00000)|lx)==0) + return one; /* +-1**+-inf is 1 (yes, weird rule) */ + if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */ + return (hy>=0) ? y : zero; + /* (|x|<1)**-,+inf = inf,0 */ + return (hy<0) ? -y : zero; } if(iy==0x3ff00000) { /* y is +-1 */ if(hy<0) return one/x; else return x; @@ -146,7 +153,7 @@ double attribute_hidden __ieee754_pow(double x, double y) if(hy==0x40000000) return x*x; /* y is 2 */ if(hy==0x3fe00000) { /* y is 0.5 */ if(hx>=0) /* x >= +0 */ - return __ieee754_sqrt(x); + return __ieee754_sqrt(x); } } diff --git a/libm/e_scalb.c b/libm/e_scalb.c index 9244cf42f..3fd78a1c0 100644 --- a/libm/e_scalb.c +++ b/libm/e_scalb.c @@ -10,9 +10,9 @@ */ /* - * __ieee754_scalb(x, fn) is provide for - * passing various standard test suite. One - * should use scalbn() instead. + * __ieee754_scalb(x, fn) is provided for + * passing various standard test suites. + * One should use scalbn() instead. */ #include "math.h" @@ -21,7 +21,6 @@ double attribute_hidden __ieee754_scalb(double x, double fn) { - return scalbn(x,fn); if (isnan(x)||isnan(fn)) return x*fn; if (!isfinite(fn)) { if(fn>0.0) return x*fn; @@ -35,9 +34,9 @@ double attribute_hidden __ieee754_scalb(double x, double fn) #if defined __UCLIBC_SUSV3_LEGACY__ /* - * wrapper scalb(double x, double fn) is provide for - * passing various standard test suite. One - * should use scalbn() instead. + * wrapper scalb(double x, double fn) is provided for + * passing various standard test suites. + * One should use scalbn() instead. */ #ifndef _IEEE_LIBM double scalb(double x, double fn) diff --git a/libm/float_wrappers.c b/libm/float_wrappers.c index a13aac93a..82b7963e1 100644 --- a/libm/float_wrappers.c +++ b/libm/float_wrappers.c @@ -293,10 +293,9 @@ long_WRAPPER1(lround) float modff (float x, float *iptr) { double y, result; - result = modf ( x, &y ); + result = modf( x, &y ); *iptr = (float)y; return (float) result; - } #endif diff --git a/libm/ldouble_wrappers.c b/libm/ldouble_wrappers.c index bf7ae15f0..118a78f64 100644 --- a/libm/ldouble_wrappers.c +++ b/libm/ldouble_wrappers.c @@ -21,6 +21,11 @@ long double func##l(long double x) \ { \ return (long double) func((double) x); \ } +#define WRAPPER2(func) \ +long double func##l(long double x, long double y) \ +{ \ + return (long double) func((double) x, (double) y); \ +} #define int_WRAPPER1(func) \ int func##l(long double x) \ { \ @@ -37,96 +42,6 @@ long long func##l(long double x) \ return func((double) x); \ } -#if defined __i386__ && defined __OPTIMIZE__ -# undef WRAPPER1 -# undef int_WRAPPER1 -# undef long_WRAPPER1 -# undef long_long_WRAPPER1 -/* gcc 4.3.1 generates really ugly code with redundant pair of store/load: - * sub $0x10,%esp - * fldt 0x14(%esp) - * fstpl 0x8(%esp) - * fldl 0x8(%esp) <-- ?? - * fstpl (%esp) <-- ?? - * call function - * add $0x10,%esp - * ret - * I can hope newer gcc will eliminate that. However, I don't think - * it will be smart enough to reuse argument stack space and use - * jump instead of call. Let's do it by hand. - * The asm below loads long double x into st(0), then stores it back - * to the same location, but as a double. At this point, stack looks - * exactly as "double func(double)" expects it to be. - * The return value is returned in st(0) per ABI in both cases (returning - * a long double or returning a double). So we can simply jump to func. - * Using __GI_func in jump to make optimized intra-library jump. - * gcc will still generate a useless "ret" after asm. Oh well... - */ -# define WRAPPER1(func) \ -long double func##l(long double x) \ -{ \ - long double st_top; \ - __asm__ ( \ - " fldt %1\n" \ - " fstpl %1\n" \ - " jmp " __stringify(__GI_##func) "\n" \ - : "=t" (st_top) \ - : "m" (x) \ - ); \ - return st_top; \ -} -# define int_WRAPPER1(func) \ -int func##l(long double x) \ -{ \ - int ret; \ - __asm__ ( \ - " fldt %1\n" \ - " fstpl %1\n" \ - " jmp " __stringify(__GI_##func) "\n" \ - : "=a" (ret) \ - : "m" (x) \ - ); \ - return ret; \ -} -# define long_WRAPPER1(func) \ -long func##l(long double x) \ -{ \ - long ret; \ - __asm__ ( \ - " fldt %1\n" \ - " fstpl %1\n" \ - " jmp " __stringify(__GI_##func) "\n" \ - : "=a" (ret) \ - : "m" (x) \ - ); \ - return ret; \ -} -# define long_long_WRAPPER1(func) \ -long long func##l(long double x) \ -{ \ - long long ret; \ - __asm__ ( \ - " fldt %1\n" \ - " fstpl %1\n" \ - " jmp " __stringify(__GI_##func) "\n" \ - : "=A" (ret) \ - : "m" (x) \ - ); \ - return ret; \ -} -#endif /* __i386__ && __OPTIMIZE__ */ - -#if defined __NO_LONG_DOUBLE_MATH -# define int_WRAPPER_C99(func) /* not needed */ -# else -# define int_WRAPPER_C99(func) \ -int func##l(long double x) \ -{ \ - return func((double) x); \ -} \ -libm_hidden_def(func##l) -#endif - /* Implement the following, as defined by SuSv3 */ #if 0 long double acoshl(long double); @@ -205,10 +120,7 @@ WRAPPER1(asin) #endif #ifdef L_atan2l -long double atan2l (long double x, long double y) -{ - return (long double) atan2( (double)x, (double)y ); -} +WRAPPER2(atan2) #endif #ifdef L_atanhl @@ -235,10 +147,7 @@ WRAPPER1(ceil) #endif #ifdef L_copysignl -long double copysignl (long double x, long double y) -{ - return (long double) copysign( (double)x, (double)y ); -} +WRAPPER2(copysign) #endif #ifdef L_coshl @@ -274,10 +183,7 @@ WRAPPER1(fabs) #endif #ifdef L_fdiml -long double fdiml (long double x, long double y) -{ - return (long double) fdim( (double)x, (double)y ); -} +WRAPPER2(fdim) #endif #ifdef L_floorl @@ -292,24 +198,15 @@ long double fmal (long double x, long double y, long double z) #endif #ifdef L_fmaxl -long double fmaxl (long double x, long double y) -{ - return (long double) fmax( (double)x, (double)y ); -} +WRAPPER2(fmax) #endif #ifdef L_fminl -long double fminl (long double x, long double y) -{ - return (long double) fmin( (double)x, (double)y ); -} +WRAPPER2(fmin) #endif #ifdef L_fmodl -long double fmodl (long double x, long double y) -{ - return (long double) fmod( (double)x, (double)y ); -} +WRAPPER2(fmod) #endif #ifdef L_frexpl @@ -320,19 +217,11 @@ long double frexpl (long double x, int *ex) #endif #ifdef L_gammal -/* WRAPPER1(gamma) won't work, tries to call __GI_xxx, - * and gamma() hasn't got one. */ -long double gammal(long double x) -{ - return (long double) gamma((double) x); -} +WRAPPER1(gamma) #endif #ifdef L_hypotl -long double hypotl (long double x, long double y) -{ - return (long double) hypot( (double)x, (double)y ); -} +WRAPPER2(hypot) #endif #ifdef L_ilogbl @@ -367,11 +256,7 @@ WRAPPER1(log1p) #endif #ifdef L_log2l -/* WRAPPER1(log2) won't work */ -long double log2l(long double x) -{ - return (long double) log2((double)x); -} +WRAPPER1(log2) #endif #ifdef L_logbl @@ -405,33 +290,21 @@ WRAPPER1(nearbyint) #endif #ifdef L_nextafterl -long double nextafterl (long double x, long double y) -{ - return (long double) nextafter( (double)x, (double)y ); -} +WRAPPER2(nextafter) #endif /* Disabled in Makefile.in */ #if 0 /* def L_nexttowardl */ -long double nexttowardl (long double x, long double y) -{ - return (long double) nexttoward( (double)x, (double)y ); -} +WRAPPER2(nexttoward) libm_hidden_def(nexttowardl) #endif #ifdef L_powl -long double powl (long double x, long double y) -{ - return (long double) pow( (double)x, (double)y ); -} +WRAPPER2(pow) #endif #ifdef L_remainderl -long double remainderl (long double x, long double y) -{ - return (long double) remainder( (double)x, (double)y ); -} +WRAPPER2(remainder) #endif #ifdef L_remquol @@ -494,34 +367,34 @@ WRAPPER1(trunc) #endif #ifdef L_significandl -/* WRAPPER1(significand) won't work, tries to call __GI_xxx, - * and significand() hasn't got one. */ -long double significandl(long double x) -{ - return (long double) significand((double) x); -} +WRAPPER1(significand) #endif -#ifdef __DO_C99_MATH__ +#if defined __DO_C99_MATH__ && !defined __NO_LONG_DOUBLE_MATH -#ifdef L___fpclassifyl -int_WRAPPER_C99(__fpclassify) -#endif +# ifdef L___fpclassifyl +int_WRAPPER1(__fpclassify) +libm_hidden_def(__fpclassifyl) +# endif -#ifdef L___finitel -int_WRAPPER_C99(__finite) -#endif +# ifdef L___finitel +int_WRAPPER1(__finite) +libm_hidden_def(__finitel) +# endif -#ifdef L___signbitl -int_WRAPPER_C99(__signbit) -#endif +# ifdef L___signbitl +int_WRAPPER1(__signbit) +libm_hidden_def(__signbitl) +# endif -#ifdef L___isnanl -int_WRAPPER_C99(__isnan) -#endif +# ifdef L___isnanl +int_WRAPPER1(__isnan) +libm_hidden_def(__isnanl) +# endif -#ifdef L___isinfl -int_WRAPPER_C99(__isinf) -#endif +# ifdef L___isinfl +int_WRAPPER1(__isinf) +libm_hidden_def(__isinfl) +# endif #endif diff --git a/libm/s_ilogb.c b/libm/s_ilogb.c index ce3029667..259ae7b55 100644 --- a/libm/s_ilogb.c +++ b/libm/s_ilogb.c @@ -10,9 +10,10 @@ */ /* ilogb(double x) - * return the binary exponent of non-zero x - * ilogb(0) = 0x80000001 - * ilogb(inf/NaN) = 0x7fffffff (no signal is raised) + * return the binary exponent of x + * ilogb(+-0) = FP_ILOGB0 + * ilogb(+-inf) = INT_MAX + * ilogb(NaN) = FP_ILOGBNAN (no signal is raised) */ #include "math.h" @@ -22,21 +23,35 @@ int ilogb(double x) { int32_t hx,lx,ix; - GET_HIGH_WORD(hx,x); + GET_HIGH_WORD(hx, x); hx &= 0x7fffffff; - if(hx<0x00100000) { - GET_LOW_WORD(lx,x); - if((hx|lx)==0) - return 0x80000001; /* ilogb(0) = 0x80000001 */ - else /* subnormal x */ - if(hx==0) { - for (ix = -1043; lx>0; lx<<=1) ix -=1; - } else { - for (ix = -1022,hx<<=11; hx>0; hx<<=1) ix -=1; + + if (hx < 0x00100000) { + GET_LOW_WORD(lx, x); + if ((hx|lx)==0) /* +-0, ilogb(0) = FP_ILOGB0 */ + return FP_ILOGB0; + /* subnormal x */ + ix = -1043; + if (hx != 0) { + ix = -1022; + lx = (hx << 11); } - return ix; + /* each leading zero mantissa bit makes exponent smaller */ + for (; lx > 0; lx <<= 1) + ix--; + return ix; } - else if (hx<0x7ff00000) return (hx>>20)-1023; - else return 0x7fffffff; + + if (hx < 0x7ff00000) /* normal x */ + return (hx>>20) - 1023; + + if (FP_ILOGBNAN != (~0U >> 1)) { + GET_LOW_WORD(lx, x); + if (hx == 0x7ff00000 && lx == 0) /* +-inf */ + return ~0U >> 1; /* = INT_MAX */ + } + + /* NAN. ilogb(NAN) = FP_ILOGBNAN */ + return FP_ILOGBNAN; } libm_hidden_def(ilogb) diff --git a/libm/s_modf.c b/libm/s_modf.c index 26c2f1cd5..f73d5fddb 100644 --- a/libm/s_modf.c +++ b/libm/s_modf.c @@ -37,10 +37,8 @@ double modf(double x, double *iptr) } else { i = (0x000fffff)>>j0; if(((i0&i)|i1)==0) { /* x is integral */ - u_int32_t high; *iptr = x; - GET_HIGH_WORD(high,x); - INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ + INSERT_WORDS(x,i0&0x80000000,0); /* return +-0 */ return x; } else { INSERT_WORDS(*iptr,i0&(~i),0); @@ -48,18 +46,17 @@ double modf(double x, double *iptr) } } } else if (j0>51) { /* no fraction part */ - u_int32_t high; *iptr = x*one; - GET_HIGH_WORD(high,x); - INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ + /* We must handle NaNs separately. */ + if (j0 == 0x400 && ((i0 & 0xfffff) | i1)) + return x*one; + INSERT_WORDS(x,i0&0x80000000,0); /* return +-0 */ return x; } else { /* fraction part in low x */ i = ((u_int32_t)(0xffffffff))>>(j0-20); if((i1&i)==0) { /* x is integral */ - u_int32_t high; *iptr = x; - GET_HIGH_WORD(high,x); - INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ + INSERT_WORDS(x,i0&0x80000000,0); /* return +-0 */ return x; } else { INSERT_WORDS(*iptr,i0,i1&(~i)); diff --git a/libm/s_rint.c b/libm/s_rint.c index 358ce76b4..06432c622 100644 --- a/libm/s_rint.c +++ b/libm/s_rint.c @@ -30,41 +30,57 @@ TWO52[2]={ double rint(double x) { - int32_t i0,j0,sx; + int32_t i0, j0, sx; u_int32_t i,i1; - double w,t; + double t; + /* We use w = x + 2^52; t = w - 2^52; trick to round x to integer. + * This trick requires that compiler does not optimize it + * by keeping intermediate result w in a register wider than double. + * Declaring w volatile assures that value gets truncated to double + * (unfortunately, it also forces store+load): + */ + volatile double w; + EXTRACT_WORDS(i0,i1,x); - sx = (i0>>31)&1; - j0 = ((i0>>20)&0x7ff)-0x3ff; - if(j0<20) { - if(j0<0) { - if(((i0&0x7fffffff)|i1)==0) return x; + /* Unbiased exponent */ + j0 = ((((u_int32_t)i0) >> 20)&0x7ff)-0x3ff; + + if (j0 > 51) { + //Why bother? Just returning x works too + //if (j0 == 0x400) /* inf or NaN */ + // return x+x; + return x; /* x is integral */ + } + + /* Sign */ + sx = ((u_int32_t)i0) >> 31; + + if (j0<20) { + if (j0<0) { /* |x| < 1 */ + if (((i0&0x7fffffff)|i1)==0) return x; i1 |= (i0&0x0fffff); i0 &= 0xfffe0000; i0 |= ((i1|-i1)>>12)&0x80000; SET_HIGH_WORD(x,i0); - w = TWO52[sx]+x; - t = w-TWO52[sx]; + w = TWO52[sx]+x; + t = w-TWO52[sx]; GET_HIGH_WORD(i0,t); SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31)); - return t; + return t; } else { i = (0x000fffff)>>j0; - if(((i0&i)|i1)==0) return x; /* x is integral */ + if (((i0&i)|i1)==0) return x; /* x is integral */ i>>=1; - if(((i0&i)|i1)!=0) { - if(j0==19) i1 = 0x40000000; else - i0 = (i0&(~i))|((0x20000)>>j0); + if (((i0&i)|i1)!=0) { + if (j0==19) i1 = 0x40000000; + else i0 = (i0&(~i))|((0x20000)>>j0); } } - } else if (j0>51) { - if(j0==0x400) return x+x; /* inf or NaN */ - else return x; /* x is integral */ } else { i = ((u_int32_t)(0xffffffff))>>(j0-20); - if((i1&i)==0) return x; /* x is integral */ + if ((i1&i)==0) return x; /* x is integral */ i>>=1; - if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20)); + if ((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20)); } INSERT_WORDS(x,i0,i1); w = TWO52[sx]+x; |