Merge commit 'origin/master' into prelink

Conflicts:

	ldso/include/dl-hash.h

Signed-off-by: Carmelo Amoroso <carmelo.amoroso@st.com>

8 files changed, 177 insertions, 267 deletions

diff --git a/libm/e_lgamma_r.c b/libm/e_lgamma_r.c
index 6d5d9c514..dbb0a0d5d 100644
--- a/libm/e_lgamma_r.c
+++ b/libm/e_lgamma_r.c
@@ -360,31 +360,35 @@ strong_alias(__ieee754_lgamma, gamma)
 #endif
 
 
-/* FIXME! Looks like someone just used __ieee754_gamma_r,
- * believing it's a "true" gamma function, but it was not!
- * Our tgamma is WRONG.
- */
-
 /* double tgamma(double x)
  * Return the Gamma function of x.
  */
 double tgamma(double x)
 {
-	double y;
-	int local_signgam;
+	int sign_of_gamma;
+	int32_t hx;
+	u_int32_t lx;
 
-	y = __ieee754_lgamma_r(x, &local_signgam); // was __ieee754_gamma_r
-	if (local_signgam < 0)
-		y = -y;
-#ifndef _IEEE_LIBM
-	if (_LIB_VERSION == _IEEE_)
-		return y;
-	if (!isfinite(y) && isfinite(x)) {
-		if (floor(x) == x && x <= 0.0)
-			return __kernel_standard(x, x, 41); /* tgamma pole */
-		return __kernel_standard(x, x, 40); /* tgamma overflow */
+	/* We don't have a real gamma implementation now.  We'll use lgamma
+	   and the exp function.  But due to the required boundary
+	   conditions we must check some values separately.  */
+
+	EXTRACT_WORDS(hx, lx, x);
+
+	if (((hx & 0x7fffffff) | lx) == 0) {
+		/* Return value for x == 0 is Inf with divide by zero exception.  */
+		return 1.0 / x;
 	}
-#endif
-	return y;
+	if (hx < 0 && (u_int32_t)hx < 0xfff00000 && rint(x) == x) {
+		/* Return value for integer x < 0 is NaN with invalid exception.  */
+		return (x - x) / (x - x);
+	}
+	if ((u_int32_t)hx == 0xfff00000 && lx == 0) {
+		/* x == -Inf.  According to ISO this is NaN.  */
+		return x - x;
+	}
+
+	x = exp(lgamma_r(x, &sign_of_gamma));
+	return sign_of_gamma >= 0 ? x : -x;
 }
 libm_hidden_def(tgamma)
diff --git a/libm/e_pow.c b/libm/e_pow.c
index 137f600c3..3be900389 100644
--- a/libm/e_pow.c
+++ b/libm/e_pow.c
@@ -21,25 +21,26 @@
  *	3. Return x**y = 2**n*exp(y'*log2)
  *
  * Special cases:
- *	1.  (anything) ** 0  is 1
- *	2.  (anything) ** 1  is itself
- *	3.  (anything) ** NAN is NAN
- *	4.  NAN ** (anything except 0) is NAN
- *	5.  +-(|x| > 1) **  +INF is +INF
- *	6.  +-(|x| > 1) **  -INF is +0
- *	7.  +-(|x| < 1) **  +INF is +0
- *	8.  +-(|x| < 1) **  -INF is +INF
- *	9.  +-1         ** +-INF is NAN
- *	10. +0 ** (+anything except 0, NAN)               is +0
- *	11. -0 ** (+anything except 0, NAN, odd integer)  is +0
- *	12. +0 ** (-anything except 0, NAN)               is +INF
- *	13. -0 ** (-anything except 0, NAN, odd integer)  is +INF
- *	14. -0 ** (odd integer) = -( +0 ** (odd integer) )
- *	15. +INF ** (+anything except 0,NAN) is +INF
- *	16. +INF ** (-anything except 0,NAN) is +0
- *	17. -INF ** (anything)  = -0 ** (-anything)
- *	18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
- *	19. (-anything except 0 and inf) ** (non-integer) is NAN
+ *	1.  +-1 ** anything  is 1.0
+ *	2.  +-1 ** +-INF     is 1.0
+ *	3.  (anything) ** 0  is 1
+ *	4.  (anything) ** 1  is itself
+ *	5.  (anything) ** NAN is NAN
+ *	6.  NAN ** (anything except 0) is NAN
+ *	7.  +-(|x| > 1) **  +INF is +INF
+ *	8.  +-(|x| > 1) **  -INF is +0
+ *	9.  +-(|x| < 1) **  +INF is +0
+ *	10  +-(|x| < 1) **  -INF is +INF
+ *	11. +0 ** (+anything except 0, NAN)               is +0
+ *	12. -0 ** (+anything except 0, NAN, odd integer)  is +0
+ *	13. +0 ** (-anything except 0, NAN)               is +INF
+ *	14. -0 ** (-anything except 0, NAN, odd integer)  is +INF
+ *	15. -0 ** (odd integer) = -( +0 ** (odd integer) )
+ *	16. +INF ** (+anything except 0,NAN) is +INF
+ *	17. +INF ** (-anything except 0,NAN) is +0
+ *	18. -INF ** (anything)  = -0 ** (-anything)
+ *	19. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
+ *	20. (-anything except 0 and inf) ** (non-integer) is NAN
  *
  * Accuracy:
  *	pow(x,y) returns x**y nearly rounded. In particular
@@ -99,8 +100,14 @@ double attribute_hidden __ieee754_pow(double x, double y)
 	u_int32_t lx,ly;
 
 	EXTRACT_WORDS(hx,lx,x);
+    /* x==1: 1**y = 1 (even if y is NaN) */
+	if (hx==0x3ff00000 && lx==0) {
+		return x;
+	}
+	ix = hx&0x7fffffff;
+
 	EXTRACT_WORDS(hy,ly,y);
-	ix = hx&0x7fffffff;  iy = hy&0x7fffffff;
+	iy = hy&0x7fffffff;
 
     /* y==zero: x**0 = 1 */
 	if((iy|ly)==0) return one;
@@ -132,13 +139,13 @@ double attribute_hidden __ieee754_pow(double x, double y)
 
     /* special value of y */
 	if(ly==0) {
-	    if (iy==0x7ff00000) {	/* y is +-inf */
-	        if(((ix-0x3ff00000)|lx)==0)
-		    return  y - y;	/* inf**+-1 is NaN */
-	        else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
-		    return (hy>=0)? y: zero;
-	        else			/* (|x|<1)**-,+inf = inf,0 */
-		    return (hy<0)?-y: zero;
+	    if (iy==0x7ff00000) {       /* y is +-inf */
+	        if (((ix-0x3ff00000)|lx)==0)
+		    return one;	        /* +-1**+-inf is 1 (yes, weird rule) */
+	        if (ix >= 0x3ff00000)   /* (|x|>1)**+-inf = inf,0 */
+		    return (hy>=0) ? y : zero;
+	        /* (|x|<1)**-,+inf = inf,0 */
+		return (hy<0) ? -y : zero;
 	    }
 	    if(iy==0x3ff00000) {	/* y is  +-1 */
 		if(hy<0) return one/x; else return x;
@@ -146,7 +153,7 @@ double attribute_hidden __ieee754_pow(double x, double y)
 	    if(hy==0x40000000) return x*x; /* y is  2 */
 	    if(hy==0x3fe00000) {	/* y is  0.5 */
 		if(hx>=0)	/* x >= +0 */
-		return __ieee754_sqrt(x);
+		    return __ieee754_sqrt(x);
 	    }
 	}
 
diff --git a/libm/e_scalb.c b/libm/e_scalb.c
index 9244cf42f..3fd78a1c0 100644
--- a/libm/e_scalb.c
+++ b/libm/e_scalb.c
@@ -10,9 +10,9 @@
  */
 
 /*
- * __ieee754_scalb(x, fn) is provide for
- * passing various standard test suite. One
- * should use scalbn() instead.
+ * __ieee754_scalb(x, fn) is provided for
+ * passing various standard test suites.
+ * One should use scalbn() instead.
  */
 
 #include "math.h"
@@ -21,7 +21,6 @@
 
 double attribute_hidden __ieee754_scalb(double x, double fn)
 {
-	return scalbn(x,fn);
 	if (isnan(x)||isnan(fn)) return x*fn;
 	if (!isfinite(fn)) {
 	    if(fn>0.0) return x*fn;
@@ -35,9 +34,9 @@ double attribute_hidden __ieee754_scalb(double x, double fn)
 
 #if defined __UCLIBC_SUSV3_LEGACY__
 /*
- * wrapper scalb(double x, double fn) is provide for
- * passing various standard test suite. One
- * should use scalbn() instead.
+ * wrapper scalb(double x, double fn) is provided for
+ * passing various standard test suites.
+ * One should use scalbn() instead.
  */
 #ifndef _IEEE_LIBM
 double scalb(double x, double fn)
diff --git a/libm/float_wrappers.c b/libm/float_wrappers.c
index a13aac93a..82b7963e1 100644
--- a/libm/float_wrappers.c
+++ b/libm/float_wrappers.c
@@ -293,10 +293,9 @@ long_WRAPPER1(lround)
 float modff (float x, float *iptr)
 {
 	double y, result;
-	result = modf ( x, &y );
+	result = modf( x, &y );
 	*iptr = (float)y;
 	return (float) result;
-
 }
 #endif
 
diff --git a/libm/ldouble_wrappers.c b/libm/ldouble_wrappers.c
index bf7ae15f0..118a78f64 100644
--- a/libm/ldouble_wrappers.c
+++ b/libm/ldouble_wrappers.c
@@ -21,6 +21,11 @@ long double func##l(long double x) \
 { \
 	return (long double) func((double) x); \
 }
+#define WRAPPER2(func) \
+long double func##l(long double x, long double y) \
+{ \
+	return (long double) func((double) x, (double) y); \
+}
 #define int_WRAPPER1(func) \
 int func##l(long double x) \
 { \
@@ -37,96 +42,6 @@ long long func##l(long double x) \
 	return func((double) x); \
 }
 
-#if defined __i386__ && defined __OPTIMIZE__
-# undef WRAPPER1
-# undef int_WRAPPER1
-# undef long_WRAPPER1
-# undef long_long_WRAPPER1
-/* gcc 4.3.1 generates really ugly code with redundant pair of store/load:
- *	sub	$0x10,%esp
- *	fldt	0x14(%esp)
- *	fstpl	0x8(%esp)
- *	fldl	0x8(%esp) <-- ??
- *	fstpl	(%esp)    <-- ??
- *	call	function
- *	add	$0x10,%esp
- *	ret
- * I can hope newer gcc will eliminate that. However, I don't think
- * it will be smart enough to reuse argument stack space and use
- * jump instead of call. Let's do it by hand.
- * The asm below loads long double x into st(0), then stores it back
- * to the same location, but as a double. At this point, stack looks
- * exactly as "double func(double)" expects it to be.
- * The return value is returned in st(0) per ABI in both cases (returning
- * a long double or returning a double). So we can simply jump to func.
- * Using __GI_func in jump to make optimized intra-library jump.
- * gcc will still generate a useless "ret" after asm. Oh well...
- */
-# define WRAPPER1(func) \
-long double func##l(long double x) \
-{ \
-	long double st_top; \
-	__asm__ ( \
-	"	fldt	%1\n" \
-	"	fstpl	%1\n" \
-	"	jmp	" __stringify(__GI_##func) "\n" \
-	: "=t" (st_top) \
-	: "m" (x) \
-	); \
-	return st_top; \
-}
-# define int_WRAPPER1(func) \
-int func##l(long double x) \
-{ \
-	int ret; \
-	__asm__ ( \
-	"	fldt	%1\n" \
-	"	fstpl	%1\n" \
-	"	jmp	" __stringify(__GI_##func) "\n" \
-	: "=a" (ret) \
-	: "m" (x) \
-	); \
-	return ret; \
-}
-# define long_WRAPPER1(func) \
-long func##l(long double x) \
-{ \
-	long ret; \
-	__asm__ ( \
-	"	fldt	%1\n" \
-	"	fstpl	%1\n" \
-	"	jmp	" __stringify(__GI_##func) "\n" \
-	: "=a" (ret) \
-	: "m" (x) \
-	); \
-	return ret; \
-}
-# define long_long_WRAPPER1(func) \
-long long func##l(long double x) \
-{ \
-	long long ret; \
-	__asm__ ( \
-	"	fldt	%1\n" \
-	"	fstpl	%1\n" \
-	"	jmp	" __stringify(__GI_##func) "\n" \
-	: "=A" (ret) \
-	: "m" (x) \
-	); \
-	return ret; \
-}
-#endif /* __i386__ && __OPTIMIZE__ */
-
-#if defined __NO_LONG_DOUBLE_MATH
-# define int_WRAPPER_C99(func) /* not needed */
-# else
-# define int_WRAPPER_C99(func) \
-int func##l(long double x) \
-{ \
-    return func((double) x); \
-} \
-libm_hidden_def(func##l)
-#endif
-
 /* Implement the following, as defined by SuSv3 */
 #if 0
 long double acoshl(long double);
@@ -205,10 +120,7 @@ WRAPPER1(asin)
 #endif
 
 #ifdef L_atan2l
-long double atan2l (long double x, long double y)
-{
-	return (long double) atan2( (double)x, (double)y );
-}
+WRAPPER2(atan2)
 #endif
 
 #ifdef L_atanhl
@@ -235,10 +147,7 @@ WRAPPER1(ceil)
 #endif
 
 #ifdef L_copysignl
-long double copysignl (long double x, long double y)
-{
-	return (long double) copysign( (double)x, (double)y );
-}
+WRAPPER2(copysign)
 #endif
 
 #ifdef L_coshl
@@ -274,10 +183,7 @@ WRAPPER1(fabs)
 #endif
 
 #ifdef L_fdiml
-long double fdiml (long double x, long double y)
-{
-	return (long double) fdim( (double)x, (double)y );
-}
+WRAPPER2(fdim)
 #endif
 
 #ifdef L_floorl
@@ -292,24 +198,15 @@ long double fmal (long double x, long double y, long double z)
 #endif
 
 #ifdef L_fmaxl
-long double fmaxl (long double x, long double y)
-{
-	return (long double) fmax( (double)x, (double)y );
-}
+WRAPPER2(fmax)
 #endif
 
 #ifdef L_fminl
-long double fminl (long double x, long double y)
-{
-	return (long double) fmin( (double)x, (double)y );
-}
+WRAPPER2(fmin)
 #endif
 
 #ifdef L_fmodl
-long double fmodl (long double x, long double y)
-{
-	return (long double) fmod( (double)x, (double)y );
-}
+WRAPPER2(fmod)
 #endif
 
 #ifdef L_frexpl
@@ -320,19 +217,11 @@ long double frexpl (long double x, int *ex)
 #endif
 
 #ifdef L_gammal
-/* WRAPPER1(gamma) won't work, tries to call __GI_xxx,
- * and gamma() hasn't got one. */
-long double gammal(long double x)
-{
-	return (long double) gamma((double) x);
-}
+WRAPPER1(gamma)
 #endif
 
 #ifdef L_hypotl
-long double hypotl (long double x, long double y)
-{
-	return (long double) hypot( (double)x, (double)y );
-}
+WRAPPER2(hypot)
 #endif
 
 #ifdef L_ilogbl
@@ -367,11 +256,7 @@ WRAPPER1(log1p)
 #endif
 
 #ifdef L_log2l
-/* WRAPPER1(log2) won't work */
-long double log2l(long double x)
-{
-	return (long double) log2((double)x);
-}
+WRAPPER1(log2)
 #endif
 
 #ifdef L_logbl
@@ -405,33 +290,21 @@ WRAPPER1(nearbyint)
 #endif
 
 #ifdef L_nextafterl
-long double nextafterl (long double x, long double y)
-{
-	return (long double) nextafter( (double)x, (double)y );
-}
+WRAPPER2(nextafter)
 #endif
 
 /* Disabled in Makefile.in */
 #if 0 /* def L_nexttowardl */
-long double nexttowardl (long double x, long double y)
-{
-	return (long double) nexttoward( (double)x, (double)y );
-}
+WRAPPER2(nexttoward)
 libm_hidden_def(nexttowardl)
 #endif
 
 #ifdef L_powl
-long double powl (long double x, long double y)
-{
-	return (long double) pow( (double)x, (double)y );
-}
+WRAPPER2(pow)
 #endif
 
 #ifdef L_remainderl
-long double remainderl (long double x, long double y)
-{
-	return (long double) remainder( (double)x, (double)y );
-}
+WRAPPER2(remainder)
 #endif
 
 #ifdef L_remquol
@@ -494,34 +367,34 @@ WRAPPER1(trunc)
 #endif
 
 #ifdef L_significandl
-/* WRAPPER1(significand) won't work, tries to call __GI_xxx,
- * and significand() hasn't got one. */
-long double significandl(long double x)
-{
-	return (long double) significand((double) x);
-}
+WRAPPER1(significand)
 #endif
 
-#ifdef __DO_C99_MATH__
+#if defined __DO_C99_MATH__ && !defined __NO_LONG_DOUBLE_MATH
 
-#ifdef L___fpclassifyl
-int_WRAPPER_C99(__fpclassify)
-#endif
+# ifdef L___fpclassifyl
+int_WRAPPER1(__fpclassify)
+libm_hidden_def(__fpclassifyl)
+# endif
 
-#ifdef L___finitel
-int_WRAPPER_C99(__finite)
-#endif
+# ifdef L___finitel
+int_WRAPPER1(__finite)
+libm_hidden_def(__finitel)
+# endif
 
-#ifdef L___signbitl
-int_WRAPPER_C99(__signbit)
-#endif
+# ifdef L___signbitl
+int_WRAPPER1(__signbit)
+libm_hidden_def(__signbitl)
+# endif
 
-#ifdef L___isnanl
-int_WRAPPER_C99(__isnan)
-#endif
+# ifdef L___isnanl
+int_WRAPPER1(__isnan)
+libm_hidden_def(__isnanl)
+# endif
 
-#ifdef L___isinfl
-int_WRAPPER_C99(__isinf)
-#endif
+# ifdef L___isinfl
+int_WRAPPER1(__isinf)
+libm_hidden_def(__isinfl)
+# endif
 
 #endif
diff --git a/libm/s_ilogb.c b/libm/s_ilogb.c
index ce3029667..259ae7b55 100644
--- a/libm/s_ilogb.c
+++ b/libm/s_ilogb.c
@@ -10,9 +10,10 @@
  */
 
 /* ilogb(double x)
- * return the binary exponent of non-zero x
- * ilogb(0) = 0x80000001
- * ilogb(inf/NaN) = 0x7fffffff (no signal is raised)
+ * return the binary exponent of x
+ * ilogb(+-0) = FP_ILOGB0
+ * ilogb(+-inf) = INT_MAX
+ * ilogb(NaN) = FP_ILOGBNAN (no signal is raised)
  */
 
 #include "math.h"
@@ -22,21 +23,35 @@ int ilogb(double x)
 {
 	int32_t hx,lx,ix;
 
-	GET_HIGH_WORD(hx,x);
+	GET_HIGH_WORD(hx, x);
 	hx &= 0x7fffffff;
-	if(hx<0x00100000) {
-	    GET_LOW_WORD(lx,x);
-	    if((hx|lx)==0)
-		return 0x80000001;	/* ilogb(0) = 0x80000001 */
-	    else			/* subnormal x */
-		if(hx==0) {
-		    for (ix = -1043; lx>0; lx<<=1) ix -=1;
-		} else {
-		    for (ix = -1022,hx<<=11; hx>0; hx<<=1) ix -=1;
+
+	if (hx < 0x00100000) {
+		GET_LOW_WORD(lx, x);
+		if ((hx|lx)==0)  /* +-0, ilogb(0) = FP_ILOGB0 */
+			return FP_ILOGB0;
+		/* subnormal x */
+		ix = -1043;
+		if (hx != 0) {
+			ix = -1022;
+			lx = (hx << 11);
 		}
-	    return ix;
+		/* each leading zero mantissa bit makes exponent smaller */
+		for (; lx > 0; lx <<= 1)
+			ix--;
+		return ix;
 	}
-	else if (hx<0x7ff00000) return (hx>>20)-1023;
-	else return 0x7fffffff;
+
+	if (hx < 0x7ff00000) /* normal x */
+		return (hx>>20) - 1023;
+
+	if (FP_ILOGBNAN != (~0U >> 1)) {
+		GET_LOW_WORD(lx, x);
+		if (hx == 0x7ff00000 && lx == 0)  /* +-inf */
+			return ~0U >> 1; /* = INT_MAX */
+	}
+
+	/* NAN. ilogb(NAN) = FP_ILOGBNAN */
+	return FP_ILOGBNAN;
 }
 libm_hidden_def(ilogb)
diff --git a/libm/s_modf.c b/libm/s_modf.c
index 26c2f1cd5..f73d5fddb 100644
--- a/libm/s_modf.c
+++ b/libm/s_modf.c
@@ -37,10 +37,8 @@ double modf(double x, double *iptr)
 	    } else {
 		i = (0x000fffff)>>j0;
 		if(((i0&i)|i1)==0) {		/* x is integral */
-		    u_int32_t high;
 		    *iptr = x;
-		    GET_HIGH_WORD(high,x);
-		    INSERT_WORDS(x,high&0x80000000,0);	/* return +-0 */
+		    INSERT_WORDS(x,i0&0x80000000,0);	/* return +-0 */
 		    return x;
 		} else {
 		    INSERT_WORDS(*iptr,i0&(~i),0);
@@ -48,18 +46,17 @@ double modf(double x, double *iptr)
 		}
 	    }
 	} else if (j0>51) {		/* no fraction part */
-	    u_int32_t high;
 	    *iptr = x*one;
-	    GET_HIGH_WORD(high,x);
-	    INSERT_WORDS(x,high&0x80000000,0);	/* return +-0 */
+	    /* We must handle NaNs separately.  */
+	    if (j0 == 0x400 && ((i0 & 0xfffff) | i1))
+	      return x*one;
+	    INSERT_WORDS(x,i0&0x80000000,0);	/* return +-0 */
 	    return x;
 	} else {			/* fraction part in low x */
 	    i = ((u_int32_t)(0xffffffff))>>(j0-20);
 	    if((i1&i)==0) { 		/* x is integral */
-	        u_int32_t high;
 		*iptr = x;
-		GET_HIGH_WORD(high,x);
-		INSERT_WORDS(x,high&0x80000000,0);	/* return +-0 */
+		INSERT_WORDS(x,i0&0x80000000,0);	/* return +-0 */
 		return x;
 	    } else {
 	        INSERT_WORDS(*iptr,i0,i1&(~i));
diff --git a/libm/s_rint.c b/libm/s_rint.c
index 358ce76b4..06432c622 100644
--- a/libm/s_rint.c
+++ b/libm/s_rint.c
@@ -30,41 +30,57 @@ TWO52[2]={
 
 double rint(double x)
 {
-	int32_t i0,j0,sx;
+	int32_t i0, j0, sx;
 	u_int32_t i,i1;
-	double w,t;
+	double t;
+	/* We use w = x + 2^52; t = w - 2^52; trick to round x to integer.
+	 * This trick requires that compiler does not optimize it
+	 * by keeping intermediate result w in a register wider than double.
+	 * Declaring w volatile assures that value gets truncated to double
+	 * (unfortunately, it also forces store+load):
+	 */
+	volatile double w;
+
 	EXTRACT_WORDS(i0,i1,x);
-	sx = (i0>>31)&1;
-	j0 = ((i0>>20)&0x7ff)-0x3ff;
-	if(j0<20) {
-	    if(j0<0) {
-		if(((i0&0x7fffffff)|i1)==0) return x;
+	/* Unbiased exponent */
+	j0 = ((((u_int32_t)i0) >> 20)&0x7ff)-0x3ff;
+
+	if (j0 > 51) {
+		//Why bother? Just returning x works too
+		//if (j0 == 0x400)  /* inf or NaN */
+		//	return x+x;
+		return x;  /* x is integral */
+	}
+
+	/* Sign */
+	sx = ((u_int32_t)i0) >> 31;
+
+	if (j0<20) {
+	    if (j0<0) { /* |x| < 1 */
+		if (((i0&0x7fffffff)|i1)==0) return x;
 		i1 |= (i0&0x0fffff);
 		i0 &= 0xfffe0000;
 		i0 |= ((i1|-i1)>>12)&0x80000;
 		SET_HIGH_WORD(x,i0);
-	        w = TWO52[sx]+x;
-	        t =  w-TWO52[sx];
+		w = TWO52[sx]+x;
+		t = w-TWO52[sx];
 		GET_HIGH_WORD(i0,t);
 		SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
-	        return t;
+		return t;
 	    } else {
 		i = (0x000fffff)>>j0;
-		if(((i0&i)|i1)==0) return x; /* x is integral */
+		if (((i0&i)|i1)==0) return x; /* x is integral */
 		i>>=1;
-		if(((i0&i)|i1)!=0) {
-		    if(j0==19) i1 = 0x40000000; else
-		    i0 = (i0&(~i))|((0x20000)>>j0);
+		if (((i0&i)|i1)!=0) {
+		    if (j0==19) i1 = 0x40000000;
+		    else i0 = (i0&(~i))|((0x20000)>>j0);
 		}
 	    }
-	} else if (j0>51) {
-	    if(j0==0x400) return x+x;	/* inf or NaN */
-	    else return x;		/* x is integral */
 	} else {
 	    i = ((u_int32_t)(0xffffffff))>>(j0-20);
-	    if((i1&i)==0) return x;	/* x is integral */
+	    if ((i1&i)==0) return x;	/* x is integral */
 	    i>>=1;
-	    if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
+	    if ((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
 	}
 	INSERT_WORDS(x,i0,i1);
 	w = TWO52[sx]+x;