home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_08_06 / 8n06122a < prev    next >
Text File  |  1990-05-02  |  14KB  |  419 lines
  1. /* Program to create "float.h" by running tests on target system
  2. ** using techniques from Kahan & Karpinski Paranoia code
  3. ** some of which are credited to W F Cody
  4. ** written by T C Prince Sept 1989
  5. ** tested on:
  6. **    SunOS 3.5 -ffpa
  7. **    SunOS 4.0 -fsoft, -f68881, -f68881 -fstore
  8. **    Multiflow 7/200 cc 1.6.2 -checkout, -O1, -O3
  9. **
  10. ** long double is simulated on non-ANSI compilers as register double
  11. */
  12. #include <math.h>
  13. #include <stdio.h>
  14. /* functions which narrow precision to float and double */
  15. #ifdef __STDC__
  16. float    storeff(double x);
  17. {
  18.     return x;
  19. }
  20.  
  21.  
  22. double    storef(long double x);
  23. {
  24.     return storeff(x);
  25. } /* narrow to float storage */
  26.  
  27.  
  28. double    store(long double x);
  29. {
  30.     return x;
  31. } /* force narrowing to double storage */
  32.  
  33.  
  34. #else
  35. long    storeff(x)
  36. long    x;
  37. {
  38.     return x;
  39. }
  40.  
  41.  
  42. double    storef(x)
  43. register double    x;
  44. {
  45.     union {
  46.         float    f;
  47.         long    l;
  48.     } xx;
  49.     xx.f = x;
  50.     xx.l = storeff(xx.l);
  51.     return (double)xx.f;
  52. }
  53.  
  54.  
  55. double    store(x)
  56. register double    x;
  57. {
  58.     return x;
  59. } /* force narrowing to double storage */
  60.  
  61.  
  62. #define SEEK_SET 0L
  63. #endif
  64. #define wp854(radix,precis) ((int)(log(radix)/log(10.)*precis+2))
  65. /* calculate width in decimal digits according to IEEE P854 */
  66. #define max(x,y) ((x)>(y)?(x):(y))
  67. #define min(x,y) ((x)<(y)?(x):(y))
  68. #define ODD(i) ((i)&1) /* like Pascal */
  69. #include <ctype.h> /* check what's there, it may not be too swift */
  70.  
  71. double    atof(dstr) char    *dstr; /*  any correctly working atof() may be used */
  72. {
  73.     int    c,
  74.     /* difference from exponent of accumulated integer to final value */
  75.     exp = 0,
  76.     neg_val = 0,/* leading negative sign ? */
  77.     e_exp = 0, /* number after e/E +/- */
  78.     neg_exp = 0;/* - after E/e ? */
  79. #ifdef __STDC__
  80.     long double x,retval = 0;
  81. #else
  82.     register double    x,retval = 0;
  83. #endif
  84.     while (isspace(c = *dstr)) /* skip leading  space */
  85.         dstr++;
  86.     switch (c) { /* optional sign */
  87.     case '-':
  88.         neg_val=1;
  89.     case '+': /* fall-through */
  90.         dstr++;
  91.     }
  92.  
  93.     while (isdigit(c = *dstr++)) {
  94.         /* add up converted values before dec point */
  95.         retval = 10 * retval + (double)(c - '0');
  96.     }
  97.     if ( c == '.')
  98.         while (isdigit(c = *dstr++)) {
  99.             --exp; /* and after dec pt */
  100.             retval = 10 * retval + (double)(c - '0');
  101.         };
  102.  
  103.     if ((c | ' ') == 'e' ) { /* found an exponent lead-in */
  104.  
  105.         switch (*dstr) { /*  sign field */
  106.         case '-':
  107.             neg_exp = 1;
  108.         case '+': /* fall-through */
  109.         case ' ': /* many FORTRAN environments generate this! */
  110.             dstr++;
  111.         default:
  112.             ;    /* fall-through */
  113.         }
  114.         if (isdigit(c = *dstr)) /* found  exponent digit */
  115.             do
  116.                 e_exp = 10 * e_exp + c - '0';
  117.             while (isdigit(c = *++dstr));
  118.     }
  119.     if (!neg_exp) {
  120.         /* we count on pow having a positive integer power case;
  121. ** a production shared library would use exp as an index into a table
  122. ** of powers of 10, but we don't know the size of the table yet;
  123. ** to reduce error and keep down the table size, we would have only
  124. ** positive powers, but partially underflowed numbers would require
  125. ** 2 steps */
  126.         if ((exp += e_exp) > 0)
  127.             retval *= pow(10., (double)exp);
  128.         else if (exp < 0 ) /* don't introduce extra error */
  129.             retval /= pow(10., (double)(-exp));
  130.     } else { /* this should avoid unnecessary underflow */
  131.         if(ODD(exp = e_exp-exp))retval /= 10;
  132.         for(x=10;exp>1;retval /= x)
  133.             do
  134.                 x *= x; /* generate 100,10000,... */
  135.             while(!ODD(exp >>= 1));
  136.     }
  137.     return (neg_val ? -retval : retval); /* apply sign */
  138. }
  139.  
  140.  
  141. #ifdef __STDC__
  142. long double    cvtest(x, wpi)
  143. int    wpi; /* ieee width  */
  144. long double    x;
  145. { /* find error in conversion to dec and back */
  146.     char    dec[50];
  147.     sprintf(dec, "%.*Lg", wpi, x);
  148.     return atof(dec);
  149. }
  150.  
  151.  
  152. #else
  153. double    cvtest(x, wpi)
  154. int    wpi; /* ieee width */
  155. double    x;
  156. {
  157.     char    dec[50];
  158.     sprintf(dec, "%.*g", wpi, x);
  159.     return atof(dec);
  160. }
  161.  
  162.  
  163. #endif
  164. /* begin main */
  165. main()
  166. {
  167.     FILE * stream; /* file handle for "float.h" output */
  168.     long    overfpos;/* position in file for DBL_MAX */
  169.     /* internal results: subtraction guarded? widths associated with precison */
  170.     int    subgrd, count, precis, wpd, wpf, wpl;
  171.     /* internal results corresponding to _RADIX, _EPSILON, _EPSILON/RADIX, 1/_EPSILON */
  172.     double    radix, ulpmd, ulppd, ulpmf, ulppf, ulppl, ulpml, w;
  173.     /* internal variables in largest accessible precision */
  174. #ifdef __STDC__
  175.     long double    y, z, v, sat;
  176. #else
  177.     register double    y, z, v, sat;
  178. #endif
  179.     /*
  180. **
  181. ** RADIX */
  182.     double    y1, y2, h, v1[16], v2[16]; /* internal test values */
  183.     stream = fopen("float.h", "w+");
  184.     /* find smallest power of 2 which is unaffected by adding 1 */
  185.     w = 1;
  186.     do
  187.         w += w;
  188.     while (store(w + 1) - w >= 1);
  189.     /* radix: smallest amount by which w may be increased */
  190.     y = 1;
  191.     while (!(radix = store(w + (y += y)) - w))/* stop as soon as non-zero */
  192.         ;
  193.     /* assume radix is same regardless of precision */
  194.     fprintf(stream, "#define FLT_RADIX\t%d\n", (int)radix);
  195.     /*
  196. **
  197. **  PRECISION */
  198.     w = 1;
  199.     precis = 0;
  200.     do {
  201.         ++precis;
  202.         w *= radix;
  203.     } while (store(w + 1) - w == 1); /* stop when error seen */
  204.     ulppd = (ulpmd = 1 / w) * radix;
  205.     wpd = wp854(radix, precis); /* corresponding decimal digits with guard */
  206.     fprintf(stream,
  207.         "#define DBL_EPSILON\t%.*g\n#define DBL_DIGITS\t%d\n#define DBL_MANT_DIG\t%d\n",
  208.         wpd, ulppd, (int)ceil(-log(ulppd) / log(10.)),precis);
  209.     w = 1;
  210.     precis = 0;
  211.     do {
  212.         ++precis;
  213.         w *= radix;
  214.     } while (storef(w + 1) - w == 1);
  215.     ulppf = (ulpmf = 1 / w) * radix;
  216.     wpf = wp854(radix, precis);
  217.     fprintf(stream,
  218.         "#define FLT_EPSILON\t%.*g\n#define FLT_DIGITS\t%d\n#define FLT_MANT_DIG\t%d\n",
  219.         wpf, ulppf, (int)ceil(-log(ulppf) / log(10.)),precis);
  220.     w = 1;
  221.     precis = 0;
  222.     do {
  223.         ++precis;
  224.         z = (w *= radix) + 1;
  225.     } while (z - w == 1);/* try to force order of operations */
  226.     ulppl = (ulpml = 1 / w) * radix;
  227.     wpl = wp854(radix, precis);
  228.     fprintf(stream,
  229.         "#define LDBL_EPSILON\t%.*g\n#define LDBL_DIGITS\t%d\n#define LDBL_MANT_DIG\t%d\n",
  230.         wpl, ulppl, (int)ceil(-log(ulppl) / log(10.)),precis);
  231.  
  232. /* ** MULTIPLICATION ROUNDING */
  233.     /* multiplication rounding test not ANSI-defined
  234. **      all rounding tests based on highest available precision
  235. **      and thus will vary with compiler options which change precision
  236. ** some known results:
  237. ** 68881 (Sun 4.0)
  238. **   multiply and divide are rounded in register precision, but compiler
  239. **   may round to double instead under certain circumstances, including 
  240. **   -fstore, they are rounded in double. add/subtract are rounded in
  241. **   register precision but suffer from    double rounding in double storage
  242. ** 8087
  243. **   all double operations suffer from double rounding
  244. **   register precision not accessible with most compilers
  245. ** Multiflow
  246. **   operations are rounded, but compiler converts v[]/y to v[]*1/y
  247. **   which suffers from double rounding
  248. **   -checkout (fc 1.6.2) fouls up addition rounding test
  249. ** Convex
  250. **   compiler reorders operations, often losing 1 bit of precision
  251. ** IBM 360 architecture
  252. **   all operations are chopped
  253. ** Cray-like architectures (including Convex)
  254. **   vector/scalar division suffers from double rounding
  255. **   or the pre-inversion itself may suffer from multiple rounding
  256. */
  257.     z = ulpml;
  258.     y = ulppl;
  259.     /* check for adequate guards to permit round up */
  260.     if (!(y1 = (1 + y) * (1 - y) - 1) & !(y2 = (1 + y + y) * (1 - y) - 1 - y))
  261.         fprintf(stream, "#define MUL_ROUNDS\t1\n");
  262.     else if (y1 < 0 & y2 == -y)
  263.         fprintf(stream, "#define MUL_ROUNDS\t0\n");
  264.     else
  265.         fprintf(stream, "#define MUL_ROUNDS\t-1\n");
  266. /*
  267. ** DIVISION ROUNDING */
  268.     /* test rounding of division, not ANSI C required */
  269.     for (count = -1; ++count < 16; ) { /* initialize vectors where loops cannot combine */
  270.         v1[count] = 1 - ulpmd - ulpmd;
  271.         v2[count] = 1 + ulppd + ulppd;
  272.     }/* now do register rounding tests */
  273.     if (!(y1 = (1 - z - z) / (1 - z) - 1 + z) & !(y2 = (1 + y + y) / (1 + y) - 1 - y) )
  274.         fprintf(stream, "#define DIV_ROUNDS\t1\n");
  275.     else if (y1 < 0 & y2 < 0)
  276.         fprintf(stream, "#define DIV_ROUNDS\t0\n");
  277.     else
  278.         fprintf(stream, "#define DIV_ROUNDS\t-1\n");
  279. /*
  280. ** VECTOR BY SCALAR DIVISION */
  281.     /* results may depend on compiler options
  282. ** same tests as above but memory to memory double arithmetic */
  283.     y1 = 1 - ulpmd;
  284.     y2 = 1 + ulppd;
  285.     for (count = -1; ++count < 16; ) {
  286.         v1[count] /= y1;
  287.         v2[count] /= y2;
  288.     }
  289.     if (!(y1 = v1[0] - y1) & !(y2 = v2[0] - y2) )
  290.         fprintf(stream, "#define VECTOR_BY_SCALAR_DIV_ROUNDS\t1\n");
  291.     else if (y1 < 0 & y2 < 0)
  292.         fprintf(stream, "#define VECTOR_BY_SCALAR_DIV_ROUNDS\t0\n");
  293.     else
  294.         fprintf(stream, "#define VECTOR_BY_SCALAR_DIV_ROUNDS\t-1\n");
  295.     /*
  296. **
  297. ** ADD/SUBTRACT ROUNDING -- ANSI REQUIRED */
  298.     /* find out if add/subtract are guarded; put in file ? */
  299.     /* K&R C allows cheating here but it probably doesn't matter
  300. ** unless the compiler fouls up like Convex which forms (1+z)-1
  301. ** resulting in z*2 thus calculating subgrd == 0 */
  302.     subgrd = 1 - (1 - z) == z & radix - (radix - y) == y;
  303.     /* test rounding of add/subtract with expressions close to threshold
  304. ** if rounding test fails, it may mean that the rounding threshold
  305. ** varies, as with fast software floating point, or it could be double
  306. ** rounding as on 8087 and 68881 */
  307.     if (1 + y * (1 - y) == 1 & 1 - z * z == 1 - z)
  308. /* much more detailed tests would be needed to distinguish FLT_ROUNDS == 3
  309. ** for obsolete machines which switch from 0 to 3 at -1 or -.5 */
  310.                 fprintf(stream, "#define FLT_ROUNDS\t0\n");
  311.     else if (1 + y == 1 + (y + .5) * y & 1 == 1 + (.5 - y) * y &
  312.         1 - z == 1 - (y + .5) * z & 1 == 1 - (.5 - y) * z & subgrd)
  313.         fprintf(stream, "#define FLT_ROUNDS\t1\n");
  314.     else
  315.         fprintf(stream, "#define FLT_ROUNDS\t-1\n");
  316. /*
  317. ** MINIMUM NORMALIZED VALUES */
  318.     /* test for partial underflow, turn off trapping */
  319.     z = (h = 1 / radix) * (1 + ulppd);
  320.     /* loop, dividing by radix, until irreversible result obtained
  321. ** there are ways to make this much faster */
  322.     do
  323.         y = z;
  324.     while (store(z *= h) * radix == y);
  325.     /* if conversion errors, adjust away from threshold
  326. ** using precautions to avoid abrupt underflow on non-IEEE machine */
  327.     y *= fabs(cvtest(y * radix, wpd) / (y * radix) - 1) + 1;
  328.     fprintf(stream,
  329.         "#define DBL_MIN\t\t%.*g\n#define DBL_MIN_10_EXP\t%d\n#define DBL_MIN_EXP\t%d\n",
  330.         wpd, y, (int)(log(y) / log(10.)),(int)(log(y)/log(radix)+.5));
  331.     z = h * (1 + ulppf);
  332.     do
  333.         y = z;
  334.     while (storef(z *= h) * radix == y);
  335.     y *= fabs(cvtest(y * radix, wpf) / (y * radix) - 1) + 1;
  336.     fprintf(stream,
  337.         "#define FLT_MIN\t\t%.*g\n#define FLT_MIN_10_EXP\t%d\n#define FLT_MIN_EXP\t%d\n",
  338.         wpf, y, (int)(log(y) / log(10.)),(int)(log(y)/log(radix)+.5));
  339.     z = h * (1 + (y = ulppl));
  340.     /* count used to obtain log base radix in case log() would underflow */
  341.     count = 0;
  342.     do {
  343.         ++count;
  344.         y = z;
  345.     } while ((z *= h) * radix == y);
  346.     y *= fabs(cvtest(y * radix, wpl) / (y * radix) - 1) + 1;
  347.     /* K&R libraries may fail to display meaningful values */
  348. #ifdef __STDC__
  349.     fprintf(stream,
  350.         "#define LDBL_MIN\t%.*Lg\n#define LDBL_MIN_10_EXP\t%d\n#define LDBL_MIN_EXP\t%d\n",
  351. #else
  352.         fprintf(stream,
  353.         "#define LDBL_MIN\t%.*g\n#define LDBL_MIN_10_EXP\t%d\n#define LDBL_MIN_EXP\t%d\n",
  354. #endif
  355.         wpl, y, (int)((1-count) / log(10.) * log(radix) ),1-count);
  356. /*
  357. ** MAXIMUM NORMALIZED VALUES */
  358.     /* keep overwriting overflow threshold in case of abort, but turn off
  359.     trapping if possible */
  360.     overfpos = ftell(stream);
  361.         z = radix * (y = -radix); /* - numbers may work best for 2's complement */
  362.     /* loop until multiplying by radix doesn't make it more negative */
  363.     do {
  364.         v = y;
  365.             if (fseek(stream, overfpos, SEEK_SET))
  366.             printf("seek failure in overflow search\n");
  367.             fprintf(stream,
  368.             "#define DBL_MAX\t\t%.17g\n#define DBL_MAX_10_EXP\t%d\n#define DBL_MAX_EXP\t%d",
  369.             -z, (int)(log(-z) / log(10.)),(int)(log(-z)/log(radix)+.5));
  370.             if (fflush(stream))printf("flush failure in overflow search\n");
  371.             z = store((y = z) * radix);         } while (z < y);
  372.         sat = (-y);
  373.         if ((z = (y = v * (radix * ulppd - radix)) + v * ((1 - radix) * ulppd)) < sat)y = z;
  374.         if (y < sat)v = y;
  375.         if (sat - fabs(v) < sat)v = sat;
  376.         /* if conversion errors, adjust away from threshold */
  377.         v -= fabs(cvtest(v * .75, wpd) - v * .75) * radix;
  378.         if (fseek(stream, overfpos, SEEK_SET))
  379.         printf("seek failure after overflow search\n");
  380.         fprintf(stream,
  381.         "#define DBL_MAX\t\t%.*g\n#define DBL_MAX_10_EXP\t%d\n#define DBL_MAX_EXP\t%d\n",
  382.         wpd, v, (int)(log(v) / log(10.) ),(int)(log(v)/log(radix)+.5));
  383.         fprintf(stream, "#ifndef HUGE_VAL\n#define HUGE_VAL\t%g\n#endif\n", sat);
  384.         z = radix * (y = -radix);
  385.         do {
  386.                v = y;
  387.             z = storef((y = z) * radix);         } while (z < y);
  388.         sat = (-y);
  389.         if ((z = (y = v * (radix * ulppf - radix)) + v * ((1 - radix) * ulppf)) < sat)y = z;
  390.         if (y < sat)v = y;
  391.         if (sat - fabs(v) < sat)v = sat;
  392.         v -= fabs(cvtest(v * .75, wpf) - v * .75) * radix;
  393.         fprintf(stream,
  394.         "#define FLT_MAX\t\t%.*g\n#define FLT_MAX_10_EXP\t%d\n#define FLT_MAX_EXP\t%d\n",
  395.         wpf, v, (int)(log(v) / log(10.) ),(int)(log(v)/log(radix)+.5));
  396.         z = radix * (y = -radix);
  397.         count = -1;
  398.         do {
  399.                     ++count;
  400.             v = y;
  401.             z = (y = z) * radix;         } while (z < y);
  402.         sat = (-y);
  403.         if ((z = (y = v * (radix * ulppl - radix)) + v * ((1 - radix) * ulppl)) < sat)y = z;
  404.         if (y < sat) v = y;
  405.         if (sat - fabs(v) < sat) {
  406.         ++count;
  407.                     v = sat;         }
  408.             v -= fabs(cvtest(v * .75, wpl) - v * .75) * radix;
  409.         /* K&R libraries may fail to display meaningful values */
  410. #ifdef __STDC__
  411.     fprintf(stream,
  412.         "#define LDBL_MAX\t%.*Lg\n#define LDBL_MAX_10_EXP\t%d\n#define LDBL_MAX_EXP\t%d\n",
  413. #else
  414.         fprintf(stream,
  415.         "#define LDBL_MAX\t%.*g\n#define LDBL_MAX_10_EXP\t%d\n#define LDBL_MAX_EXP\t%d\n",
  416. #endif
  417.         wpl, v, (int)((count+2) * log(radix) / log(10.) ),count+2);
  418.         exit(0); }
  419.