Usenet 1994 January/usenetsourcesnewsgroupsinfomagicjanuary1994.iso/sources/misc/volume8/interval.c++

C/C++ Source or Header | 1989-10-28 | 11KB | 421 lines

Newsgroups: comp.sources.misc
subject: v08i098: Interval Arithmetic in C++
From: [e-mail] (Brandon S. Allbery - comp.sources.misc)
Reply-To: @uunet.uu.net:[e-mail]

Posting-number: Volume 8, Issue 98
Submitted-by: @uunet.uu.net:[e-mail]
Archive-name: interval.c++

    Enclosed is a shar archive for a simple interval
arithmetic package written in C++.  There is no
restriction on use or re-distribution.  Comments and
improvements may be addressed to:
  From U.K.: [e-mail]
  Elsewhere: [e-mail]
------------------------ cut here ----------------------------
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create:
#    Makefile
#    README
#    Tinterval.c
#    interval.c
#    interval.h
# This archive created: Thu Aug 24 13:13:29 1989
export PATH; PATH=/bin:/usr/bin:$PATH
echo shar: "extracting 'Makefile'" '(715 characters)'
if test -f 'Makefile'
then
    echo shar: "will not over-write existing file 'Makefile'"
else
cat << \SHAR_EOF > 'Makefile'
HDRS          = interval.h
CC          = CC
LINKER          = CC
MAKEFILE      = Makefile
OBJS          = Tinterval.o \
        interval.o
PRINT          = pr
PROGRAM          = Tinterval
SRCS          = Tinterval.c \
        interval.c

all:        $(PROGRAM)

$(PROGRAM):     $(OBJS) $(LIBS)
        @echo -n "Linking $(PROGRAM) ... "
        @$(LINKER) $(LDFLAGS) $(OBJS) $(LIBS) -o $(PROGRAM)
        @echo "done"

clean:;        @rm -f $(OBJS) core a.out

depend:;    @mkmf -f $(MAKEFILE) PROGRAM=$(PROGRAM) DEST=$(DEST)

index:;        @ctags -wx $(HDRS) $(SRCS)

print:;        @$(PRINT) $(HDRS) $(SRCS)

tags:           $(HDRS) $(SRCS); @ctags $(HDRS) $(SRCS)

###
Tinterval.o: /usr/include/CC/stdio.h /usr/include/CC/string.h interval.h
interval.o: /usr/include/CC/stdio.h interval.h
SHAR_EOF
if test 715 -ne "`wc -c < 'Makefile'`"
then
    echo shar: "error transmitting 'Makefile'" '(should have been 715 characters)'
fi
fi
echo shar: "extracting 'README'" '(498 characters)'
if test -f 'README'
then
    echo shar: "will not over-write existing file 'README'"
else
cat << \SHAR_EOF > 'README'
    This is a simple implementation of an interval
arithmetic package I wrote while trying to learn C++.
The accuracy is dubious (since proper attention has not
been paid to rounding) and there is no attempt to recover
from errors.  The result is not too robust, but fun, useful
and (in my case), instructive.
    Be sure to modify the makefile to use your c++ compiler.
You may also need to change the pathnames of include files
(e.g., <CC/stdio.h>) in interval.c and Tinterval.c for your
system.
SHAR_EOF
if test 498 -ne "`wc -c < 'README'`"
then
    echo shar: "error transmitting 'README'" '(should have been 498 characters)'
fi
fi
echo shar: "extracting 'Tinterval.c'" '(1020 characters)'
if test -f 'Tinterval.c'
then
    echo shar: "will not over-write existing file 'Tinterval.c'"
else
cat << \SHAR_EOF > 'Tinterval.c'
// A simple test program for the interval class

#include <CC/stdio.h>
#include <CC/string.h>
#include "interval.h"

void Usage(char *name)
{
    fprintf(stderr, "usage: %s real real +|-|*|/ real real\n", name);
}

main(int argc, char **argv)
{
    float a,b,c,d;

    if (argc != 6) {
    Usage(argv[0]);
    exit(1);
    }

    a = atof(argv[1]);
    b = atof(argv[2]);
    c = atof(argv[4]);
    d = atof(argv[5]);
    {
    interval I1(a,b);
    interval I2(c,d);
    interval Result, Result2;

    switch (argv[3][0]) {
    case '+': Result = I1 + I2; break;
    case '-': Result = I1 - I2; break;
    case '*': Result = I1 * I2; break;
    case '/': Result = I1 / I2; break;
    default:
        Usage(argv[0]);
        exit(2);
    }

    printf("[%f, %f] %s [%f, %f] == [%f, %f]\n",
        I1.lo(), I1.hi(), argv[3], I2.lo(), I2.hi(),
        Result.lo(), Result.hi());

    // Now check out += and integer-to-interval conversion
    Result2 = Result;
    Result2 += 4;
    printf("[%f, %f] + 4 == [%f, %f]\n",
        Result.lo(), Result.hi(), Result2.lo(), Result2.hi());
    }
}
SHAR_EOF
if test 1020 -ne "`wc -c < 'Tinterval.c'`"
then
    echo shar: "error transmitting 'Tinterval.c'" '(should have been 1020 characters)'
fi
fi
echo shar: "extracting 'interval.c'" '(4825 characters)'
if test -f 'interval.c'
then
    echo shar: "will not over-write existing file 'interval.c'"
else
cat << \SHAR_EOF > 'interval.c'
// This software is in the public domain.
// Permission is granted to use this software for any
// purpose commercial or non-commercial.
// The implementer/distributer assume no liability for
// any damages, incidental, consequential or otherwise,
// arising from the use of this software.

// Implementation of class interval
//
// Most operators are trivial and could be inlined.
//
// Multiply is more interesting.  The complicated testing
// may be more time consuming than the floating-point
// operations it attempts to avoid.
// How could this testing be streamlined?
//
// Note that divide blurts out error/warning messages using fprintf :-(
// Strictly speaking, division is undefined if the divisor interval
// contains zero.  This code is sloppy in allowing it, but it warns.
//
// No attempt is made to control rounding so accuracy is dubious.
//
// No attempt is made to detect/recover from overflow or underflow.
// This is a serious deficiency.  Is there a 'portable' around it?

#include <CC/stdio.h>
#include "interval.h"

interval interval::operator+(interval I)
{
    interval temp(lo_bound, hi_bound);
    return temp += I;
}

interval interval::operator-(interval I)
{
    interval temp(lo_bound, hi_bound);
    return temp -= I;
}

interval interval::operator*(interval I)
{
    interval temp(lo_bound, hi_bound);
    return temp *= I;
}

interval interval::operator/(interval I)
{
    interval temp(lo_bound, hi_bound);
    return temp /= I;
}

interval &interval::operator+=(interval I)
{
    lo_bound += I.lo_bound;
    hi_bound += I.hi_bound;
    return *this;
}

interval &interval::operator-=(interval I)
{
    lo_bound -= I.lo_bound;
    hi_bound -= I.hi_bound;
    return *this;
}

// Multiply is a bit complicated.
// A naive version simply takes the minimum of all combinations
// as the new lo_bound and the maximum as the new hi_bound.
// But this involves four double precision floating multiples.
// A quick examination of the signs of the intervals reveals
// that there are nine possible cases of which only one
// requires all four multiplications.  The other cases only
// require two multiplications.  We assume that it is faster to
// compute and dispatch to the right case than to perform
// the extra multiplies.  If not, use the naive scheme.
//
// The cases can be described as follows:
// Each interval is either non-negative (lo bound >= 0),
// non-positive (hi bound <= 0), or crosses(includes) zero,
// including some numbers on either side.
// If we label the intervals, A and B, corresponding to
// "this" and "I", we get the following matrix:
//
//        Bcross        Bnonneg        Bnonpos
//
//    Across
//
//    Anonneg
//
//    Anonpos

interval &interval::operator*=(interval I)
{
    enum possibility {
    AcrossBcross, AcrossBnn, AcrossBnp,
    AnnBcross, AnnBnn, AnnBnp,
    AnpBcross, AnpBnn, AnpBnp
    } choice;

    if (lo_bound >= 0.0) choice = AnnBcross;
    else if (hi_bound <= 0.0) choice = AnpBcross;
    else choice = AcrossBcross;

    if (I.lo_bound >= 0.0) choice += 1;
    else if (I.hi_bound <= 0.0) choice += 2;

    switch (choice)
    {
    case AcrossBcross: {
            double  HL = hi_bound*I.lo_bound,
                HH = hi_bound*I.hi_bound,
                LL = lo_bound*I.lo_bound,
                LH = lo_bound*I.hi_bound;
            lo_bound = HL<LH?HL:LH;
            hi_bound = LL>HH?LL:HH;
        }
        break;

    case AcrossBnn:
        lo_bound *= I.hi_bound;
        hi_bound *= I.hi_bound;
        break;

    case AcrossBnp: {
        double new_hi_bound = lo_bound * I.lo_bound;

        lo_bound = hi_bound*I.lo_bound;
        hi_bound = new_hi_bound;
        }
        break;

    case AnnBcross:
        lo_bound = hi_bound*I.lo_bound;
        hi_bound *= I.hi_bound;
        break;

    case AnnBnn:
        lo_bound *= I.lo_bound;
        hi_bound *= I.hi_bound;
        break;

    case AnnBnp: {
        double new_hi_bound = lo_bound * I.hi_bound;

        lo_bound = hi_bound*I.lo_bound;
        hi_bound = new_hi_bound;
        }
        break;

    case AnpBcross:
        hi_bound = lo_bound*I.lo_bound;
        lo_bound *= I.hi_bound;
        break;

    case AnpBnn:
        lo_bound *= I.hi_bound;
        hi_bound *= I.lo_bound;
        break;

    case AnpBnp: {
        double new_hi_bound = lo_bound * I.lo_bound;

        lo_bound = hi_bound*I.hi_bound;
        hi_bound = new_hi_bound;
        }
        break;

    default:
        fprintf(stderr,
            "Bad interval (low bound > hi bound) [%f, %f] * [%f, %f]\n",
            lo_bound, hi_bound, I.lo_bound, I.hi_bound);
        lo_bound = hi_bound = 0.0;
        break;
    }

    return *this;
}

interval &interval::operator/=(interval I)
{
    if (I.lo_bound == 0.0 && I.hi_bound == 0)
    fprintf(stderr, "Error: Division by zero attempted - [%f,%f]/[f,%f]\n",
        lo_bound, hi_bound, I.lo_bound, I.hi_bound);
    else {
    if (I.lo_bound < 0 && I.hi_bound > 0)
        fprintf(stderr, "Warning: Division by zero possible - [%f,%f]/[%f,%f]\n",
        lo_bound, hi_bound, I.lo_bound, I.hi_bound);
    interval inv(1.0/I.lo_bound, 1.0/I.hi_bound);
    return inv *= *this;
    }
}
SHAR_EOF
if test 4825 -ne "`wc -c < 'interval.c'`"
then
    echo shar: "error transmitting 'interval.c'" '(should have been 4825 characters)'
fi
fi
echo shar: "extracting 'interval.h'" '(1684 characters)'
if test -f 'interval.h'
then
    echo shar: "will not over-write existing file 'interval.h'"
else
cat << \SHAR_EOF > 'interval.h'
// Class interval
//    A definition of intervals for interval arithmetic.
//
//    The representation of the interval is two doubles
//    but conversion operators are defined for int and float.
//
//    Error handling is weak.  In an IEEE compliant implementation,
//    quiet NANs could be used to indicate bounds that were
//    invalid.  Otherwise, more state (and checking) would
//    need to be added.

// To get the definition of MAXDOUBLE
#include <values.h>

class interval {
    double lo_bound, hi_bound;
public:
    double lo(void)        { return lo_bound; }
    double hi(void)        { return hi_bound; }
    double width(void)        { return hi_bound - lo_bound; }
    interval()            { lo_bound = -MAXDOUBLE; hi_bound = MAXDOUBLE; }
    interval(double l, double h){ lo_bound = l<h?l:h; hi_bound = l<h?h:l; }
    interval(double f)        { lo_bound = f; hi_bound = f; }
    interval(float l, float h)    { lo_bound = l<h?l:h; hi_bound = l<h?h:l; }
    interval(float f)        { lo_bound = f; hi_bound = f; }
    interval(int l, int h)    { lo_bound = l<h?l:h; hi_bound = l<h?h:l; }
    interval(int i)        { lo_bound = i; hi_bound = i; }
    interval operator+(interval I);
    interval operator-(interval I);
    interval operator*(interval I);
    interval operator/(interval I);
    interval &operator+=(interval I);
    interval &operator-=(interval I);
    interval &operator*=(interval I);
    interval &operator/=(interval I);
    int contains(interval I)    { return lo_bound <= I.lo_bound &&
                     hi_bound >= I.hi_bound; }
    int overlaps(interval I)    { return lo_bound <= I.hi_bound &&
                     hi_bound >= I.lo_bound; }
    int equal(interval I)    { return lo_bound == I.lo_bound &&
                     hi_bound == I.hi_bound; }
};
SHAR_EOF
if test 1684 -ne "`wc -c < 'interval.h'`"
then
    echo shar: "error transmitting 'interval.h'" '(should have been 1684 characters)'
fi
fi
exit 0
#    End of shell archive
#
#
#