home *** CD-ROM | disk | FTP | other *** search
/ AP Professional Graphics CD-ROM Library / AP Professional Graphics CD-ROM Library.iso / pc / unix / appendix / gemsiii / interval.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-12-14  |  5.0 KB  |  207 lines
  1. /* The following code implements the four basic interval arithmetic 
  2. operations +, -, * and / within the C++ class interval as provided 
  3. by the AT & T C++ Release 2.0 for the Sparc architecture. There is 
  4. no checking that the intervals are valid intervals, i.e. if the 
  5. interval is [lo.hi] then it is valid iff lo <= hi.
  6.  
  7. It is assumed that a mechanism exists to direct the result of a 
  8. floating point operation so that it is chopped. This corresponds
  9. to "tozero" in ANSI/IEEE Std. 754-1985 and it is performed by the 
  10. routine
  11.                ieee_flags
  12. in this code.
  13.  
  14. The first executable statement of any program using this code must 
  15. therefore be 
  16.                setmode(1);
  17. or a corresponding call in in the architecture used in order that 
  18. the rounding mode will be the chopped mode which is assumed in the 
  19. class. Changes have to be made to the routine help_ if another
  20. rounding mode is used.
  21. */
  22.  
  23.  
  24. #define    min(a, b)        (((a) < (b)) ? (a) : (b))
  25. #define    MIN(a, b, c, d)        min(min(min(a, b), c), d)
  26. #define    max(a, b)        (((a) > (b)) ? (a) : (b))
  27. #define    MAX(a, b, c, d)        max(max(max(a, b), c), d)
  28.  
  29. #include <stdio.h>
  30. #include <stream.h>
  31. #include <sys/ieeefp.h>
  32.  
  33. extern "C" {
  34.     int ieee_flags(char *, char *, char *in, char **);
  35.       }
  36. void setmode_(int *mode)
  37. {
  38. char *out;
  39.  
  40.     ieee_flags("set", "direction", *mode ? "tozero" : "nearest", &out);
  41.     ieee_flags("get", "direction", NULL, &out);
  42. }
  43.  
  44.  
  45.  
  46. /*
  47.  
  48. This routine adds 1 to the last digit of a double data item which is here
  49. assumed to have 8 bytes where
  50.  
  51. seeeeeee eeeemmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmmm mmmmmmmm
  52.  
  53. s=sign of number
  54. e=exponent            01111111111 is an exponent of 0
  55. m=mantissa bits with assumed leading 1
  56.  
  57. The algorithm works by adding one to the last two bytes of the mantissa.
  58. If the result of this addition is zero then it overflows to the next two
  59. bytes. This means that the next two bytes must be checked for overflow 
  60. until the mantissa bits are exhausted. It should then be noted that the
  61. exponent will be automatically incremented and the mantissa will have 
  62. been set to all zeros.
  63.  
  64. No check is made to see if the exponent overflows. Negative numbers will 
  65. turn into the next larger negative number.
  66.  
  67. This routine assumes that the number of bytes in a double is and integer
  68. multiple of the number of bytes in an unsigned integer.
  69. */
  70.  
  71.  
  72. void helpx_(double *a)
  73. {
  74. unsigned char *b;
  75. unsigned int *c,i;
  76.  
  77.     b=(unsigned char *)a;
  78.     i=sizeof(double)-sizeof(unsigned int);
  79.     do {
  80.         c=(unsigned int *)(b+i);
  81.         (*c)++;
  82.         i-=sizeof(unsigned int);
  83.     } while ((*c)==0 && i>=0);
  84. }
  85.  
  86. /* A printing utility for printing the octal contents of a double 
  87. data item */
  88.  
  89. void print_bin(double a)
  90. {
  91. char *b;
  92. int i,j,k;
  93.  
  94.     k=0;
  95.     b=(char *)(&a);
  96.     printf("%f ",a);
  97.     for (i=0;i<sizeof(double);i++) {
  98.         for(j=0;j<8;j++) {
  99.             k++;
  100.             if (k==2 || k==13) putchar(' ');
  101.             printf("%d",(b[i]&0x80)==128);
  102.             b[i]<<=1;
  103.         }
  104.     }
  105.     printf("\n");
  106. }
  107.  
  108. /*
  109. The basic arithmetic operations are embedded in the class interval. In
  110. each case the pair of double data items representing the interval [lo,hi]
  111. is checked via the routine help_. This routine accesses helpx_ whenever
  112.           lo<0
  113. or
  114.           hi>0.
  115. In both cases the next double representable item is calculated in helpx_
  116. and the result is returned to the calling operator. In this manner it is
  117. guaranteed that the machine interval result will contain the interval
  118. result that would have been computed using infinite (real) interval
  119. arithmetic.
  120. */
  121.  
  122. class    interval {
  123.         double lo,hi;
  124.  
  125.  
  126. public:
  127.     interval(double lo = 0, double hi = 0) {
  128.         this->lo = lo;
  129.         this->hi = hi;
  130.     };
  131.  
  132. /* An interval printing utility */
  133.         void print() {
  134.         printf("[ %3lx, %3lx ]\n", lo, hi);
  135.     };
  136.  
  137.     friend interval help_(interval a) {
  138.       if(a.lo < 0 ) {
  139.         helpx_(&a.lo);
  140.       }
  141.       if(a.hi > 0 ){
  142.         helpx_(&a.hi);
  143.       }
  144.       return(a);
  145.     }    
  146.  
  147.  
  148.     friend interval operator+(interval a, interval b){
  149.         return(help_(interval(a.lo+b.lo,a.hi+b.hi)));
  150.     };
  151.     friend interval operator-(interval a, interval b){
  152.         return(help_(interval(a.lo-b.hi,a.hi-b.lo)));
  153.     };
  154.     friend interval operator*(interval a, interval b){
  155.                 double ac,ad,bc,bd;
  156.                 ac=a.lo * b.lo;
  157.         ad=a.lo * b.hi;
  158.         bc=a.hi * b.lo;
  159.         bd=a.hi * b.hi;
  160.         return(help_(interval(MIN(ac,ad,bc,bd),
  161.                                       MAX(ac,ad,bc,bd))));
  162.     };
  163.     friend interval operator/(interval a, interval b){
  164.         if( b.lo == 0.0 || b.hi == 0.0 ){
  165.             cout << form("\n bad interval for division");
  166.         }
  167.         return( a * help_(interval(1.0 / b.hi, 1.0 / b.lo)) );
  168.     };
  169.     friend double lower_bound(interval a) {
  170.       return(a.lo);
  171.     };
  172.     friend double upper_bound(interval a){
  173.       return(a.hi);
  174.     };
  175. };
  176.  
  177.  
  178.  
  179. main()
  180. {
  181. int mode = 1;
  182. interval x(1.0, 2.0);
  183. interval y(4.0, 5.0);
  184. interval a(-1.0,1.0);
  185. interval b,z;
  186.  
  187.         setmode_(&mode);
  188.         a.print();
  189.         print_bin(lower_bound(a));
  190.         print_bin(upper_bound(a));
  191.         b=help_(a);
  192.  
  193.         print_bin(lower_bound(b));
  194.         print_bin(upper_bound(b));
  195.         print_bin(-lower_bound(b));
  196.         b.print();
  197.     z = x + y;
  198.     z.print();
  199.     z = x - y;
  200.     z.print();
  201.     z = x * y;
  202.     z.print();
  203.     z = x / y;
  204.     z.print();
  205. }
  206.  
  207.