X-Git-Url: http://nsz.repo.hu/git/?a=blobdiff_plain;f=src%2Fmath%2Fmodfl.c;h=a47b1924f75d1125b712e3ccf5b0f39af27fd8f2;hb=c5f4b2dfea320356f69445dc1adf8f73596a3c36;hp=bbfcdb8a8c37c4094da847d1e5ef864add7affd8;hpb=fe2df4ad3ea6080d01dea8d836cd4bbc896160fd;p=musl

diff --git a/src/math/modfl.c b/src/math/modfl.c
index bbfcdb8a..a47b1924 100644
--- a/src/math/modfl.c
+++ b/src/math/modfl.c
@@ -1,100 +1,53 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/s_modfl.c */
-/*-
- * Copyright (c) 2007 David Schultz <das@FreeBSD.ORG>
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- *
- * Derived from s_modf.c, which has the following Copyright:
- * ====================================================
- * 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"
 
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 long double modfl(long double x, long double *iptr)
 {
-	return modf(x, (double *)iptr);
+	double d;
+	long double r;
+
+	r = modf(x, &d);
+	*iptr = d;
+	return r;
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
 
-#if LDBL_MANL_SIZE > 32
-#define MASK    ((uint64_t)-1)
-#else
-#define MASK    ((uint32_t)-1)
-#endif
-/* Return the last n bits of a word, representing the fractional part. */
-#define GETFRAC(bits, n)        ((bits) & ~(MASK << (n)))
-/* 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.0, -0.0 };
+static const long double toint = 1/LDBL_EPSILON;
 
 long double modfl(long double x, long double *iptr)
 {
-	union IEEEl2bits ux;
-	int e;
+	union ldshape u = {x};
+	int e = (u.i.se & 0x7fff) - 0x3fff;
+	int s = u.i.se >> 15;
+	long double absx;
+	long double y;
 
-	ux.e = x;
-	e = ux.bits.exp - LDBL_MAX_EXP + 1;
-	if (e < HIBITS) {                       /* Integer part is in manh. */
-		if (e < 0) {                    /* |x|<1 */
-			*iptr = zero[ux.bits.sign];
-			return x;
-		}
-		if ((GETFRAC(ux.bits.manh, HIBITS - 1 - e)|ux.bits.manl) == 0) {
-			/* x is an integer. */
-			*iptr = x;
-			return zero[ux.bits.sign];
-		}
-		/* Clear all but the top e+1 bits. */
-		ux.bits.manh >>= HIBITS - 1 - e;
-		ux.bits.manh <<= HIBITS - 1 - e;
-		ux.bits.manl = 0;
-		*iptr = ux.e;
-		return x - ux.e;
-	} else if (e >= LDBL_MANT_DIG - 1) {    /* x has no fraction part. */
+	/* no fractional part */
+	if (e >= LDBL_MANT_DIG-1) {
 		*iptr = x;
-		if (e == LDBL_MAX_EXP && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)) /* nan */
+		if (isnan(x))
 			return x;
-		return zero[ux.bits.sign];
-	} else {                                /* Fraction part is in manl. */
-		if (GETFRAC(ux.bits.manl, LDBL_MANT_DIG - 1 - e) == 0) {
-			/* x is integral. */
-			*iptr = x;
-			return zero[ux.bits.sign];
-		}
-		/* Clear all but the top e+1 bits. */
-		ux.bits.manl >>= LDBL_MANT_DIG - 1 - e;
-		ux.bits.manl <<= LDBL_MANT_DIG - 1 - e;
-		*iptr = ux.e;
-		return x - ux.e;
+		return s ? -0.0 : 0.0;
+	}
+
+	/* no integral part*/
+	if (e < 0) {
+		*iptr = s ? -0.0 : 0.0;
+		return x;
+	}
+
+	/* raises spurious inexact */
+	absx = s ? -x : x;
+	y = absx + toint - toint - absx;
+	if (y == 0) {
+		*iptr = x;
+		return s ? -0.0 : 0.0;
 	}
+	if (y > 0)
+		y -= 1;
+	if (s)
+		y = -y;
+	*iptr = x + y;
+	return -y;
 }
 #endif