modf: replace MUSL version with Sun one for portability

diff --git a/Library/libs/modf.c b/Library/libs/modf.c
index 1128aed..083e967 100644 (file)
--- a/Library/libs/modf.c
+++ b/Library/libs/modf.c
@@ -1,40 +1,75 @@
-/* From MUSL */
-
 #include <math.h>
 #include "libm.h"
 
-double modf(double x, double *iptr)
-{
-       union {double x; uint64_t n;} u = {x};
-       uint64_t mask;
-       int e;
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
 
-       e = (int)(u.n>>52 & 0x7ff) - 0x3ff;
 
-       /* no fractional part */
-       if (e >= 52) {
-               *iptr = x;
-               if (e == 0x400 && u.n<<12 != 0) /* nan */
-                       return x;
-               u.n &= (uint64_t)1<<63;
-               return u.x;
-       }
+/*
+ * modf(double x, double *iptr)
+ * return fraction part of x, and return x's integral part in *iptr.
+ * Method:
+ *     Bit twiddling.
+ *
+ * Exception:
+ *     No exception.
+ */
 
-       /* no integral part*/
-       if (e < 0) {
-               u.n &= (uint64_t)1<<63;
-               *iptr = u.x;
-               return x;
-       }
+static const double one = 1.0;
 
-       mask = (uint64_t)-1>>12 >> e;
-       if ((u.n & mask) == 0) {
+double
+modf(double x, double *iptr)
+{
+       int32_t i0,i1,j0;
+       uint32_t i;
+       EXTRACT_WORDS(i0,i1,x);
+       j0 = ((i0>>20)&0x7ff)-0x3ff;    /* exponent of x */
+       if(j0<20) {                     /* integer part in high x */
+           if(j0<0) {                  /* |x|<1 */
+               INSERT_WORDS(*iptr,i0&0x80000000,0);    /* *iptr = +-0 */
+               return x;
+           } else {
+               i = (0x000fffff)>>j0;
+               if(((i0&i)|i1)==0) {            /* x is integral */
+                   uint32_t high;
+                   *iptr = x;
+                   GET_HIGH_WORD(high,x);
+                   INSERT_WORDS(x,high&0x80000000,0);  /* return +-0 */
+                   return x;
+               } else {
+                   INSERT_WORDS(*iptr,i0&(~i),0);
+                   return x - *iptr;
+               }
+           }
+       } else if (j0>51) {             /* no fraction part */
+           uint32_t high;
+           if (j0 == 0x400) {          /* inf/NaN */
                *iptr = x;
-               u.n &= (uint64_t)1<<63;
-               return u.x;
+               return 0.0 / x;
+           }
+           *iptr = x*one;
+           GET_HIGH_WORD(high,x);
+           INSERT_WORDS(x,high&0x80000000,0);  /* return +-0 */
+           return x;
+       } else {                        /* fraction part in low x */
+           i = ((uint32_t)(0xffffffff))>>(j0-20);
+           if((i1&i)==0) {             /* x is integral */
+               uint32_t high;
+               *iptr = x;
+               GET_HIGH_WORD(high,x);
+               INSERT_WORDS(x,high&0x80000000,0);      /* return +-0 */
+               return x;
+           } else {
+               INSERT_WORDS(*iptr,i0,i1&(~i));
+               return x - *iptr;
+           }
        }
-       u.n &= ~mask;
-       *iptr = u.x;
-       STRICT_ASSIGN(double, x, x - *iptr);
-       return x;
 }