diff options
| author | nsz <nsz@port70.net> | 2012-03-19 23:27:45 +0100 |
|---|---|---|
| committer | nsz <nsz@port70.net> | 2012-03-19 23:27:45 +0100 |
| commit | 75483499dad38b97f5dabb710635e6a8bbbb1c84 (patch) | |
| tree | 521da945aa4281d4c26ca22ead8ee2b30355ef0b /src/math | |
| parent | 2786c7d21611b9fa3b2fe356542cf213e7dd0ba4 (diff) | |
| download | musl-75483499dad38b97f5dabb710635e6a8bbbb1c84.tar.gz musl-75483499dad38b97f5dabb710635e6a8bbbb1c84.tar.xz musl-75483499dad38b97f5dabb710635e6a8bbbb1c84.zip | |
new modff.c code, fix nan handling in modfl
Diffstat (limited to 'src/math')
| -rw-r--r-- | src/math/modff.c | 66 | ||||
| -rw-r--r-- | src/math/modfl.c | 4 |
2 files changed, 26 insertions, 44 deletions
diff --git a/src/math/modff.c b/src/math/modff.c index d535314c..bf6e4ced 100644 --- a/src/math/modff.c +++ b/src/math/modff.c @@ -1,51 +1,33 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_modff.c */ -/* - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - #include "libm.h" -static const float one = 1.0; - float modff(float x, float *iptr) { - int32_t i0,j0; - uint32_t i; + uint32_t u, mask; + int e; - GET_FLOAT_WORD(i0, x); - j0 = ((i0>>23) & 0xff) - 0x7f; /* exponent of x */ - if (j0 < 23) { /* integer part in x */ - if (j0 < 0) { /* |x| < 1 */ - SET_FLOAT_WORD(*iptr, i0 & 0x80000000); /* *iptr = +-0 */ - return x; - } - i = 0x007fffff >> j0; - if ((i0&i) == 0) { /* x is integral */ - uint32_t ix; - *iptr = x; - GET_FLOAT_WORD(ix, x); - SET_FLOAT_WORD(x, ix & 0x80000000); /* return +-0 */ - return x; - } - SET_FLOAT_WORD(*iptr, i0&~i); - return x - *iptr; - } else { /* no fraction part */ - uint32_t ix; - *iptr = x*one; - if (x != x) /* NaN */ + GET_FLOAT_WORD(u, x); + e = (int)(u>>23 & 0xff) - 0x7f; + + /* no fractional part */ + if (e >= 23) { + *iptr = x; + if (e == 0x80 && u<<9 != 0) /* nan */ return x; - GET_FLOAT_WORD(ix, x); - SET_FLOAT_WORD(x, ix & 0x80000000); /* return +-0 */ + SET_FLOAT_WORD(x, u & 0x80000000); + return x; + } + /* no integral part */ + if (e < 0) { + SET_FLOAT_WORD(*iptr, u & 0x80000000); + return x; + } + + mask = 0x007fffff>>e; + if ((u & mask) == 0) { + *iptr = x; + SET_FLOAT_WORD(x, u & 0x80000000); return x; } + SET_FLOAT_WORD(*iptr, u & ~mask); + return x - *iptr; } diff --git a/src/math/modfl.c b/src/math/modfl.c index 2ca67b11..6520a1c2 100644 --- a/src/math/modfl.c +++ b/src/math/modfl.c @@ -54,7 +54,7 @@ long double modfl(long double x, long double *iptr) /* The number of fraction bits in manh, not counting the integer bit */ #define HIBITS (LDBL_MANT_DIG - LDBL_MANL_SIZE) -static const long double zero[] = { 0.0L, -0.0L }; +static const long double zero[] = { 0.0, -0.0 }; long double modfl(long double x, long double *iptr) { @@ -81,7 +81,7 @@ long double modfl(long double x, long double *iptr) return x - ux.e; } else if (e >= LDBL_MANT_DIG - 1) { /* x has no fraction part. */ *iptr = x; - if (x != x) /* Handle NaNs. */ + if (e == LDBL_MAX_EXP && (ux.bits.manh|ux.bits.manl)) /* nan */ return x; return zero[ux.bits.sign]; } else { /* Fraction part is in manl. */ |