summaryrefslogtreecommitdiff
path: root/libm
diff options
context:
space:
mode:
authorCarmelo Amoroso <carmelo.amoroso@st.com>2010-11-28 21:09:47 +0100
committerCarmelo Amoroso <carmelo.amoroso@st.com>2010-11-28 21:09:47 +0100
commit57f6cce72c889f2e7b342e87a565a064320a6ee2 (patch)
tree5466406472918d3adb63398f5fb9ea1511a5e73f /libm
parent4ec89b87bc0eea8d9ca6b50564d12eeb3b1b0119 (diff)
parent5dea871ac73ef4608022c058b50adc946917f9b9 (diff)
Merge commit 'origin/master' into prelink
Conflicts: ldso/include/dl-hash.h Signed-off-by: Carmelo Amoroso <carmelo.amoroso@st.com>
Diffstat (limited to 'libm')
-rw-r--r--libm/e_lgamma_r.c42
-rw-r--r--libm/e_pow.c63
-rw-r--r--libm/e_scalb.c13
-rw-r--r--libm/float_wrappers.c3
-rw-r--r--libm/ldouble_wrappers.c207
-rw-r--r--libm/s_ilogb.c47
-rw-r--r--libm/s_modf.c15
-rw-r--r--libm/s_rint.c54
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;