nsz Git - musl/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 double modf(double x, double *iptr)
  25 {
  26         int32_t i0,i1,j0;
  27         uint32_t i;
  28
  29         EXTRACT_WORDS(i0, i1, x);
  30         j0 = ((i0>>20) & 0x7ff) - 0x3ff; /* exponent of x */
  31         if (j0 < 20) {  /* integer part in high x */
  32                 if (j0 < 0) {  /* |x| < 1 */
  33                         INSERT_WORDS(*iptr, i0 & 0x80000000, 0); /* *iptr = +-0 */
  34                         return x;
  35                 }
  36                 i = 0x000fffff >> j0;
  37                 if (((i0&i)|i1) == 0) {  /* x is integral */
  38                         uint32_t high;
  39                         *iptr = x;
  40                         GET_HIGH_WORD(high, x);
  41                         INSERT_WORDS(x, high & 0x80000000, 0);  /* return +-0 */
  42                         return x;
  43                 }
  44                 INSERT_WORDS(*iptr, i0&~i, 0);
  45                 return x - *iptr;
  46         } else if (j0 > 51) {  /* no fraction part */
  47                 uint32_t high;
  48                 if (j0 == 0x400) {  /* inf/NaN */
  49                         *iptr = x;
  50                         return 0.0 / x;
  51                 }
  52                 *iptr = x;
  53                 GET_HIGH_WORD(high, x);
  54                 INSERT_WORDS(x, high & 0x80000000, 0);  /* return +-0 */
  55                 return x;
  56         } else {               /* fraction part in low x */
  57                 i = (uint32_t)0xffffffff >> (j0 - 20);
  58                 if ((i1&i) == 0) {  /* x is integral */
  59                         uint32_t high;
  60                         *iptr = x;
  61                         GET_HIGH_WORD(high, x);
  62                         INSERT_WORDS(x, high & 0x80000000, 0);  /* return +-0 */
  63                         return x;
  64                 }
  65                 INSERT_WORDS(*iptr, i0, i1&~i);
  66                 return x - *iptr;
  67         }
  68 }