diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_logb.c | 37 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c | 13 |
2 files changed, 30 insertions, 20 deletions
diff --git a/sysdeps/ieee754/dbl-64/s_logb.c b/sysdeps/ieee754/dbl-64/s_logb.c index 2382fbb414..baa35e14d8 100644 --- a/sysdeps/ieee754/dbl-64/s_logb.c +++ b/sysdeps/ieee754/dbl-64/s_logb.c @@ -10,10 +10,6 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_logb.c,v 1.8 1995/05/10 20:47:50 jtc Exp $"; -#endif - /* * double logb(x) * IEEE 754 logb. Included to pass IEEE test suite. Not recommend. @@ -23,20 +19,29 @@ static char rcsid[] = "$NetBSD: s_logb.c,v 1.8 1995/05/10 20:47:50 jtc Exp $"; #include <math.h> #include <math_private.h> -double __logb(double x) +double +__logb (double x) { - int32_t lx,ix; - EXTRACT_WORDS(ix,lx,x); - ix &= 0x7fffffff; /* high |x| */ - if((ix|lx)==0) return -1.0/fabs(x); - if(ix>=0x7ff00000) return x*x; - if((ix>>=20)==0) /* IEEE 754 logb */ - return -1022.0; - else - return (double) (ix-1023); + int32_t lx, ix, rix; + + EXTRACT_WORDS (ix, lx, x); + ix &= 0x7fffffff; /* high |x| */ + if ((ix | lx) == 0) + return -1.0 / fabs (x); + if (ix >= 0x7ff00000) + return x * x; + if (__builtin_expect ((rix = ix >> 20) == 0, 0)) + { + /* POSIX specifies that denormal number is treated as + though it were normalized. */ + int m1 = (ix == 0) ? 0 : __builtin_clz (ix); + int m2 = (lx == 0) ? 0 : __builtin_clz (lx); + int ma = (m1 == 0) ? m2 + 32 : m1; + return -1022.0 + (double)(11 - ma); + } + return (double) (rix - 1023); } weak_alias (__logb, logb) #ifdef NO_LONG_DOUBLE -strong_alias (__logb, __logbl) -weak_alias (__logb, logbl) +strong_alias (__logb, __logbl) weak_alias (__logb, logbl) #endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c index 2ad6c7ddbd..474eeef36b 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c @@ -1,5 +1,5 @@ /* Compute radix independent exponent. - Copyright (C) 2011 Free Software Foundation, Inc. + Copyright (C) 2011, 2012 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. @@ -25,16 +25,21 @@ double __logb (double x) { - int64_t ix; + int64_t ix, ex; EXTRACT_WORDS64 (ix, x); ix &= UINT64_C(0x7fffffffffffffff); if (ix == 0) return -1.0 / fabs (x); - unsigned int ex = ix >> 52; + ex = ix >> 52; if (ex == 0x7ff) return x * x; - return ex == 0 ? -1022.0 : (double) (ex - 1023); + if (__builtin_expect (ex == 0, 0)) + { + int m = (ix == 0) ? 0 : __builtin_clzl (ix); + return -1022.0 + (double)(11 -m); + } + return (double) (ex - 1023); } weak_alias (__logb, logb) #ifdef NO_LONG_DOUBLE |