From 89c9aa491a7cee97bf78a29cddbf0a25c902a671 Mon Sep 17 00:00:00 2001 From: Adhemerval Zanella Date: Thu, 10 May 2012 15:11:55 -0500 Subject: Fix for logb/logbf/logbl (bugs 13954/13955/13956) POSIX 2008 states that if the input for 'logb[f|l]' is a subnormal number it should be treated as if it were normalized. This means the implementation should calculate the log2 of the mantissa and add it to the subnormal exponent (-126 for float and -1022 for double and IBM long double). This patch takes care of that. --- sysdeps/ieee754/dbl-64/s_logb.c | 37 ++++++++++++++++------------- sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c | 13 ++++++---- sysdeps/ieee754/flt-32/s_logbf.c | 32 ++++++++++++++----------- sysdeps/ieee754/ldbl-128/s_logbl.c | 31 ++++++++++++++++-------- sysdeps/ieee754/ldbl-128ibm/s_logbl.c | 35 ++++++++++++++++----------- sysdeps/ieee754/ldbl-96/s_logbl.c | 35 ++++++++++++++++----------- 6 files changed, 111 insertions(+), 72 deletions(-) (limited to 'sysdeps/ieee754') diff --git a/sysdeps/ieee754/dbl-64/s_logb.c b/sysdeps/ieee754/dbl-64/s_logb.c index 2382fbb..baa35e1 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 #include -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 2ad6c7d..474eeef 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 , 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 diff --git a/sysdeps/ieee754/flt-32/s_logbf.c b/sysdeps/ieee754/flt-32/s_logbf.c index b6aa0f0..025c70d 100644 --- a/sysdeps/ieee754/flt-32/s_logbf.c +++ b/sysdeps/ieee754/flt-32/s_logbf.c @@ -13,23 +13,27 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_logbf.c,v 1.4 1995/05/10 20:47:51 jtc Exp $"; -#endif - #include #include -float __logbf(float x) +float +__logbf (float x) { - int32_t ix; - GET_FLOAT_WORD(ix,x); - ix &= 0x7fffffff; /* high |x| */ - if(ix==0) return (float)-1.0/fabsf(x); - if(ix>=0x7f800000) return x*x; - if((ix>>=23)==0) /* IEEE 754 logb */ - return -126.0; - else - return (float) (ix-127); + int32_t ix, rix; + + GET_FLOAT_WORD (ix, x); + ix &= 0x7fffffff; /* high |x| */ + if (ix == 0) + return (float) -1.0 / fabsf (x); + if (ix >= 0x7f800000) + return x * x; + if (__builtin_expect ((rix = ix >> 23) == 0, 0)) + { + /* POSIX specifies that denormal number is treated as + though it were normalized. */ + int m = (ix == 0) ? 0 : __builtin_clz (ix); + return -126.0 + (float)(8 - m); + } + return (float) (rix - 127); } weak_alias (__logbf, logbf) diff --git a/sysdeps/ieee754/ldbl-128/s_logbl.c b/sysdeps/ieee754/ldbl-128/s_logbl.c index 0b09b28..cf6003e 100644 --- a/sysdeps/ieee754/ldbl-128/s_logbl.c +++ b/sysdeps/ieee754/ldbl-128/s_logbl.c @@ -26,16 +26,27 @@ static char rcsid[] = "$NetBSD: $"; #include #include -long double __logbl(long double x) +long double +__logbl (long double x) { - int64_t lx,hx; - GET_LDOUBLE_WORDS64(hx,lx,x); - hx &= 0x7fffffffffffffffLL; /* high |x| */ - if((hx|lx)==0) return -1.0/fabs(x); - if(hx>=0x7fff000000000000LL) return x*x; - if((hx>>=48)==0) /* IEEE 754 logb */ - return -16382.0; - else - return (long double) (hx-0x3fff); + int64_t lx, hx, ex; + + GET_LDOUBLE_WORDS64 (hx, lx, x); + hx &= 0x7fffffffffffffffLL; /* high |x| */ + if ((hx | lx) == 0) + return -1.0 / fabs (x); + if (hx >= 0x7fff000000000000LL) + return x * x; + if ((ex = hx >> 48) == 0) /* IEEE 754 logb */ + { + /* POSIX specifies that denormal number is treated as + though it were normalized. */ + int m1 = (hx == 0) ? 0 : __builtin_clzll (hx); + int m2 = (lx == 0) ? 0 : __builtin_clzll (lx); + int ma = (m1 == 0) ? m2 + 64 : m1; + return -16382.0 + (long double)(15 - ma); + } + return (long double) (ex - 16383); } + weak_alias (__logbl, logbl) diff --git a/sysdeps/ieee754/ldbl-128ibm/s_logbl.c b/sysdeps/ieee754/ldbl-128ibm/s_logbl.c index f38b129..678b6ca 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_logbl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_logbl.c @@ -13,10 +13,6 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: $"; -#endif - /* * long double logbl(x) * IEEE 754 logb. Included to pass IEEE test suite. Not recommend. @@ -27,16 +23,27 @@ static char rcsid[] = "$NetBSD: $"; #include #include -long double __logbl(long double x) +long double +__logbl (long double x) { - int64_t lx,hx; - GET_LDOUBLE_WORDS64(hx,lx,x); - hx &= 0x7fffffffffffffffLL; /* high |x| */ - if((hx|(lx&0x7fffffffffffffffLL))==0) return -1.0/fabs(x); - if(hx>=0x7ff0000000000000LL) return x*x; - if((hx>>=52)==0) /* IEEE 754 logb */ - return -1022.0; - else - return (long double) (hx-0x3ff); + int64_t lx, hx, rhx; + + GET_LDOUBLE_WORDS64 (hx, lx, x); + hx &= 0x7fffffffffffffffLL; /* high |x| */ + if ((hx | (lx & 0x7fffffffffffffffLL)) == 0) + return -1.0 / fabs (x); + if (hx >= 0x7ff0000000000000LL) + return x * x; + if (__builtin_expect ((rhx = hx >> 52) == 0, 0)) + { + /* POSIX specifies that denormal number is treated as + though it were normalized. */ + int m1 = (hx == 0) ? 0 : __builtin_clzll (hx); + int m2 = (lx == 0) ? 0 : __builtin_clzll (lx); + int ma = (m1 == 0) ? m2 + 64 : m1; + return -1022.0 + (long double)(11 - ma); + } + return (long double) (rhx - 1023); } + long_double_symbol (libm, __logbl, logbl); diff --git a/sysdeps/ieee754/ldbl-96/s_logbl.c b/sysdeps/ieee754/ldbl-96/s_logbl.c index 95b644c..d8ad4bc 100644 --- a/sysdeps/ieee754/ldbl-96/s_logbl.c +++ b/sysdeps/ieee754/ldbl-96/s_logbl.c @@ -14,10 +14,6 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: $"; -#endif - /* * long double logbl(x) * IEEE 754 logb. Included to pass IEEE test suite. Not recommend. @@ -27,16 +23,27 @@ static char rcsid[] = "$NetBSD: $"; #include #include -long double __logbl(long double x) +long double +__logbl (long double x) { - int32_t es,lx,ix; - GET_LDOUBLE_WORDS(es,ix,lx,x); - es &= 0x7fff; /* exponent */ - if((es|ix|lx)==0) return -1.0/fabs(x); - if(es==0x7fff) return x*x; - if(es==0) /* IEEE 754 logb */ - return -16382.0; - else - return (long double) (es-0x3fff); + int32_t es, lx, ix; + + GET_LDOUBLE_WORDS (es, ix, lx, x); + es &= 0x7fff; /* exponent */ + if ((es | ix | lx) == 0) + return -1.0 / fabs (x); + if (es == 0x7fff) + return x * x; + if (es == 0) /* IEEE 754 logb */ + { + /* 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 -16382.0 - (long double)(ma); + } + return (long double) (es - 16383); } + weak_alias (__logbl, logbl) -- cgit v1.1