nsz Git - libm/blob - src/math/modf.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/s_modf.c */
   2 /*
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunPro, a Sun Microsystems, Inc. business.
   7  * Permission to use, copy, modify, and distribute this
   8  * software is freely granted, provided that this notice
   9  * is preserved.
  10  * ====================================================
  11  */
  12 /*
  13  * modf(double x, double *iptr)
  14  * return fraction part of x, and return x's integral part in *iptr.
  15  * Method:
  16  *      Bit twiddling.
  17  *
  18  * Exception:
  19  *      No exception.
  20  */
  21
  22 #include "libm.h"
  23
  24 static const double one = 1.0;
  25
  26 double modf(double x, double *iptr)
  27 {
  28         int32_t i0,i1,j0;
  29         uint32_t i;
  30
  31         EXTRACT_WORDS(i0, i1, x);
  32         j0 = ((i0>>20) & 0x7ff) - 0x3ff; /* exponent of x */
  33         if (j0 < 20) {  /* integer part in high x */
  34                 if (j0 < 0) {  /* |x| < 1 */
  35                         INSERT_WORDS(*iptr, i0 & 0x80000000, 0); /* *iptr = +-0 */
  36                         return x;
  37                 }
  38                 i = 0x000fffff >> j0;
  39                 if (((i0&i)|i1) == 0) {  /* x is integral */
  40                         uint32_t high;
  41                         *iptr = x;
  42                         GET_HIGH_WORD(high, x);
  43                         INSERT_WORDS(x, high & 0x80000000, 0);  /* return +-0 */
  44                         return x;
  45                 }
  46                 INSERT_WORDS(*iptr, i0&~i, 0);
  47                 return x - *iptr;
  48         } else if (j0 > 51) {  /* no fraction part */
  49                 uint32_t high;
  50                 if (j0 == 0x400) {  /* inf/NaN */
  51                         *iptr = x;
  52                         return 0.0 / x;
  53                 }
  54                 *iptr = x*one;
  55                 GET_HIGH_WORD(high, x);
  56                 INSERT_WORDS(x, high & 0x80000000, 0);  /* return +-0 */
  57                 return x;
  58         } else {               /* fraction part in low x */
  59                 i = (uint32_t)0xffffffff >> (j0 - 20);
  60                 if ((i1&i) == 0) {  /* x is integral */
  61                         uint32_t high;
  62                         *iptr = x;
  63                         GET_HIGH_WORD(high, x);
  64                         INSERT_WORDS(x, high & 0x80000000, 0);  /* return +-0 */
  65                         return x;
  66                 }
  67                 INSERT_WORDS(*iptr, i0, i1&~i);
  68                 return x - *iptr;
  69         }
  70 }