sysdeps/ieee754/dbl-64/s_modf.c

   1 /*
   2  * ====================================================
   3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   4  *
   5  * Developed at SunPro, a Sun Microsystems, Inc. business.
   6  * Permission to use, copy, modify, and distribute this
   7  * software is freely granted, provided that this notice
   8  * is preserved.
   9  * ====================================================
  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 <math.h>
  23 #include <math_private.h>
  24
  25 static const double one = 1.0;
  26
  27 double
  28 __modf (double x, double *iptr)
  29 {
  30   int32_t i0, i1, j0;
  31   u_int32_t i;
  32   EXTRACT_WORDS (i0, i1, x);
  33   j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;    /* exponent of x */
  34   if (j0 < 20)                          /* integer part in high x */
  35     {
  36       if (j0 < 0)                       /* |x|<1 */
  37         {
  38           INSERT_WORDS (*iptr, i0 & 0x80000000, 0);     /* *iptr = +-0 */
  39           return x;
  40         }
  41       else
  42         {
  43           i = (0x000fffff) >> j0;
  44           if (((i0 & i) | i1) == 0)             /* x is integral */
  45             {
  46               *iptr = x;
  47               INSERT_WORDS (x, i0 & 0x80000000, 0);     /* return +-0 */
  48               return x;
  49             }
  50           else
  51             {
  52               INSERT_WORDS (*iptr, i0 & (~i), 0);
  53               return x - *iptr;
  54             }
  55         }
  56     }
  57   else if (__glibc_unlikely (j0 > 51))              /* no fraction part */
  58     {
  59       *iptr = x * one;
  60       /* We must handle NaNs separately.  */
  61       if (j0 == 0x400 && ((i0 & 0xfffff) | i1))
  62         return x * one;
  63       INSERT_WORDS (x, i0 & 0x80000000, 0);     /* return +-0 */
  64       return x;
  65     }
  66   else                                  /* fraction part in low x */
  67     {
  68       i = ((u_int32_t) (0xffffffff)) >> (j0 - 20);
  69       if ((i1 & i) == 0)                /* x is integral */
  70         {
  71           *iptr = x;
  72           INSERT_WORDS (x, i0 & 0x80000000, 0);         /* return +-0 */
  73           return x;
  74         }
  75       else
  76         {
  77           INSERT_WORDS (*iptr, i0, i1 & (~i));
  78           return x - *iptr;
  79         }
  80     }
  81 }
  82 weak_alias (__modf, modf)
  83 #ifdef NO_LONG_DOUBLE
  84 strong_alias (__modf, __modfl)
  85 weak_alias (__modf, modfl)
  86 #endif