OpenE2K/glibc
2
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Fix for logb/logbf/logbl (bugs 13954/13955/13956)

Browse Source 
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.
This commit is contained in:
Adhemerval Zanella 2012-05-10 15:11:55 -05:00 committed by Ryan S. Arnold
parent 021db4be6f
commit 89c9aa491a
9 changed files with 143 additions and 74 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
ChangeLog
View File

@ -1,3 +1,16 @@
2012-05-09 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.
2012-05-09 Andreas Schwab <schwab@linux-m68k.org>
Andreas Jaeger <aj@suse.de>

5
NEWS
View File

@ -23,8 +23,9 @@ Version 2.16
13852, 13854, 13871, 13872, 13873, 13879, 13883, 13884, 13885, 13886,
13892, 13895, 13908, 13910, 13911, 13912, 13913, 13914, 13915, 13916,
13917, 13918, 13919, 13920, 13921, 13922, 13923, 13924, 13926, 13927,
13928, 13938, 13941, 13942, 13963, 13967, 13970, 13973, 13979, 13983,
14027, 14033, 14034, 14040, 14049, 14053, 14055, 14064, 14080, 14083
13928, 13938, 13941, 13942, 13954, 13955, 13956, 13963, 13967, 13970,
13973, 13979, 13983, 14027, 14033, 14034, 14040, 14049, 14053, 14055,
14064, 14080, 14083
* ISO C11 support:

16
math/libm-test.inc
View File

@ -5376,6 +5376,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);
}

37
sysdeps/ieee754/dbl-64/s_logb.c
View File

@ -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

13
sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c
View File

@ -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

32
sysdeps/ieee754/flt-32/s_logbf.c
View File

@ -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)

31
sysdeps/ieee754/ldbl-128/s_logbl.c
View File

@ -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_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)

35
sysdeps/ieee754/ldbl-128ibm/s_logbl.c
View File

@ -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_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);

35
sysdeps/ieee754/ldbl-96/s_logbl.c
View File

@ -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)