lib/libm/src/s_modf.c

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