home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / altsrcs / 3 / 3137 / modf.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-03-27  |  3.3 KB  |  145 lines
  1. /*
  2.  * double modf(value, iptr)
  3.  * double value;
  4.  * double *iptr;
  5.  *
  6.  * returns fractional part of value
  7.  *       in *iptr returns the integral part
  8.  *       such that (*iptr + fractional) == value
  9.  *
  10.  * ++jrb    bammi@dsrgsun.ces.cwru.edu
  11.  *
  12.  * special thanks to peter housel  housel@ecn.purdue.edu
  13.  * this (improved) version is based on his X86 code
  14.  *
  15.  */
  16. #include "flonum.h"
  17.  
  18. #define BIAS 1023
  19. #define B1   1022    /* BIAS - 1 */
  20. #define B2   1075    /* (BIAS - 1) + 53 == max bias for integ part    */
  21. #define B3   1011    /* (BIAS - 1) - 11 == the max bias for a frac. # */
  22.  
  23. struct ldouble {
  24.     unsigned long hi, lo;
  25. };                /* another alias for double */
  26.  
  27. double modf(value, iptr)
  28. double value, *iptr;
  29. {
  30.     struct bitdouble *bd = (struct bitdouble *)&value;
  31.     unsigned int expo = bd->exp;
  32. #ifdef __STDC__
  33.     double norm( double, int, int, int );
  34. #else
  35.     extern double norm();
  36. #endif
  37.     
  38.     if(expo <= B1)    /* all frac, no integ */
  39.     {
  40.     *iptr = 0.0;
  41.     return value;
  42.     }
  43.     if(expo >= B2)    /* all integ, no frac */
  44.     {
  45.     *iptr = value;
  46.     return 0.0;
  47.     }
  48.  
  49.     /* both integ and frac */
  50.     {
  51.     register unsigned long i0, i1;    /* integral   part accumulator */
  52.     register unsigned long f0, f1;    /* fractional part accumulator */
  53.     unsigned int sign  = bd->sign;
  54.     
  55.     i0 = bd->mant1 | 0x00100000L;
  56.     i1 = bd->mant2;
  57.     f0 = f1 = 0L;
  58.  
  59.     do
  60.     {
  61.         /* shift R integ/frac, with bit falling into frac acc */
  62.         asm volatile(" lsrl  #1,%1;
  63.                 roxrl #1,%0"
  64.         : "=d" (i1), "=d" (i0)
  65.         : "0"  (i1), "1"  (i0) );
  66.  
  67.         asm volatile(" roxrl #1,%1;
  68.                 roxrl #1,%0"
  69.         : "=d" (f1), "=d" (f0)
  70.         : "0"  (f1), "1"  (f0) );
  71.  
  72.     } while(++expo < B2);
  73.  
  74.     /* stuff the results, and normalize */
  75.     ((struct ldouble *)&value)->hi = i0;    /* integral part */
  76.     ((struct ldouble *)&value)->lo = i1;
  77.     *iptr = norm(value, expo, sign, 0);
  78.  
  79.     /* dont call norm if frac is all zero, as B3 exp will screw it up */
  80.     if((f0 == 0) && (f1 == 0))
  81.         return 0.0;
  82.  
  83.     ((struct ldouble *)&value)->hi = f0;    /* fractional part */
  84.     ((struct ldouble *)&value)->lo = f1;
  85.     return norm(value, B3, sign, 0);
  86.     }
  87. }
  88.  
  89.  
  90. #ifdef TEST
  91. #include <stdio.h>
  92. #ifdef __MSHORT__
  93. #error "please run this test in 32 bit int mode"
  94. #endif
  95.  
  96. #define NTEST    100000L
  97. #define MAXRAND 0x7fffffffL    /* assuming 32 bit ints */
  98. #define ABS(X)    ( (X < 0) ? (-X) : X )
  99. extern long rand();
  100.  
  101. int main()
  102. {
  103.     double frac, integ, e, expected, result, max_abs_err;
  104.     register long i;
  105.     register long errs = 0;
  106.     register int s;
  107.  
  108.     max_abs_err = 0.0;
  109.     for(i = 0; i < NTEST; i++)
  110.     {
  111.     expected = ((double)(s = rand()) + (double)rand())/(double)(MAXRAND);
  112.     if(s > (MAXRAND >> 1)) expected = -expected;
  113.     frac = modf(expected, &integ);
  114.     if(ABS(frac) >= 1.0)
  115.     {
  116.         printf("|frac| >= 1, %.6e:  integ %.6e  frac %.6e\n",
  117.            expected, integ, frac);
  118.         errs++;
  119.         continue;
  120.     }
  121.     if( (integ != 0) && (ABS(integ) < 1.0) )
  122.     {
  123.         printf("|integ| < 1, %.6e:  integ %.6e  frac %.6e\n",
  124.            expected, integ, frac);
  125.         errs++;
  126.         continue;
  127.     }
  128.  
  129.     result = integ + frac;
  130.     e = (expected == 0.0) ? result : (result - expected)/expected;
  131.     if(e < 0) e = (-e);
  132.     if(e > 1.0e-6)
  133.     {
  134.         printf("%.4e: I %.4e F %.4e R %.4e E %.8e\n",
  135.            expected, integ, frac, result, e);
  136.         errs++;
  137.     }
  138.     if (e > max_abs_err) max_abs_err = e;
  139.     }
  140.     
  141.     printf("%ld Error(s), Max abs err %.14e\n", errs, max_abs_err);
  142.     return errs;
  143. }
  144. #endif /* TEST */
  145.