/* 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 architecture. 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
#include
#include
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;i0.
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();
}