SUBROUTINE ADD (A,B,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C***********************************************************************
C
C Called by --
C
C  any routine requiring interval addition.
C
C***********************************************************************
C
C Function --
C
C  This routine adds the interval A and the interval A.  It
C  simulates directed roundings with the routine RNDOUT;  the interval
C  result should contain the interval which would have been obtained
C  with exact interval arithmetic.  However, in general it will not
C  be the smallest possible machine-representable such containing
C  interval.  See the documentation in subroutine RNDOUT for more
C  detailed information.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), B(2), RESULT(2)
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C  A        is the first operand to the addition.
C           (INPUT)
C
C  B        is the second operand to the addition.
C           (INPUT)
C
C  RESULT   is the interval-arithmetic sum of A and B.
C           (OUTPUT)
C
C***********************************************************************
C
C Internal variable declarations --
C
      LOGICAL RNDDWN, RNDUP
C
C***********************************************************************
C
C Internal variable descriptions --
C
C RNDDWN    is set to .TRUE. if RNDOUT has to round down, and is set
C           to .FALSE. otherwise.
C
C RNDUP     is set to .TRUE. if RNDOUT has to round up, and is set
C           to .FALSE. otherwise.
C
C***********************************************************************
C
C Common block declarations -- none
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines -- none
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C RNDOUT
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C I/O functions --  none
C
C***********************************************************************
C
C Internal constant declarations -- none
C
C***********************************************************************
C
C Beginning of executable statements --
C
      RNDDWN = (A(1).NE.0D0).AND.(B(1).NE.0D0)
      RNDUP = (A(2).NE.0D0).AND.(B(2).NE.0D0)
C
      RESULT(1) = A(1) + B(1)
      RESULT(2) = A(2) + B(2)
C
      CALL RNDOUT(RESULT,RNDDWN,RNDUP)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE MULT(A,B,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C***********************************************************************
C
C Called by --
C
C  Any routine requiring interval multiplication.
C
C***********************************************************************
C
C Function --
C
C  This routine multiplies the interval A and the interval B.  It
C  simulates directed roundings with the routine RNDOUT;  the interval
C  result should contain the interval which would have been obtained
C  with exact interval arithmetic.  However, in general it will not
C  be the smallest possible machine-representable such containing
C  interval.  See the documentation in subroutine RNDOUT for more
C  detailed information.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), B(2), RESULT(2)
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C  A        is the first operand to the multiplication.
C           (INPUT)
C
C  B        is the second operand to the multiplication.
C           (INPUT)
C
C  RESULT   is the interval-arithmetic sum of A and B.
C           (OUTPUT)
C
C***********************************************************************
C
C Internal variable declarations --
C
        DOUBLE PRECISION R1, R2, T1(2), T2(2)
C
C***********************************************************************
C
C Internal variable descriptions --
C
C RNDDWN    is set to .TRUE. if RNDOUT has to round down, and is set
C           to .FALSE. otherwise.
C
C RNDUP     is set to .TRUE. if RNDOUT has to round up, and is set
C           to .FALSE. otherwise.
C
C***********************************************************************
C
C Common block declarations -- none
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines --
C
C  MAX, MIN
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  SCLMLT
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C I/O functions --  none
C
C***********************************************************************
C
C Internal constant declarations -- none
C
C***********************************************************************
C
C Beginning of executable statements --
C
      R1 = A(1)
      R2 = A(2)
      T1(1) = B(1)
      T1(2) = B(2)
C
      CALL SCLMLT(R1,T1,T2)
      CALL SCLMLT(R2,T1,T1)
      RESULT(1) = MIN(T1(1),T2(1))
      RESULT(2) = MAX(T1(2),T2(2))
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE POWER(A,N,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C***********************************************************************
C
C Called by --
C
C  Any routine requiring computation of a positive integer power of
C  an interval.
C
C***********************************************************************
C
C Function --
C
C  This routine computes the N-th power, of the interval A,
C  where N is a nonnegative integer.  It simulates directed
C  roundings with the routine RNDOUT;  the interval result
C  should contain the interval which would have been obtained
C  with exact interval arithmetic.  However, in general it will not
C  be the smallest possible machine-representable such containing
C  interval.  See the documentation in subroutine RNDOUT for more
C  detailed information.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), RESULT(2)
      INTEGER N
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C  A        is the base (an interval).
C           (INPUT)
C
C  N        is the power (a nonnegative integer).
C           (INPUT)
C
C  RESULT   is the interval-arithmetic value of A ** N.
C           (OUTPUT)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION B(2)
      INTEGER M
      LOGICAL RNDDWN, RNDUP
C
C***********************************************************************
C
C Internal variable descriptions --
C
C  B        is a temporary interval value.
C
C RNDDWN    is set to .TRUE. if RNDOUT has to round down, and is set
C           to .FALSE. otherwise.
C
C RNDUP     is set to .TRUE. if RNDOUT has to round up, and is set
C           to .FALSE. otherwise.
C
C***********************************************************************
C
C Common block declarations -- none
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines --
C
C  ABS, MAX, MOD
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  MULT, RNDOUT
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C I/O functions --  none
C
C***********************************************************************
C
C Internal constant declarations -- none
C
C***********************************************************************
C
C Beginning of executable statements --
C
C
C  If N is less than or equal to 0 then return the point interval
C  [ 1 , 1 ].
C
      IF (N.LE.0) THEN
         RESULT(1) = 1D0
         RESULT(2) = 1D0
         RETURN
      END IF
C
C     If N = 1 then return A.
C
      IF (N.EQ.1) THEN
         RESULT(1) = A(1)
         RESULT(2) = A(2)
         RETURN
      END IF
C
C  Let M be the largest even integer less than or equal to N,
C  and set B = A ** M.  This is to take advantage of the fact that
C  an interval to an even power must have a left endpoint which
C  is greater than or equal to zero.
C
      M = N
      IF (MOD(N,2).EQ.1) M = M - 1
C
      RNDDWN = .TRUE.
      RNDUP = .TRUE.
      IF (A(1).GT.0D0) THEN
         B(1) = A(1)**M
         B(2) = A(2)**M
      ELSE
         IF (A(2).LT.0) THEN
            B(1) = A(2)**M
            B(2) = A(1)**M
         ELSE
            B(1) = 0D0
            RNDDWN = .FALSE.
            B(2) = MAX(ABS(A(1)),ABS(A(2)))**M
         END IF
      END IF
      CALL RNDOUT(B,RNDDWN,RNDUP)
C
C  If N is even, then the result is B;  otherwise, it is B * A.
C
      IF (M.EQ.N) THEN
         RESULT(1) = B(1)
         RESULT(2) = B(2)
      ELSE
         CALL MULT(A,B,RESULT)
      END IF
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE RNDOUT(X,RNDDWN,RNDUP)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C Part of the generalized bisection package
C (interval arithmetic subpackage)
C
C***********************************************************************
C
C Called by -- ADD, SUB, POWER, SCLMLT, SCLADD, XDIV, XSCLSB
C
C***********************************************************************
C
C Function --
C
C  This routine is intended to simulate directed roundings in a
C  reasonably transportable way.  It is called for each elementary
C  operation involving intervals.  The endpoints of the result interval
C  are first computed with the machine's usual floating point
C  arithmetic.
C
C  If RNDDWN = .TRUE., then this routine decreases the left
C  endpoint of that approximate result by the absolute value of
C  that endpoint times a rigorous estimate for the maximum relative
C  error in an elementary operation.
C
C  If RNDUP = .TRUE., then this routine increases the right
C  endpoint of that approximate result by the absolute value of
C  that endpoint times a rigorous estimate for the maximum relative
C  error in an elementary operation.
C
C  For this routine to work properly, a machine-dependent parameter
C  must be installed in the routine SIMINI.  See the documentation in
C  that routine for details.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(2)
      LOGICAL RNDDWN, RNDUP
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C X         is the interval to be adjusted.
C           (I/O)
C
C RNDDWN    is set to .TRUE. if the left endpoint is to be adjusted, and
C           is set to .FALSE. otherwise.
C           (INPUT)
C
C RNDUP     is set to .TRUE. if the right endpoint is to be adjusted,
C           and is set to .FALSE. otherwise.
C           (INPUT)
C
C***********************************************************************
C
C Internal variable declarations -- none
C
C***********************************************************************
C
C Common block declarations --
C
      DOUBLE PRECISION MXULP, TTINY2, TOL0
      COMMON /MACH1/ MXULP, TTINY2, TOL0
C
C  This common block holds machine parameters which are set in
C  SIMINI and used here.
C
C  Variable descriptions
C
C  MXULP       (machine epsilon)
C                * (maximum error in ULP's of the floating pt. op's)
C
C  TTINY2      2 * (smallest representable positive machine number)
C                * (maximum error in ULP's of the floating pt. op's)
C
C  TOL0        TTINY2 / MXULP
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines -- none
C
C***********************************************************************
C
C Package-supplied functions and subroutines -- none
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C I/O functions -- none
C
C***********************************************************************
C
C Internal constant declarations -- none
C
C***********************************************************************
C
C Beginning of executable statements --
C
      IF (RNDDWN) THEN
         IF (X(1).GE.TOL0) THEN
            X(1) = X(1) - MXULP * X(1)
         ELSE IF (X(1).LE.-TOL0) THEN
            X(1) = X(1) + MXULP * X(1)
         ELSE IF ((X(1).GE.TTINY2).OR.(X(1).LE.0D0)) THEN
            X(1) = X(1) - TTINY2
         ELSE
            X(1) = 0D0
         END IF
      END IF
C
      IF (RNDUP) THEN
         IF (X(2).GE.TOL0) THEN
            X(2) = X(2) + MXULP * X(2)
         ELSE IF (X(2).LE.-TOL0) THEN
            X(2) = X(2) - MXULP * X(2)
         ELSE IF ((X(2).LE.-TTINY2).OR.(X(2).GE.0D0)) THEN
            X(2) = X(2) + TTINY2
         ELSE
            X(2) = 0D0
         END IF
      END IF
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SCLADD (R,A,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C***********************************************************************
C
C Called by --
C
C  Any routine requiring addition of a point value to an interval.
C
C***********************************************************************
C
C FUNCTION --
C
C  This routine adds the interval A to the point R.  It simulates
C  directed roundings with the routine RNDOUT;  the interval
C  result should contain the interval which would have been obtained
C  with exact interval arithmetic.  However, in general it will not
C  be the smallest possible machine-representable such containing
C  interval.  See the documentation in subroutine RNDOUT for more
C  detailed information.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION R, A(2), RESULT(2)
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C  R        is the point to be added to the interval.
C           (INPUT)
C
C  A        is the interval to be added to the point.
C           (INPUT)
C
C  RESULT   is the interval-arithmetic sum of R and B.
C           (OUTPUT)
C
C***********************************************************************
C
C Internal variable declarations --
C
      LOGICAL RNDDWN, RNDUP
C
C***********************************************************************
C
C Internal variable descriptions --
C
C RNDDWN    is set to .TRUE. if RNDOUT has to round down, and is set
C           to .FALSE. otherwise.
C
C RNDUP     is set to .TRUE. if RNDOUT has to round up, and is set
C           to .FALSE. otherwise.
C
C***********************************************************************
C
C Common block declarations -- none
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines -- none
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  RNDOUT
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C I/O functions --  none
C
C***********************************************************************
C
C Internal constant declarations -- none
C
C***********************************************************************
C
C Beginning of executable statements --
C
      RNDDWN = (R.NE.0D0).AND.(A(1).NE.0D0)
      RNDUP = (R.NE.0D0).AND.(A(2).NE.0D0)
C
      RESULT(1) = R + A(1)
      RESULT(2) = R + A(2)
C
      CALL RNDOUT(RESULT,RNDDWN,RNDUP)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SCLMLT (R,A,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C***********************************************************************
C
C Called by --
C
C  Any routine requiring multiplication of an interval and a point
C  value.
C
C***********************************************************************
C
C Function --
C
C  This routine multiplies the interval A and the point R.  It
C  simulates directed roundings with the routine RNDOUT;  the interval
C  result should contain the interval which would have been obtained
C  with exact interval arithmetic.  However, in general it will not
C  be the smallest possible machine-representable such containing
C  interval.  See the documentation in subroutine RNDOUT for more
C  detailed information.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION R, A(2), RESULT(2)
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C  R        is the point to be multiplied to the interval.
C           (INPUT)
C
C  A        is the interval to be multiplied to the point.
C           (INPUT)
C
C  RESULT   is the interval-arithmetic product R * B.
C           (OUTPUT)
C
C***********************************************************************
C
C Internal variable declarations --
C
      LOGICAL RNDDWN, RNDUP
      DOUBLE PRECISION T1, T2
C
C***********************************************************************
C
C Internal variable descriptions --
C
C RNDDWN    is set to .TRUE. if RNDOUT is to round down, and is set
C           to .FALSE. otherwise.
C
C RNDUP     is set to .TRUE. if RNDOUT is to round up, and is set
C           to .FALSE. otherwise.
C
C T1 and T2 are temporary variables.
C
C***********************************************************************
C
C Common block declarations -- none
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines --
C
C  MAX, MIN
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  RNDOUT
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C I/O functions --  none
C
C***********************************************************************
C
C Internal constant declarations -- none
C
C***********************************************************************
C
C Beginning of executable statements --
C
      IF ((R.EQ.0D0).OR.((A(1).EQ.0D0).AND.(A(2).EQ.0D0))) THEN
         RESULT(1) = 0D0
         RESULT(2) = 0D0
         RETURN
      END IF
C
      T1 = A(1)
      T2 = A(2)
      RNDDWN = .TRUE.
      RNDUP = .TRUE.
C
      IF (T1.EQ.0D0) THEN
         IF (R.LT.0D0) THEN
            RESULT(1) = R * T2
            RESULT(2) = 0D0
            RNDUP = .FALSE.
         ELSE
            RESULT(1) = 0D0
            RESULT(2) = R * T2
            RNDDWN = .FALSE.
         END IF
      ELSE
         IF (T2.EQ.0D0) THEN
            IF (R.LT.0D0) THEN
               RESULT(1) = 0D0
               RESULT(2) = R * T1
               RNDDWN = .FALSE.
             ELSE
               RESULT(1) = R * T1
               RESULT(2) = 0D0
               RNDUP = .FALSE.
             END IF
         ELSE
            T1 = R * T1
            T2 = R * T2
            RESULT(1) = MIN(T1,T2)
            RESULT(2) = MAX(T1,T2)
         END IF
      END IF
C
      CALL RNDOUT(RESULT,RNDDWN,RNDUP)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SIMINI
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C Part of the generalized bisection package
C (interval arithmetic subpackage)
C
C***********************************************************************
C
C Called by -- GENBIS
C
C***********************************************************************
C
C Function --
C
C  This routine sets certain machine parameters used to simulate
C  directed roundings in a reasonably transportable way.  In
C  particular, it sets the amount by which to decrease the left endpoint
C  and increase the right endpoint of an interval computed using usual
C  floating point arithmetic to guarantee that the resulting interval
C  will contain the result which would have been obtained with true
C  interval arithmetic.
C
C  This routine assumes that the four elementary floating point
C  operations and taking an integer power will give results with a
C  maximum error of one ULP (unit in the last place).  If this is not
C  so, change the value of MAXERR in the data statement below to the
C  maximum number of ULP's by which a floating point result can differ
C  from the true result (for '+', '-', '*', '/', and '** N').
C
C  When determining the maximum error of the result A op B, where
C  A and B are floating point numbers, we assume that A and B are
C  represented exactly.  For example, if A and B are almost equal,
C  then it is not unreasonable to assume that A - B, where the
C  subtraction is a floating point subtraction, is within a few
C  units of the last place of the true result.
C
C  If SIMINI is installed correctly, then the conclusions this
C  package prints out will have mathematical rigor.
C
C***********************************************************************
C
C Common block declarations --
C
      DOUBLE PRECISION MXULP, TTINY2, TOL0
      COMMON /MACH1/ MXULP, TTINY2, TOL0
C
C  This common block holds machine parameters which are set here
C  and used in RNDOUT.
C
C  Variable descriptions
C
C  MXULP       (machine epsilon)
C                * (maximum error in ULP's of the floating pt. op's)
C
C  TTINY2      2 * (smallest representable positive machine number)
C                * (maximum error in ULP's of the floating pt. op's)
C
C  TOL0        TTINY2 / MXULP
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines -- DBLE
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  D1MACH (the SLATEC routine for machine constants)
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C I/O functions -- none
C
C***********************************************************************
C
C Internal constant declarations --
C
      INTEGER MAXERR
      DATA MAXERR/1/
C
C***********************************************************************
C
C Internal constant descriptions --
C
C  MAXERR      is the maximum number of ULP's (units in the last
C              place) by which a result of one of the floating point
C              operations (+, -, *, /, ** N) can differ from the
C              true result.  (See explanation above.)
C
C********** WARNING: The value of MAXERR is machine dependent and
C                    must be manually set.
C
C***********************************************************************
C
C Beginning of executable statements --
C
      MXULP = DBLE(MAXERR) * D1MACH(4)
      TTINY2 = 2D0 * DBLE(MAXERR) * D1MACH(1)
      TOL0 = TTINY2 / MXULP
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SUB (A,B,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C***********************************************************************
C
C Called by --
C
C  Any routine requiring interval subtraction.
C
C***********************************************************************
C
C Function --
C
C  This routine subtracts the interval B from the interval A.  It
C  simulates directed roundings with the routine RNDOUT;  the interval
C  result should contain the interval which would have been obtained
C  with exact interval arithmetic.  However, in general it will not
C  be the smallest possible machine-representable such containing
C  interval.  See the documentation in subroutine RNDOUT for more
C  detailed information.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), B(2), RESULT(2)
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C  A        is the first operand to the subtraction.
C           (INPUT)
C
C  B        is the second operand to the subtraction
C           (INPUT)
C
C  RESULT   is the interval-arithmetic value of A - B.
C           (OUTPUT)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TA1, TA2, TB1, TB2
      LOGICAL RNDDWN, RNDUP
C
C***********************************************************************
C
C Internal variable descriptions --
C
C TA1, TA2, TB1, and TB2 are temporaries.
C
C RNDDWN    is set to .TRUE. if RNDOUT has to round down, and is set
C           to .FALSE. otherwise.
C
C RNDUP     is set to .TRUE. if RNDOUT has to round up, and is set
C           to .FALSE. otherwise.
C
C***********************************************************************
C
C Common block declarations -- none
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines -- none
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C RNDOUT
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C I/O functions --  none
C
C***********************************************************************
C
C Internal constant declarations -- none
C
C***********************************************************************
C
C Beginning of executable statements --
C
      TA1 = A(1)
      TA2 = A(2)
      TB1 = B(1)
      TB2 = B(2)
C
      RNDDWN = (TB2.NE.0D0)
      RNDUP = (TB1.NE.0D0)
C
      RESULT(1) = TA1 - TB2
      RESULT(2) = TA2 - TB1
C
      CALL RNDOUT(RESULT,RNDDWN,RNDUP)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE XDIV (XCASE,A,B,R1,R2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C***********************************************************************
C
C Called by -- DVSBIN
C
C***********************************************************************
C
C Function --
C
C  This routine performs the extended interval arithmetic division
C  A / B, and places the result(s) in R1 (and R2).  The value of
C  XCASE is set according to whether the result is one finite
C  interval, one or two semi-infinite intervals, or the real line.
C  This routine does not use directed roundings, but the true interval
C  result should be contained in the returned interval(s) because the
C  routine RNDOUT widens the floating point intervals slightly.
C  This scheme does not in general produce the smallest machine
C  representable intervals containing the result, but is reasonably
C  portable.  For more detailed information, see the documentation in
C  the routine RNDOUT.
C
C***********************************************************************
C
C Argument declarations --
C
      INTEGER XCASE
      DOUBLE PRECISION A(2), B(2), R1(2), R2(2)
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C XCASE     is set to
C
C              (1)   if the result is a single finite interval
C                    (stored in R1);
C
C              (2)   if the result is a single semi-infinite interval
C                    (stored in R1) of the form
C                    [ finite , + infinity ];
C
C              (3)   if the result is a single semi-infinite interval
C                    (stored in R1) of the form
C                    [ - infinity , finite ];
C
C              (4)   if the result is the real line;
C
C              (5)   if the result is a union of two semi-infinite
C                    intervals R1 = [ - infinity, finite ] and
C                    R2 = [ finite, + infinity ].
C
C           (OUTPUT)
C
C A         is the numerator of the quotient.
C           (INPUT)
C
C B         is the denominator of the quotient.
C           (INPUT)
C
C R1        is the interval quotient (or one element of it, when it
C           consists of more than one semi-infinite interval).
C           (OUTPUT)
C
C R2        is the second element of the interval quotient if it is
C           a union of two semi-infinite intervals.
C           (OUTPUT)
C
C           When R1 or R2 contains semi-infinite intervals,
C           only the finite endpoints are stored.
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION C(2)
      LOGICAL RNDDWN, RNDUP
C
C***********************************************************************
C
C Internal variable descriptions --
C
C  C        is a temporary interval.
C
C RNDDWN    is set to .TRUE. if RNDOUT has to round down, and is set
C           to .FALSE. otherwise.
C
C RNDUP     is set to .TRUE. if RNDOUT has to round up, and is set
C           to .FALSE. otherwise.
C
C***********************************************************************
C
C Common block declarations -- none
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines -- none
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  MULT, RNDOUT
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C I/O functions -- none
C
C***********************************************************************
C
C Internal constant declarations -- none
C
C***********************************************************************
C
C Beginning of executable statements --
C
      RNDDWN = .TRUE.
      RNDUP = .TRUE.
C
C  Do usual interval division if zero is not in the denominator.
C
      IF ((B(1).GT.0D0).OR.(B(2).LT.0D0)) THEN
         XCASE = 1
         C(1) = 1D0/B(2)
         C(2) = 1D0/B(1)
C
         CALL RNDOUT(C,RNDDWN,RNDUP)
C
         CALL MULT(A,C,R1)
         RETURN
      END IF
C
C  If the denominator is equal to zero or if the numerator contains
C  zero, then set the result to the real line.
C
      IF ((B(1).EQ.B(2)).OR.((A(1).LT.0D0).AND.(A(2).GT.0D0))) THEN
         XCASE = 4
         RETURN
      END IF
C
C  This is the case when the left endpoint of the denominator is zero.
C
      IF (B(1).EQ.0D0) THEN
         IF (A(1).GE.0D0) THEN
            XCASE = 2
            RNDDWN = (A(1).NE.0D0)
            RNDUP = .FALSE.
            R1(1) = A(1) / B(2)
            CALL RNDOUT(R1,RNDDWN,RNDUP)
            RETURN
         ELSE
            XCASE = 3
            RNDUP = (A(2).NE.0D0)
            RNDDWN = .FALSE.
            R1(2) = A(2) / B(2)
            CALL RNDOUT(R1,RNDDWN,RNDUP)
            RETURN
         END IF
      END IF
C
C  This is the case when the right endpoint of the denominator is zero.
C
      IF (B(2).EQ.0D0) THEN
         IF (A(2).LE.0D0) THEN
            XCASE = 2
            RNDDWN = (A(2).NE.0D0)
            RNDUP = .FALSE.
            R1(1) = A(2) / B(1)
            CALL RNDOUT(R1,RNDDWN,RNDUP)
            RETURN
         ELSE
            XCASE = 3
            RNDUP = (A(1).NE.0D0)
            RNDDWN = .FALSE.
            R1(2) = A(1) / B(1)
            CALL RNDOUT(R1,RNDDWN,RNDUP)
            RETURN
         END IF
      END IF
C
C  This is the case when the denominator contains zero but the
C  numerator does not.  The result is a union of two semi-infinite
C  intervals.
C
      IF (A(1).GE.0D0) THEN
         RNDUP = (A(1).NE.0D0)
         C(2) = A(1) / B(1)
         RNDDWN = RNDUP
         C(1) = A(1) / B(2)
      ELSE
         RNDUP = (A(2).NE.0D0)
         C(2) = A(2) / B(2)
         RNDDWN = RNDUP
         C(1) = A(2) / B(1)
      END IF
      IF (C(1).GT.C(2)) THEN
         XCASE = 5
         R1(2) = C(2)
         R2(1) = C(1)
      ELSE
         XCASE = 4
      END IF
      CALL RNDOUT(R2,RNDDWN,RNDUP)
C
      CALL RNDOUT(R1,RNDDWN,RNDUP)
      RETURN
C
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE XINT(CASE,A,B,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C***********************************************************************
C
C Called by -- DVSBIN
C
C***********************************************************************
C
C Function --
C
C  This routine intersects a finite interval A with an extended-value
C  interval B whose type is given by CASE.  It places the resulting
C  interval in RESULT, if it is not null;  the routine resets the
C  variable CASE to indicate the type of RESULT.
C
C  All intervals are represented by 2-vectors.
C
C***********************************************************************
C
C Argument declarations --
C
      INTEGER CASE
      DOUBLE PRECISION A(2), B(2), RESULT(2)
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C CASE      indicates the type of A on entry, and the type of
C           RESULT on return.  It is set to
C
C              (0)   if the quantity is the empty interval;
C
C              (1)   if the quantity is a single finite interval;
C
C              (2)   if the quantity is a single semi-infinite interval
C                    of the form [finite, +infinity];
C
C              (3)   if the quantity is a single semi-infinite interval
C                    of the form [ - infinity , finite ];
C
C              (4)   if the quantity is the real line.
C
C           (I/O)
C
C A         is a finite interval.
C           (INPUT)
C
C B         is an extended-value interval (possibly infinite or
C           semi-infinite) whose type on entry is given by CASE.
C           (INPUT)
C
C RESULT    is set to the intersection of A and B on return, if
C           that intersection is not null.
C           (OUTPUT)
C
C***********************************************************************
C
C Internal variable declarations -- none
C
C***********************************************************************
C
C Common block declarations -- none
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines --
C
C  MAX, MIN
C
C***********************************************************************
C
C Package-supplied functions and subroutines -- none
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C I/O functions -- none
C
C***********************************************************************
C
C Internal constant declarations -- none
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Finite case
C
      IF (CASE.EQ.1) THEN
         IF ((A(2).LT.B(1)).OR.(A(1).GT.B(2))) THEN
            CASE = 0
            RETURN
         ELSE
            RESULT(1) = MAX(A(1),B(1))
            RESULT(2) = MIN(A(2),B(2))
            RETURN
         END IF
      END IF
C
C  Semi-infinite case
C
      IF (CASE.EQ.2) THEN
         IF (A(2).LT.B(1)) THEN
            CASE = 0
            RETURN
         ELSE
            CASE = 1
            RESULT(1) = MAX(A(1),B(1))
            RESULT(2) = A(2)
            RETURN
         END IF
      END IF
C
      IF (CASE.EQ.3) THEN
         IF (A(1).GT.B(2)) THEN
            CASE = 0
            RETURN
         ELSE
            CASE = 1
            RESULT(1) = A(1)
            RESULT(2) = MIN(A(2),B(2))
            RETURN
         END IF
      END IF
C
C  Infinite case
C
      IF (CASE.EQ.4) THEN
         CASE = 1
         RESULT(1) = A(1)
         RESULT(2) = A(2)
         RETURN
      END IF
C
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE XSCLSB(XCASE,R,A,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C***********************************************************************
C
C Called by -- DVSBIN
C
C***********************************************************************
C
C Function --
C
C  This routine subtracts an extended-value interval A, whose type
C  is indicated by XCASE, from a point R.  The result is an extended-
C  value interval RESULT, whose type is stored in XCASE.
C  This routine does not use directed roundings, but the true interval
C  result should be contained in the returned interval(s) because the
C  routine RNDOUT widens the floating point intervals slightly.
C  This scheme does not in general produce the smallest machine
C  representable intervals containing the result, but is reasonably
C  portable.  For more detailed information, see the documentation in
C  the routine RNDOUT.
C
C***********************************************************************
C
C Argument declarations --
C
      INTEGER XCASE
      DOUBLE PRECISION R, A(2), RESULT(2)
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C XCASE     indicates the type of A on entry, and the type of
C           RESULT on return.  It is set to
C
C              (1)   if the quantity is a single finite interval;
C
C              (2)   if the quantity is a single semi-infinite interval
C                    of the form [finite, +infinity];
C
C              (3)   if the quantity is a single semi-infinite interval
C                    of the form [ - infinity , finite ];
C
C              (4)   if the quantity is the real line;
C
C           (I/O)
C
C R         is the point from which the interval is subtracted.
C           (INPUT)
C
C A         is the extended-value interval to be subtracted.
C           On entry, its type is given by XCASE.
C           (INPUT)
C
C RESULT    is the extended-value interval result  R - A.
C           On return, its type is given by XCASE.
C           (OUTPUT)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION T1, T2
      LOGICAL RNDDWN, RNDUP
C
C***********************************************************************
C
C Internal variable descriptions --
C
C  T1, and T2 are temporary storage.
C
C RNDDWN    is set to .TRUE. if RNDOUT has to round down, and is set
C           to .FALSE. otherwise.
C
C RNDUP     is set to .TRUE. if RNDOUT has to round up, and is set
C           to .FALSE. otherwise.
C
C***********************************************************************
C
C Common block declarations -- none
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines -- none
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  RNDOUT
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C I/O functions -- none
C
C***********************************************************************
C
C Internal constant declarations -- none
C
C***********************************************************************
C
C Beginning of executable statements --
C
      T1 = A(1)
      T2 = A(2)
      RNDDWN = (T2.NE.0D0)
      RNDUP = (T1.NE.0D0)
C
C     FINITE CASE
C
      IF (XCASE.EQ.1) THEN
         RESULT(1) = R - T2
         RESULT(2) = R - T1
         CALL RNDOUT(RESULT,RNDDWN,RNDUP)
         RETURN
      END IF
C
C  The semi-infinite cases follow.
C
      IF (XCASE.EQ.2) THEN
         XCASE = 3
         RESULT(2) = R - T1
         CALL RNDOUT(RESULT,.FALSE.,RNDUP)
         RETURN
      END IF
      IF (XCASE.EQ.3) THEN
         XCASE = 2
         RESULT(1) = R - T2
         CALL RNDOUT(RESULT,RNDDWN,.FALSE.)
         RETURN
      END IF
C
C  If we get this far, then the interval is the real line, and
C  hence so is the result. We therefore return without changing XCASE.
C
      RETURN
      END