lib/libc/gen/modf.c

   1 /* @(#)s_modf.c 5.1 93/09/24 */
   2 /* $FreeBSD: head/lib/libc/gen/modf.c 226606 2011-10-21 06:40:36Z 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 <sys/types.h>
  25 #include <machine/endian.h>
  26 #include <math.h>
  27
  28 /* Bit fiddling routines copied from msun/src/math_private.h,v 1.15 */
  29
  30 #if BYTE_ORDER == BIG_ENDIAN
  31
  32 typedef union
  33 {
  34   double value;
  35   struct
  36   {
  37     u_int32_t msw;
  38     u_int32_t lsw;
  39   } parts;
  40 } ieee_double_shape_type;
  41
  42 #endif
  43
  44 #if BYTE_ORDER == LITTLE_ENDIAN
  45
  46 typedef union
  47 {
  48   double value;
  49   struct
  50   {
  51     u_int32_t lsw;
  52     u_int32_t msw;
  53   } parts;
  54 } ieee_double_shape_type;
  55
  56 #endif
  57
  58 /* Get two 32 bit ints from a double.  */
  59
  60 #define EXTRACT_WORDS(ix0,ix1,d)                                \
  61 do {                                                            \
  62   ieee_double_shape_type ew_u;                                  \
  63   ew_u.value = (d);                                             \
  64   (ix0) = ew_u.parts.msw;                                       \
  65   (ix1) = ew_u.parts.lsw;                                       \
  66 } while (0)
  67
  68 /* Get the more significant 32 bit int from a double.  */
  69
  70 #define GET_HIGH_WORD(i,d)                                      \
  71 do {                                                            \
  72   ieee_double_shape_type gh_u;                                  \
  73   gh_u.value = (d);                                             \
  74   (i) = gh_u.parts.msw;                                         \
  75 } while (0)
  76
  77 /* Set a double from two 32 bit ints.  */
  78
  79 #define INSERT_WORDS(d,ix0,ix1)                                 \
  80 do {                                                            \
  81   ieee_double_shape_type iw_u;                                  \
  82   iw_u.parts.msw = (ix0);                                       \
  83   iw_u.parts.lsw = (ix1);                                       \
  84   (d) = iw_u.value;                                             \
  85 } while (0)
  86
  87 static const double one = 1.0;
  88
  89 double
  90 modf(double x, double *iptr)
  91 {
  92         int32_t i0,i1,j0;
  93         u_int32_t i;
  94         EXTRACT_WORDS(i0,i1,x);
  95         j0 = ((i0>>20)&0x7ff)-0x3ff;    /* exponent of x */
  96         if(j0<20) {                     /* integer part in high x */
  97             if(j0<0) {                  /* |x|<1 */
  98                 INSERT_WORDS(*iptr,i0&0x80000000,0);    /* *iptr = +-0 */
  99                 return x;
 100             } else {
 101                 i = (0x000fffff)>>j0;
 102                 if(((i0&i)|i1)==0) {            /* x is integral */
 103                     u_int32_t high;
 104                     *iptr = x;
 105                     GET_HIGH_WORD(high,x);
 106                     INSERT_WORDS(x,high&0x80000000,0);  /* return +-0 */
 107                     return x;
 108                 } else {
 109                     INSERT_WORDS(*iptr,i0&(~i),0);
 110                     return x - *iptr;
 111                 }
 112             }
 113         } else if (j0>51) {             /* no fraction part */
 114             u_int32_t high;
 115             if (j0 == 0x400) {          /* inf/NaN */
 116                 *iptr = x;
 117                 return 0.0 / x;
 118             }
 119             *iptr = x*one;
 120             GET_HIGH_WORD(high,x);
 121             INSERT_WORDS(x,high&0x80000000,0);  /* return +-0 */
 122             return x;
 123         } else {                        /* fraction part in low x */
 124             i = ((u_int32_t)(0xffffffff))>>(j0-20);
 125             if((i1&i)==0) {             /* x is integral */
 126                 u_int32_t high;
 127                 *iptr = x;
 128                 GET_HIGH_WORD(high,x);
 129                 INSERT_WORDS(x,high&0x80000000,0);      /* return +-0 */
 130                 return x;
 131             } else {
 132                 INSERT_WORDS(*iptr,i0,i1&(~i));
 133                 return x - *iptr;
 134             }
 135         }
 136 }