home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / math / s_modf.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-16  |  1.8 KB  |  81 lines
  1.  
  2. /* @(#)s_modf.c 1.3 95/01/18 */
  3. /*
  4.  * ====================================================
  5.  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  6.  *
  7.  * Developed at SunSoft, 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 "fdlibm.h"
  25.  
  26. #ifdef __STDC__
  27. static const double one = 1.0;
  28. #else
  29. static double one = 1.0;
  30. #endif
  31.  
  32. #ifdef __STDC__
  33.     double modf(double x, double *iptr)
  34. #else
  35.     double modf(x, iptr)
  36.     double x,*iptr;
  37. #endif
  38. {
  39.     int i0,i1,j0;
  40.     unsigned i;
  41.     i0 =  __HI(x);        /* high x */
  42.     i1 =  __LO(x);        /* low  x */
  43.     j0 = ((i0>>20)&0x7ff)-0x3ff;    /* exponent of x */
  44.     if(j0<20) {            /* integer part in high x */
  45.         if(j0<0) {            /* |x|<1 */
  46.         __HIp(iptr) = i0&0x80000000;
  47.         __LOp(iptr) = 0;        /* *iptr = +-0 */
  48.         return x;
  49.         } else {
  50.         i = (0x000fffff)>>j0;
  51.         if(((i0&i)|i1)==0) {        /* x is integral */
  52.             *iptr = x;
  53.             __HI(x) &= 0x80000000;
  54.             __LO(x)  = 0;    /* return +-0 */
  55.             return x;
  56.         } else {
  57.             __HIp(iptr) = i0&(~i);
  58.             __LOp(iptr) = 0;
  59.             return x - *iptr;
  60.         }
  61.         }
  62.     } else if (j0>51) {        /* no fraction part */
  63.         *iptr = x*one;
  64.         __HI(x) &= 0x80000000;
  65.         __LO(x)  = 0;    /* return +-0 */
  66.         return x;
  67.     } else {            /* fraction part in low x */
  68.         i = ((unsigned)(0xffffffff))>>(j0-20);
  69.         if((i1&i)==0) {         /* x is integral */
  70.         *iptr = x;
  71.         __HI(x) &= 0x80000000;
  72.         __LO(x)  = 0;    /* return +-0 */
  73.         return x;
  74.         } else {
  75.         __HIp(iptr) = i0;
  76.         __LOp(iptr) = i1&(~i);
  77.         return x - *iptr;
  78.         }
  79.     }
  80. }
  81.