lib/libm/src/s_modf.c

   1 /* @(#)s_modf.c 5.1 93/09/24 */
   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 #include <sys/cdefs.h>
  14 #if defined(LIBM_SCCS) && !defined(lint)
  15 __RCSID("$NetBSD: s_modf.c,v 1.14 2010/01/27 14:07:41 drochner Exp $");
  16 #endif
  17
  18 /*
  19  * modf(double x, double *iptr)
  20  * return fraction part of x, and return x's integral part in *iptr.
  21  * Method:
  22  *      Bit twiddling.
  23  *
  24  * Exception:
  25  *      No exception.
  26  */
  27
  28 #include "math.h"
  29 #include "math_private.h"
  30
  31 static const double one = 1.0;
  32
  33 double
  34 modf(double x, double *iptr)
  35 {
  36         int32_t i0,i1,jj0;
  37         u_int32_t i;
  38         EXTRACT_WORDS(i0,i1,x);
  39         jj0 = (((uint32_t)i0>>20)&0x7ff)-0x3ff; /* exponent of x */
  40         if(jj0<20) {                    /* integer part in high x */
  41             if(jj0<0) {                 /* |x|<1 */
  42                 INSERT_WORDS(*iptr,i0&0x80000000,0);    /* *iptr = +-0 */
  43                 return x;
  44             } else {
  45                 i = (0x000fffff)>>jj0;
  46                 if(((i0&i)|i1)==0) {            /* x is integral */
  47                     u_int32_t high;
  48                     *iptr = x;
  49                     GET_HIGH_WORD(high,x);
  50                     INSERT_WORDS(x,high&0x80000000,0);  /* return +-0 */
  51                     return x;
  52                 } else {
  53                     INSERT_WORDS(*iptr,i0&(~i),0);
  54                     return x - *iptr;
  55                 }
  56             }
  57         } else if (jj0>51) {            /* no fraction part */
  58             u_int32_t high;
  59             *iptr = x*one;
  60             if (jj0 == 0x400)           /* +-inf or NaN */
  61                 return 0.0 / x;         /* +-0 or NaN */
  62             GET_HIGH_WORD(high,x);
  63             INSERT_WORDS(x,high&0x80000000,0);  /* return +-0 */
  64             return x;
  65         } else {                        /* fraction part in low x */
  66             i = ((u_int32_t)(0xffffffff))>>(jj0-20);
  67             if((i1&i)==0) {             /* x is integral */
  68                 u_int32_t high;
  69                 *iptr = x;
  70                 GET_HIGH_WORD(high,x);
  71                 INSERT_WORDS(x,high&0x80000000,0);      /* return +-0 */
  72                 return x;
  73             } else {
  74                 INSERT_WORDS(*iptr,i0,i1&(~i));
  75                 return x - *iptr;
  76             }
  77         }
  78 }