This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Re: [PATCH] Fix for logb/logbf/logbl (bz 13954/13955/13956)
- From: Andreas Jaeger <aj at suse dot com>
- To: libc-alpha at sourceware dot org, David Miller <davem at davemloft dot net>,Andreas Krebbel <Andreas dot Krebbel at de dot ibm dot com>
- Cc: Adhemerval Zanella <azanella at linux dot vnet dot ibm dot com>
- Date: Wed, 2 May 2012 14:15:51 +0200
- Subject: Re: [PATCH] Fix for logb/logbf/logbl (bz 13954/13955/13956)
- References: <4F9EEB79.20408@linux.vnet.ibm.com> <4F9FD94A.3080807@linux.vnet.ibm.com> <4FA04264.30708@linux.vnet.ibm.com>
On Tuesday, May 01, 2012 22:07:00 Adhemerval Zanella wrote:
> On 05/01/2012 09:38 AM, Adhemerval Zanella wrote:
> >> What about the dbl-64/wordsize-64, ldbl-96 and ldbl-128 routines? Do
> >> they have the same problem and need fixing as well?
> >> Your code change looks fine. Could you check the other
> >> implementations as well and resubmit then?
> >
> > I didn't check on above implementation but I'll do.
>
> I have added fixes for 'dbl-64/wordsize-64', 'ldbl-96' and 'ldbl-128',
> although for 'ldbl-128' I could not test (anyone to help me on this?).
> For 'dbl-64/wordsize-64' I have tested on PPC64 and for 'ldbl-96' on
> i386 and x86_64.
grepping Implies for ldbl-128 reveals:
So, let's ask Dave and Andreas for testing of this patch,
Andreas
> --
>
> 2012-05-01 Adhemerval Zanella <azanella@linux.vnet.ibm.com>
>
> [BZ #13954]
> [BZ #13955]
> [BZ #13956]
> * sysdeps/ieee754/dbl-64/s_logb.c (__logb): Fix for subnormal number.
> * sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c (__logb): Likewise.
> * sysdeps/ieee754/flt-32/s_logbf.c (__logf): Likewise.
> * sysdeps/ieee754/ldbl-128/s_logbl.c (__logbl): Likewise.
> * sysdeps/ieee754/ldbl-128ibm/s_logbl.c (__logbl): Likewise.
> * sysdeps/ieee754/ldbl-96/s_logbl.c (__logbl): Likewise.
> * math/libm-test.inc (logb_test) : Additional logb tests.
>
>
> diff --git a/math/libm-test.inc b/math/libm-test.inc
> index d643bad..8ce844b 100644
> --- a/math/libm-test.inc
> +++ b/math/libm-test.inc
> @@ -5364,6 +5364,22 @@ logb_test (void)
> TEST_f_f (logb, 1024, 10);
> TEST_f_f (logb, -2000, 10);
>
> + TEST_f_f (logb, 0x0.1p-127, -131);
> + TEST_f_f (logb, 0x0.01p-127, -135);
> + TEST_f_f (logb, 0x0.011p-127, -135);
> +#ifndef TEST_FLOAT
> + TEST_f_f (logb, 0x0.8p-1022, -1023);
> + TEST_f_f (logb, 0x0.1p-1022, -1026);
> + TEST_f_f (logb, 0x0.00111p-1022, -1034);
> + TEST_f_f (logb, 0x0.00001p-1022, -1042);
> + TEST_f_f (logb, 0x0.000011p-1022, -1042);
> + TEST_f_f (logb, 0x0.0000000000001p-1022, -1074);
> +#endif
> +#if defined TEST_LDOUBLE && LDBL_MIN_EXP - LDBL_MANT_DIG <= -16400
> + TEST_f_f (logb, 0x1p-16400L, -16400);
> + TEST_f_f (logb, 0x.00000000001p-16382L, -16426);
> +#endif
> +
> END (logb);
> }
>
> 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 <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 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 <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
> diff --git a/sysdeps/ieee754/flt-32/s_logbf.c
> b/sysdeps/ieee754/flt-32/s_logbf.c index b6aa0f0..166d1a4 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 <math.h>
> #include <math_private.h>
>
> -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..654b95e 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 <math.h>
> #include <math_private.h>
>
> -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_clzl (hx);
> + int m2 = (lx == 0) ? 0 : __builtin_clz; (lx);
> + int ma = (m1 == 0) ? m2 + 64 : m1;
> + return -16382.0 + (long double)(15 - ma);
> + }
> + return (long double) (hx - 0x3fff);
> }
> +
> weak_alias (__logbl, logbl)
> diff --git a/sysdeps/ieee754/ldbl-128ibm/s_logbl.c
> b/sysdeps/ieee754/ldbl-128ibm/s_logbl.c index f38b129..accd2c2 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 <math_private.h>
> #include <math_ldbl_opt.h>
>
> -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_clzl (hx);
> + int m2 = (lx == 0) ? 0 : __builtin_clzl (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..c556cdd 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 <math.h>
> #include <math_private.h>
>
> -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)
--
Andreas Jaeger aj@{suse.com,opensuse.org} Twitter/Identica: jaegerandi
SUSE LINUX Products GmbH, Maxfeldstr. 5, 90409 Nürnberg, Germany
GF: Jeff Hawn,Jennifer Guild,Felix Imendörffer,HRB16746 (AG Nürnberg)
GPG fingerprint = 93A3 365E CE47 B889 DF7F FED1 389A 563C C272 A126