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 #include <libm-alias-double.h>
  25
  26 static const double one = 1.0;
  27
  28 double
  29 __modf (double x, double *iptr)
  30 {
  31   int32_t i0, i1, j0;
  32   uint32_t i;
  33   EXTRACT_WORDS (i0, i1, x);
  34   j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;    /* exponent of x */
  35   if (j0 < 20)                          /* integer part in high x */
  36     {
  37       if (j0 < 0)                       /* |x|<1 */
  38         {
  39           INSERT_WORDS (*iptr, i0 & 0x80000000, 0);     /* *iptr = +-0 */
  40           return x;
  41         }
  42       else
  43         {
  44           i = (0x000fffff) >> j0;
  45           if (((i0 & i) | i1) == 0)             /* x is integral */
  46             {
  47               *iptr = x;
  48               INSERT_WORDS (x, i0 & 0x80000000, 0);     /* return +-0 */
  49               return x;
  50             }
  51           else
  52             {
  53               INSERT_WORDS (*iptr, i0 & (~i), 0);
  54               return x - *iptr;
  55             }
  56         }
  57     }
  58   else if (__glibc_unlikely (j0 > 51))              /* no fraction part */
  59     {
  60       *iptr = x * one;
  61       /* We must handle NaNs separately.  */
  62       if (j0 == 0x400 && ((i0 & 0xfffff) | i1))
  63         return x * one;
  64       INSERT_WORDS (x, i0 & 0x80000000, 0);     /* return +-0 */
  65       return x;
  66     }
  67   else                                  /* fraction part in low x */
  68     {
  69       i = ((uint32_t) (0xffffffff)) >> (j0 - 20);
  70       if ((i1 & i) == 0)                /* x is integral */
  71         {
  72           *iptr = x;
  73           INSERT_WORDS (x, i0 & 0x80000000, 0);         /* return +-0 */
  74           return x;
  75         }
  76       else
  77         {
  78           INSERT_WORDS (*iptr, i0, i1 & (~i));
  79           return x - *iptr;
  80         }
  81     }
  82 }
  83 #ifndef __modf
  84 libm_alias_double (__modf, modf)
  85 #endif