* sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c (nextafterl): Remove

unused ily variable.  Fix nextafterl on +-__LDBL_MAX__ and +-Inf.
	Remove unreachable code at the end.

2007-06-01  Steven Munroe  <sjmunroe@us.ibm.com>

	* sysdeps/ieee754/ldbl-128ibm/s_fpclassifyl.c: Correct description of
	ldbl-128ibm in comment.
	(fpclassifyl): Correct classification of denormals.
	* sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c (nextafterl): Correct
	return value for MIN denormal. Rewrite using long double math too
	correctly handle denormals and canonicalize the results.

2007-06-05  Jakub Jelinek  <jakub@redhat.com>

	* sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c
	(__mpn_construct_long_double): Fix conversion where result ought
	to be smaller than __LDBL_MIN__, or the low double should be
	denormal.  Fix decision where to negate low double - honor round
	to even rules.
	* stdio-common/tst-sprintf2.c: Include string.h.
	(COMPARE_LDBL): Define.
	(TEST): Also test whether a string hexadecimal float representation
	can be parsed back to the number.
	(main): Add a couple of further tests.

2007-06-04  Jakub Jelinek  <jakub@redhat.com>
This commit is contained in:
Ulrich Drepper 2007-06-08 03:08:45 +00:00
parent 835abc5c0d
commit 7e3706ea25
5 changed files with 193 additions and 69 deletions

View file

@ -6,6 +6,34 @@
* sysdeps/x86_64/ldbl2mpn.c: New file.
* sysdeps/ia64/ldbl2mpn.c: New file.
2007-06-04 Jakub Jelinek <jakub@redhat.com>
* sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c (nextafterl): Remove
unused ily variable. Fix nextafterl on +-__LDBL_MAX__ and +-Inf.
Remove unreachable code at the end.
2007-06-01 Steven Munroe <sjmunroe@us.ibm.com>
* sysdeps/ieee754/ldbl-128ibm/s_fpclassifyl.c: Correct description of
ldbl-128ibm in comment.
(fpclassifyl): Correct classification of denormals.
* sysdeps/ieee754/ldbl-128ibm/s_nextafterl.c (nextafterl): Correct
return value for MIN denormal. Rewrite using long double math too
correctly handle denormals and canonicalize the results.
2007-06-05 Jakub Jelinek <jakub@redhat.com>
* sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c
(__mpn_construct_long_double): Fix conversion where result ought
to be smaller than __LDBL_MIN__, or the low double should be
denormal. Fix decision where to negate low double - honor round
to even rules.
* stdio-common/tst-sprintf2.c: Include string.h.
(COMPARE_LDBL): Define.
(TEST): Also test whether a string hexadecimal float representation
can be parsed back to the number.
(main): Add a couple of further tests.
2007-06-04 Jakub Jelinek <jakub@redhat.com>
* sysdeps/ieee754/ldbl-128ibm/printf_fphex.c

View file

@ -1,14 +1,22 @@
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
int
main (void)
{
volatile union { long double l; long long x[2]; } u;
volatile union { long double l; long long x[2]; } u, v;
char buf[64];
int result = 0;
#if LDBL_MANT_DIG == 106 || LDBL_MANT_DIG == 113
# define COMPARE_LDBL(u, v) \
((u).l == (v).l && (u).x[0] == (v).x[0] && (u).x[1] == (v).x[1])
#else
# define COMPARE_LDBL(u, v) ((u).l == (v).l)
#endif
#define TEST(val) \
do \
{ \
@ -19,6 +27,12 @@ main (void)
printf ("Error on line %d: %s != %s\n", __LINE__, buf, #val); \
result = 1; \
} \
if (sscanf (#val, "%La", &v.l) != 1 || !COMPARE_LDBL (u, v)) \
{ \
printf ("Error sscanf on line %d: %La != %La\n", __LINE__, \
u.l, v.l); \
result = 1; \
} \
/* printf ("%s %La %016Lx %016Lx\n", #val, u.l, u.x[0], u.x[1]); */ \
} \
while (0)
@ -54,6 +68,15 @@ main (void)
TEST (0x1.23456789abcdef123456789abc8p+64L);
TEST (0x1.23456789abcde7123456789abc8p+64L);
TEST (0x1.123456789abcdef123456789p-969L);
# if LDBL_MANT_DIG == 106
TEST (-0x1.2d71957cc1263bbbeb1d365f1e8p-969L);
TEST (0x1.23456789abcdef0123456789abp-970L);
TEST (0x1.579bde02468acp-1001L);
TEST (0x0.abcdef0123456p-1022L);
TEST (0x1.abcdef0123456p-1022L);
TEST (0x1.abcdef012345678p-1014L);
TEST (0x1.abcdef0123456f8p-1014L);
# endif
#endif
return result;
}

View file

@ -1,4 +1,4 @@
/* Copyright (C) 1995, 1996, 1997, 1998, 1999, 2002, 2003, 2004, 2006
/* Copyright (C) 1995, 1996, 1997, 1998, 1999, 2002, 2003, 2004, 2006, 2007
Free Software Foundation, Inc.
This file is part of the GNU C Library.
@ -31,19 +31,20 @@ long double
__mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
{
union ibm_extended_long_double u;
unsigned long hidden2, lzcount;
unsigned long lzcount;
unsigned long long hi, lo;
int exponent2;
u.ieee.negative = sign;
u.ieee.negative2 = sign;
u.ieee.exponent = expt + IBM_EXTENDED_LONG_DOUBLE_BIAS;
u.ieee.exponent2 = expt - 53 + IBM_EXTENDED_LONG_DOUBLE_BIAS;
u.ieee.exponent2 = 0;
exponent2 = expt - 53 + IBM_EXTENDED_LONG_DOUBLE_BIAS;
#if BITS_PER_MP_LIMB == 32
/* The low order 53 bits (52 + hidden) go into the lower double */
lo = frac_ptr[0];
lo |= (frac_ptr[1] & ((1LL << (53 - 32)) - 1)) << 32;
hidden2 = (frac_ptr[1] >> (52 - 32)) & ((mp_limb_t) 1);
/* The high order 53 bits (52 + hidden) go into the upper double */
hi = (frac_ptr[1] >> (53 - 32)) & ((1 << 11) - 1);
hi |= ((unsigned long long) frac_ptr[2]) << 11;
@ -51,7 +52,6 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
#elif BITS_PER_MP_LIMB == 64
/* The low order 53 bits (52 + hidden) go into the lower double */
lo = frac_ptr[0] & (((mp_limb_t) 1 << 53) - 1);
hidden2 = (frac_ptr[0] >> 52) & ((mp_limb_t) 1);
/* The high order 53 bits (52 + hidden) go into the upper double */
hi = (frac_ptr[0] >> 53) & (((mp_limb_t) 1 << 11) - 1);
hi |= (frac_ptr[1] << 11);
@ -59,14 +59,62 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
#error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
#endif
if ((hi & (1LL << 52)) == 0 && (hi | lo) != 0)
{
/* denormal number */
unsigned long long val = hi ? hi : lo;
if (sizeof (val) == sizeof (long))
lzcount = __builtin_clzl (val);
else if ((val >> 32) != 0)
lzcount = __builtin_clzl ((long) (val >> 32));
else
lzcount = __builtin_clzl ((long) val) + 32;
if (hi)
lzcount = lzcount - 11;
else
lzcount = lzcount + 42;
if (lzcount > u.ieee.exponent)
{
lzcount = u.ieee.exponent;
u.ieee.exponent = 0;
exponent2 -= lzcount;
}
else
{
u.ieee.exponent -= (lzcount - 1);
exponent2 -= (lzcount - 1);
}
if (lzcount <= 53)
{
hi = (hi << lzcount) | (lo >> (53 - lzcount));
lo = (lo << lzcount) & ((1LL << 53) - 1);
}
else
{
hi = lo << (lzcount - 53);
lo = 0;
}
}
if (lo != 0L)
{
/* hidden2 bit of low double controls rounding of the high double.
If hidden2 is '1' then round up hi and adjust lo (2nd mantissa)
If hidden2 is '1' and either the explicit mantissa is non-zero
or hi is odd, then round up hi and adjust lo (2nd mantissa)
plus change the sign of the low double to compensate. */
if (hidden2)
if ((lo & (1LL << 52)) != 0
&& ((hi & 1) != 0 || (lo & ((1LL << 52) - 1))))
{
hi++;
if ((hi & ((1LL << 52) - 1)) == 0)
{
if ((hi & (1LL << 53)) != 0)
hi -= 1LL << 52;
u.ieee.exponent++;
}
u.ieee.negative2 = !sign;
lo = (1LL << 53) - lo;
}
@ -85,17 +133,18 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
if (lzcount > 0)
{
lo = lo << lzcount;
u.ieee.exponent2 = u.ieee.exponent2 - lzcount;
exponent2 = exponent2 - lzcount;
}
if (exponent2 > 0)
u.ieee.exponent2 = exponent2;
else
lo >>= 1 - exponent2;
}
else
{
u.ieee.negative2 = 0;
u.ieee.exponent2 = 0;
}
u.ieee.negative2 = 0;
u.ieee.mantissa3 = lo & 0xffffffffLL;
u.ieee.mantissa2 = (lo >> 32) & 0xffffff;
u.ieee.mantissa2 = (lo >> 32) & 0xfffff;
u.ieee.mantissa1 = hi & 0xffffffffLL;
u.ieee.mantissa0 = (hi >> 32) & ((1LL << (LDBL_MANT_DIG - 86)) - 1);

View file

@ -1,5 +1,5 @@
/* Return classification value corresponding to argument.
Copyright (C) 1997,1999,2002,2004,2006 Free Software Foundation, Inc.
Copyright (C) 1997,1999,2002,2004,2006,2007 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
Jakub Jelinek <jj@ultra.linux.cz>, 1999.
@ -30,14 +30,16 @@
* -NaN fffn nnnn nnnn nnnn xxxx xxxx xxxx xxxx
* +Inf 7ff0 0000 0000 0000 xxxx xxxx xxxx xxxx
* -Inf fff0 0000 0000 0000 xxxx xxxx xxxx xxxx
* +0 0000 0000 0000 0000
* -0 8000 0000 0000 0000
* +normal 001n nnnn nnnn nnnn (smallest)
* -normal 801n nnnn nnnn nnnn (smallest)
* +normal 7fen nnnn nnnn nnnn (largest)
* -normal ffen nnnn nnnn nnnn (largest)
* +denorm 000n nnnn nnnn nnnn
* -denorm 800n nnnn nnnn nnnn
* +0 0000 0000 0000 0000 xxxx xxxx xxxx xxxx
* -0 8000 0000 0000 0000 xxxx xxxx xxxx xxxx
* +normal 0360 0000 0000 0000 0000 0000 0000 0000 (smallest)
* -normal 8360 0000 0000 0000 0000 0000 0000 0000 (smallest)
* +normal 7fef ffff ffff ffff 7c8f ffff ffff fffe (largest)
* +normal ffef ffff ffff ffff fc8f ffff ffff fffe (largest)
* +denorm 0360 0000 0000 0000 8000 0000 0000 0001 (largest)
* -denorm 8360 0000 0000 0000 0000 0000 0000 0001 (largest)
* +denorm 000n nnnn nnnn nnnn xxxx xxxx xxxx xxxx
* -denorm 800n nnnn nnnn nnnn xxxx xxxx xxxx xxxx
*/
int
@ -59,12 +61,23 @@ ___fpclassifyl (long double x)
/* +/-zero or +/- normal or +/- denormal */
if (hx & 0x7fffffffffffffffULL) {
/* +/- normal or +/- denormal */
if ((hx & 0x7ff0000000000000ULL) >= 0x0360000000000000ULL) {
if ((hx & 0x7ff0000000000000ULL) > 0x0360000000000000ULL) {
/* +/- normal */
retval = FP_NORMAL;
} else {
/* +/- denormal */
retval = FP_SUBNORMAL;
if ((hx & 0x7ff0000000000000ULL) == 0x0360000000000000ULL) {
if ((lx & 0x7fffffffffffffff) /* lower is non-zero */
&& ((lx^hx) & 0x8000000000000000ULL)) { /* and sign differs */
/* +/- denormal */
retval = FP_SUBNORMAL;
} else {
/* +/- normal */
retval = FP_NORMAL;
}
} else {
/* +/- denormal */
retval = FP_SUBNORMAL;
}
}
} else {
/* +/- zero */

View file

@ -35,7 +35,7 @@ static char rcsid[] = "$NetBSD: $";
long double x,y;
#endif
{
int64_t hx,hy,ihx,ihy,ilx,ily;
int64_t hx,hy,ihx,ihy,ilx;
u_int64_t lx,ly;
GET_LDOUBLE_WORDS64(hx,lx,x);
@ -43,7 +43,6 @@ static char rcsid[] = "$NetBSD: $";
ihx = hx&0x7fffffffffffffffLL; /* |hx| */
ilx = lx&0x7fffffffffffffffLL; /* |lx| */
ihy = hy&0x7fffffffffffffffLL; /* |hy| */
ily = ly&0x7fffffffffffffffLL; /* |ly| */
if((((ihx&0x7ff0000000000000LL)==0x7ff0000000000000LL)&&
((ihx&0x000fffffffffffffLL)!=0)) || /* x is nan */
@ -54,54 +53,66 @@ static char rcsid[] = "$NetBSD: $";
return y; /* x=y, return y */
if(ihx == 0 && ilx == 0) { /* x == 0 */
long double u;
SET_LDOUBLE_WORDS64(x,hy&0x8000000000000000ULL,1);/* return +-minsubnormal */
u = math_opt_barrier (u);
hy = (hy & 0x8000000000000000ULL) | 1;
SET_LDOUBLE_WORDS64(x,hy,0ULL);/* return +-minsubnormal */
u = math_opt_barrier (x);
u = u * u;
math_force_eval (u); /* raise underflow flag */
return x;
}
if(ihx>=0) { /* x > 0 */
if(ihx>ihy||((ihx==ihy)&&(ilx>ily))) { /* x > y, x -= ulp */
if(ilx==0)
hx--;
else
lx--;
} else { /* x < y, x += ulp */
if((hx==0x7fefffffffffffffLL)&&(lx==0x7c8ffffffffffffeLL))
{
SET_LDOUBLE_WORDS64(x,0x7ff0000000000000,0x8000000000000000);
return x;
}
else if((hx==0xffefffffffffffffLL)&&(lx==0xfc8ffffffffffffeLL))
{
SET_LDOUBLE_WORDS64(x,0xfff0000000000000,0x8000000000000000);
return x;
}
else if((lx&0x7fffffffffffffff)==0) hx++;
else
lx++;
long double u;
if(x > y) { /* x > y, x -= ulp */
if((hx==0xffefffffffffffffLL)&&(lx==0xfc8ffffffffffffeLL))
return x+x; /* overflow, return -inf */
if (hx >= 0x7ff0000000000000LL) {
SET_LDOUBLE_WORDS64(u,0x7fefffffffffffffLL,0x7c8ffffffffffffeLL);
return u;
}
} else { /* x < 0 */
if(ihy>=0||ihx>ihy||((ihx==ihy)&&(ilx>ily))){/* x < y, x -= ulp */
if((lx&0x7fffffffffffffff)==0)
hx--;
else
lx--;
} else { /* x > y, x += ulp */
if((lx&0x7fffffffffffffff)==0) hx++;
else
lx++;
if(ihx <= 0x0360000000000000LL) { /* x <= LDBL_MIN */
u = math_opt_barrier (x);
x -= __LDBL_DENORM_MIN__;
if (ihx < 0x0360000000000000LL
|| (hx > 0 && (int64_t) lx <= 0)
|| (hx < 0 && (int64_t) lx > 1)) {
u = u * u;
math_force_eval (u); /* raise underflow flag */
}
return x;
}
if (ihx < 0x06a0000000000000LL) { /* ulp will denormal */
SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL),0ULL);
u *= 0x1.0000000000000p-105L;
} else
SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL)-0x0690000000000000LL,0ULL);
return x - u;
} else { /* x < y, x += ulp */
if((hx==0x7fefffffffffffffLL)&&(lx==0x7c8ffffffffffffeLL))
return x+x; /* overflow, return +inf */
if ((u_int64_t) hx >= 0xfff0000000000000ULL) {
SET_LDOUBLE_WORDS64(u,0xffefffffffffffffLL,0xfc8ffffffffffffeLL);
return u;
}
if(ihx <= 0x0360000000000000LL) { /* x <= LDBL_MIN */
u = math_opt_barrier (x);
x += __LDBL_DENORM_MIN__;
if (ihx < 0x0360000000000000LL
|| (hx > 0 && (int64_t) lx < 0 && lx != 0x8000000000000001LL)
|| (hx < 0 && (int64_t) lx >= 0)) {
u = u * u;
math_force_eval (u); /* raise underflow flag */
}
if (x == 0.0L) /* handle negative __LDBL_DENORM_MIN__ case */
x = -0.0L;
return x;
}
if (ihx < 0x06a0000000000000LL) { /* ulp will denormal */
SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL),0ULL);
u *= 0x1.0000000000000p-105L;
} else
SET_LDOUBLE_WORDS64(u,(hx&0x7ff0000000000000LL)-0x0690000000000000LL,0ULL);
return x + u;
}
hy = hx&0x7ff0000000000000LL;
if(hy==0x7ff0000000000000LL) return x+x;/* overflow */
if(hy==0) {
long double u = x * x; /* underflow */
math_force_eval (u); /* raise underflow flag */
}
SET_LDOUBLE_WORDS64(x,hx,lx);
return x;
}
strong_alias (__nextafterl, __nexttowardl)
long_double_symbol (libm, __nextafterl, nextafterl);