home *** CD-ROM | disk | FTP | other *** search
/ Garbo / Garbo.cdr / pc / source / byteunix.lzh / byte.2 / precision.c < prev    next >
C/C++ Source or Header  |  1990-05-11  |  6KB  |  195 lines

  1. /*******************************************************************************
  2.  *  The BYTE UNIX Benchmarks - Release 2
  3.  *          Module: precision.c   SID: 2.4 4/17/90 16:45:35
  4.  *          
  5.  *******************************************************************************
  6.  * Bug reports, patches, comments, suggestions should be sent to:
  7.  *
  8.  *    Ben Smith or Rick Grehan at BYTE Magazine
  9.  *    bensmith@bixpb.UUCP    rick_g@bixpb.UUCP
  10.  *
  11.  *******************************************************************************
  12.  *  Modification Log:
  13.  *  $Header: precision.c,v 3.5 87/08/06 08:10:59 kenj Exp $
  14.  *
  15.  ******************************************************************************/
  16. char SCCSid[] = "@(#) @(#)precision.c:2.4 -- 4/17/90 16:45:35";
  17. /* Program to determine properties of the arithmetic available. */
  18. /* Makes certain assumptions about the likely format of numbers */
  19. /*                                */
  20. /* Author: Steven Pemberton, CWI, Amsterdam. steven@mcvax    */
  21. /*                                */
  22. /* If your C system is not unix but does have signal/setjmp,    */
  23. /*    add a #define unix                    */
  24. /* You may also need to change the #include <sys/signal.h> line */
  25. /*    and add some calls to signal().                */
  26. /*
  27.  */
  28.  
  29. #ifdef unix
  30.  
  31. #define SIGNAL
  32.  
  33. #include <sys/signal.h>
  34. #include <setjmp.h>
  35.  
  36. jmp_buf lab;
  37. overflow(sig) int sig; { /*what to do on overflow*/
  38.     signal(sig, overflow);
  39.     longjmp(lab, 1);
  40. }
  41.  
  42. #endif
  43.  
  44. int tenlog(v) double v; {
  45.     /*The largest power of ten less than v*/
  46.     int p=0;
  47.     while (v>10) { p++; v/=10; }
  48.     return p;
  49. }
  50.  
  51. int two(v) int v; {
  52.     /* (the closest power of two to v)-1 */
  53.     int t=1, s;
  54.     while (t<v) t=t*2+1;
  55.     s=(t-1)/2;
  56.     if ((v-s)>(t-v)) return(t);
  57.     return(s);
  58. }
  59.  
  60. double twopower(n, e) int n, *e; {
  61.     /* Calculate 2**n without overflow worries */
  62.     /* Result is r*10**e */
  63.     double r=1.0; *e=0;
  64.     while (n-- > 0) {
  65.         r*=2.0;
  66.         if (r>10.0) { r/=10.0; (*e)++; }
  67.     }
  68.     return(r);
  69. }
  70.  
  71. main() {
  72.     short maxshort, newshort;
  73.     int maxint, newint, i, maxfexp, maxdexp, bits,
  74.         fmantis, dmantis, ddmantis,
  75.         shortpower, intpower, longpower,
  76.         fpower, dpower, fipower, dipower, ddipower, lpower, base;
  77.     long maxlong, newlong;
  78.     float maxfloat, newfloat, sum, f, maxifloat;
  79.     double maxdouble, newdouble, maxidouble, maxiexpr,
  80.            d, incr, dsum;
  81.  
  82. #ifdef SIGNAL
  83.     signal(SIGFPE, overflow); /*signal(SIGOVER, overflow);*/
  84. #endif
  85.  
  86. /****** Calculate max short *********************************************/
  87. /*      Calculate 2**n-1 until overflow - then use the previous value    */
  88.  
  89.     newshort=1; maxshort=0;
  90. #ifdef SIGNAL
  91.     if (setjmp(lab)==0)
  92. #endif
  93.     for(shortpower=0; newshort>maxshort; shortpower++) {
  94.         maxshort=newshort;
  95.         newshort=newshort*2+1;
  96.     }
  97.     bits= (shortpower+1)/sizeof(short);
  98.     printf("maxshort=%d (=2**%d-1)\n", maxshort, shortpower);
  99.  
  100. /****** Calculate max int by the same method ***************************/
  101.  
  102.     newint=1; maxint=0;
  103. #ifdef SIGNAL
  104.     if (setjmp(lab)==0)
  105. #endif
  106.     for(intpower=0; newint>maxint; intpower++) {
  107.         maxint=newint;
  108.         newint=newint*2+1;
  109.     }
  110.     printf("maxint=%d (=2**%d-1)\n", maxint, intpower);
  111.  
  112. /****** Calculate max long by the same method ***************************/
  113.  
  114.     newlong=1; maxlong=0;
  115. #ifdef SIGNAL
  116.     if (setjmp(lab)==0)
  117. #endif
  118.     for(longpower=0; newlong>maxlong; longpower++) {
  119.         maxlong=newlong;
  120.         newlong=newlong*2+1;
  121.     }
  122.     printf("maxlong=%ld (=2**%d-1)\n", maxlong, longpower);
  123.  
  124. /****** Calculate max float, assuming it's a power of two ***************/
  125. /*    Calculate 2**i until it overflows, and then use the nearest    */
  126. /*    power of two (some machines overflow early, some late)        */
  127.  
  128.     newfloat=1; maxfloat=0;
  129. #ifdef SIGNAL
  130.     if (setjmp(lab)==0)
  131. #endif
  132.     for(i=0;newfloat>maxfloat;i++) {
  133.         maxfloat=newfloat;
  134.         newfloat=newfloat*2;
  135.     }
  136.     maxfexp=two(i); maxfloat=twopower(maxfexp, &fpower);
  137.     printf("maxfloat=~%fE%d (=~2**%d) (%d bits)\n",
  138.             maxfloat, fpower, maxfexp, sizeof(float)*bits);
  139.  
  140. /****** Calculate max double, assuming it's a power of two **************/
  141.  
  142.     newdouble=1; maxdouble=0;
  143. #ifdef SIGNAL
  144.     if (setjmp(lab)==0)
  145. #endif
  146.     for(i=0;newdouble>maxdouble;i++) {
  147.         maxdouble=newdouble;
  148.         newdouble*=2;
  149.     }
  150. #ifdef SIGNAL
  151.     if (setjmp(lab)!=0) { printf("\nUnexpected overflow\n"); exit(1); }
  152. #endif
  153.     maxdexp=two(i); maxdouble=twopower(maxdexp, &dpower);
  154.     printf("maxdouble=~%fE%d (=~2**%d) (%d bits)\n",
  155.             maxdouble, dpower, maxdexp, sizeof(double)*bits);
  156.  
  157. /****** Calculate the accuracy for float, double, and expressions *******/
  158. /*    maxintfloat and maxintdouble are the largest values that can    */
  159. /*    still be integer values; ie such that (x+1)-x=1.        */
  160. /*    Some systems really do use extra precision in expressions    */
  161.  
  162.     f=2.0; incr=1.0; sum=f+incr;
  163.     for (fmantis=0; sum>2.0; fmantis++) { incr/=2; sum=f+incr; }
  164.     printf("max float exp=%d mantissa bits=%d\n", maxfexp, fmantis);
  165.  
  166.     d=2.0; incr=1.0; dsum=d+incr;
  167.     for (dmantis=0; dsum>2.0; dmantis++) { incr/=2; dsum=d+incr; }
  168.     printf("max double exp=%d mantissa bits=%d\n", maxdexp, dmantis);
  169.  
  170.     d=2.0; incr=1.0;
  171.     for (ddmantis=0; d+incr>2.0; ddmantis++) incr/=2;
  172.  
  173.     maxifloat=twopower(fmantis, &fipower);
  174.     printf("maxintfloat=~%fE%d (=2**%d-1) (%d digit precision)\n",
  175.             maxifloat, fipower, fmantis, fipower);
  176.  
  177.     maxidouble=twopower(dmantis, &dipower);
  178.     printf("maxintdouble=~%fE%d (=2**%d-1) (%d digit precision)\n",
  179.             maxidouble, dipower, dmantis, dipower);
  180.  
  181.     maxiexpr=twopower(ddmantis, &ddipower);
  182.     printf("maxint for expressions=~%fE%d (=2**%d-1) (%d digit precision)\n",
  183.             maxiexpr, ddipower, ddmantis, ddipower);
  184.  
  185. /****** BASE is the largest power of ten such that BASE*BASE can be    */
  186. /*    computed exactly as a double, and BASE+BASE as a long, useful    */
  187. /*    for multi-length arithmetic                    */
  188.  
  189.     lpower= tenlog((double)(maxlong/2));
  190.     base= (dipower/2)>lpower?lpower:(dipower/2);
  191.     printf("BASE=1E%d\n", base);
  192.     
  193.     exit(0);
  194. }
  195.