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