home *** CD-ROM | disk | FTP | other *** search
/ Serving the Web / ServingTheWeb1995.disc1of1.iso / linux / slacksrce / d / libc / libc-4.6 / libc-4 / libc-linux / sysdeps / i386 / modf.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-02-21  |  1.7 KB  |  78 lines
  1. #include <math.h>
  2.  
  3. #ifndef SOFT_387
  4.  
  5. double
  6. modf(double x, double *iptr)
  7. {
  8.   double tmp;
  9.   volatile short cw, cwtmp;
  10.  
  11.   __asm__ volatile ("fnstcw %0" : "=m" (cw) : );
  12.   cwtmp = cw | 0xc00;
  13.   __asm__ volatile ("fldcw %0" : : "m" (cwtmp));
  14.   __asm__ volatile ("frndint" : "=t" (tmp) : "0" (x));
  15.   __asm__ volatile ("fldcw %0" : : "m" (cw));
  16.   *iptr = tmp;
  17.  
  18.   return (x - tmp);
  19. }
  20.  
  21. #else
  22.  
  23. #include <ieee754.h>
  24.  
  25. #define shiftleft(dp,n)    {    /* n = 0 to 32 */ \
  26.     dp->mant1 = ((dp->mant1 << (n)) + (dp->mant2 >> (32-(n)))) \
  27.         & 0x0FFFFF; dp->mant2 <<= (n); dp->exp -= (n); }
  28.  
  29.  
  30. /*  Returns fractional part of d, stores integer part in *integ
  31.  */
  32. double modf(double d, double *integ)
  33. {
  34.     struct bitdouble *dp = (struct bitdouble *)&d;
  35.     struct bitdouble *ip = (struct bitdouble *)integ;
  36.     int e = dp->exp - BIAS;
  37.  
  38.     if (e < 0) {        /* no integer part */
  39.     *integ = 0;
  40.     return d;
  41.     }
  42.  
  43.     /* compute integer: clear fractional part where necessary 
  44.      */
  45.     *integ = d;
  46.     if (e <= 20) {
  47.     ip->mant1 &= (-1L << (20-e));        /* split in mant1... */
  48.     ip->mant2 = 0;
  49.     }
  50.     else 
  51.       if (e <= 52) 
  52.     ip->mant2 &= (-1L << (52-e));        /* split in mant2 */
  53.       else return 0;                /* no fractional part */
  54.  
  55.     /* compute fractional part: shift left over integer part
  56.      */
  57.     if (e)
  58.       if (e <= 32)
  59.     shiftleft(dp,e)
  60.       else {
  61.     dp->mant1 = (dp->mant2 << (e-32)) & 0x0FFFFF;
  62.     dp->mant2 = 0;
  63.     dp->exp -= e;
  64.       }
  65.  
  66.     /* adjust fractional part shifting left... 
  67.      */
  68.     if (dp->mant1==0 && dp->mant2==0)    /* fraction is zero */
  69.     return 0;
  70.  
  71.     while (!(dp->mant1 & 0x080000))     /* stack to the left */
  72.     shiftleft(dp,1);
  73.  
  74.     shiftleft(dp,1);            /* lose 'invisible bit' */
  75.     return d;
  76. }
  77. #endif
  78.