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