Photo CD Demo 1/Demo.bin/gems/gemsiii/interval.c

C/C++ Source or Header | 1992-03-16 | 5KB | 207 lines

/* The following code implements the four basic interval arithmetic
operations +, -, * and / within the C++ class interval as provided
by the AT & T C++ Release 2.0 for the Sparc architechture. There is
no checking that the intervals are valid intervals, i.e. if the
interval is [lo.hi] then it is valid iff lo <= hi.

It is assumed that a mechanism exists to direct the result of a
floating point operation so that it is chopped. This corresponds
to "tozero" in ANSI/IEEE Std. 754-1985 and it is performed by the
routine
               ieee_flags
in this code.

The first executable statement of any program using this code must
therefore be
               setmode(1);
or a corresponding call in in the architecture used in order that
the rounding mode will be the chopped mode which is assumed in the
class. Changes have to be made to the routine help_ if another
rounding mode is used.
*/


#define    min(a, b)        (((a) < (b)) ? (a) : (b))
#define    MIN(a, b, c, d)        min(min(min(a, b), c), d)
#define    max(a, b)        (((a) > (b)) ? (a) : (b))
#define    MAX(a, b, c, d)        max(max(max(a, b), c), d)

#include <stdio.h>
#include <stream.h>
#include <sys/ieeefp.h>

extern "C" {
    int ieee_flags(char *, char *, char *in, char **);
      }
void setmode_(int *mode)
{
char *out;

    ieee_flags("set", "direction", *mode ? "tozero" : "nearest", &out);
    ieee_flags("get", "direction", NULL, &out);
}



/*

This routine adds 1 to the last digit of a double data item which is here
assumed to have 8 bytes where

seeeeeee eeeemmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmmm mmmmmmmm

s=sign of number
e=exponent            01111111111 is an exponent of 0
m=mantissa bits with assumed leading 1

The algorithm works by adding one to the last two bytes of the mantissa.
If the result of this addition is zero then it overflows to the next two
bytes. This means that the next two bytes must be checked for overflow
until the mantissa bits are exhausted. It should then be noted that the
exponent will be automatically incremented and the mantissa will have
been set to all zeros.

No check is made to see if the exponent overflows. Negative numbers will
turn into the next larger negative number.

This routine assumes that the number of bytes in a double is and integer
multiple of the number of bytes in an unsigned integer.
*/


void helpx_(double *a)
{
unsigned char *b;
unsigned int *c,i;

    b=(unsigned char *)a;
    i=sizeof(double)-sizeof(unsigned int);
    do {
        c=(unsigned int *)(b+i);
        (*c)++;
        i-=sizeof(unsigned int);
    } while ((*c)==0 && i>=0);
}

/* A printing utility for printing the octal contents of a double
data item */

void print_bin(double a)
{
char *b;
int i,j,k;

    k=0;
    b=(char *)(&a);
    printf("%f ",a);
    for (i=0;i<sizeof(double);i++) {
        for(j=0;j<8;j++) {
            k++;
            if (k==2 || k==13) putchar(' ');
            printf("%d",(b[i]&0x80)==128);
            b[i]<<=1;
        }
    }
    printf("\n");
}

/*
The basic arithmetic operations are embedded in the class interval. In
each case the pair of double data items representing the interval [lo,hi]
is checked via the routine help_. This routine accesses helpx_ whenever
          lo<0
or
          hi>0.
In both cases the next double representable item is calculated in helpx_
and the result is returned to the calling operator. In this manner it is
guaranteed that the machine interval result will contain the interval
result that would have been computed using infinite (real) interval
arithmetic.
*/

class    interval {
        double lo,hi;


public:
    interval(double lo = 0, double hi = 0) {
        this->lo = lo;
        this->hi = hi;
    };

/* An interval printing utility */
        void print() {
        printf("[ %3lx, %3lx ]\n", lo, hi);
    };

    friend interval help_(interval a) {
      if(a.lo < 0 ) {
        helpx_(&a.lo);
      }
      if(a.hi > 0 ){
        helpx_(&a.hi);
      }
      return(a);
    }


    friend interval operator+(interval a, interval b){
        return(help_(interval(a.lo+b.lo,a.hi+b.hi)));
    };
    friend interval operator-(interval a, interval b){
        return(help_(interval(a.lo-b.hi,a.hi-b.lo)));
    };
    friend interval operator*(interval a, interval b){
                double ac,ad,bc,bd;
                ac=a.lo * b.lo;
        ad=a.lo * b.hi;
        bc=a.hi * b.lo;
        bd=a.hi * b.hi;
        return(help_(interval(MIN(ac,ad,bc,bd),
                                      MAX(ac,ad,bc,bd))));
    };
    friend interval operator/(interval a, interval b){
        if( b.lo == 0.0 || b.hi == 0.0 ){
            cout << form("\n bad interval for division");
        }
        return( a * help_(interval(1.0 / b.hi, 1.0 / b.lo)) );
    };
    friend double lower_bound(interval a) {
      return(a.lo);
    };
    friend double upper_bound(interval a){
      return(a.hi);
    };
};



main()
{
int mode = 1;
interval x(1.0, 2.0);
interval y(4.0, 5.0);
interval a(-1.0,1.0);
interval b,z;

        setmode_(&mode);
        a.print();
        print_bin(lower_bound(a));
        print_bin(upper_bound(a));
        b=help_(a);

        print_bin(lower_bound(b));
        print_bin(upper_bound(b));
        print_bin(-lower_bound(b));
        b.print();
    z = x + y;
    z.print();
    z = x - y;
    z.print();
    z = x * y;
    z.print();
    z = x / y;
    z.print();
}