C*** basicops.f
      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 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 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 CANCEL (A,B,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   March 13, 1990
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C***********************************************************************
C
C Called by --
C
C  Any routine requiring an interval addend to be canceled (removed)
C  from a previously accumulated interval sum.
C
C***********************************************************************
C
C Function --
C
C  Given an interval B, and a previously accumulated interval sum A
C  for which B was an addend, this routine returns an interval which
C  contains, and is hopefully close to, the sum of the other addends
C  for A.  Directed roundings are simulated with the routine RNDOUT;
C  the interval result should contain the interval which would have
C  been obtained with exact interval arithmetic.  However, in general
C  it will not be the smallest possible machine-representable such
C  containing interval.  See the documentation in subroutine RNDOUT for
C  more 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
      LOGICAL RNDDWN, RNDUP
C
C***********************************************************************
C
C Internal variable descriptions --
C
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 Package-supplied functions and subroutines --
C
C RNDOUT
C
C***********************************************************************
C
C Beginning of executable statements --
C
      RNDDWN = (B(1).NE.0D0)
      RNDUP = (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 IDIV (AA,BB,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   May 21, 1992
C
C Part of the interval arithmetic elementary function package.
C
C***********************************************************************
C
C Function --
C
C  This routine performs interval division.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION AA(2), BB(2), RESULT(2)
C
C***********************************************************************
C
C Internal variable declarations --

      DOUBLE PRECISION C(2)
C
C***********************************************************************
C
C  Common block declarations --
C
      DOUBLE PRECISION MAXX, MXLGM1, NEGINF, POSINF, A(2), PI(2), E(2),
     *   PI2(2), PI3(2), PI4(2), PI6(2), PI8(2), E14(2), A3(2), ONE(2),
     *   OD2F(2), OD3F(2), OD4F(2), OD5F(2), OD6F(2), OD7F(2), OD8F(2),
     *   OD9F(2), OD10F(2), OD11F(2), OD12F(2), OD14F(2), EIGHT(2),
     *   ZERO(2), TWO(2), THREE(2), FOUR(2), SXTNTH(2), NINE(2), TEN(2),
     *   TWOT7(2), SQT10(2), ESXTNT(2), SXTEEN(2), THIRD(2), FOURTH(2),
     *   FIFTH(2), SIXTH(2), SEVNTH(2), EIGHTH(2), NINTH(2), TENTH(2),
     *   ELEVTH(2), TWLVTH(2), THRTTH(2), FORTTH(2), TT7TH(2), TINY2,
     *   CBTEP, SQT3(2), ODSQT3(2)
C
      COMMON /MTHCNS/  MAXX, MXLGM1, NEGINF, POSINF, A, PI, E, PI2, PI3,
     *   PI4, PI6, PI8, ONE, OD2F, OD3F, E14, A3, OD4F, OD5F, OD6F,
     *   OD7F, OD8F, OD9F, OD10F, OD11F, OD12F, OD14F, EIGHT, ZERO, TWO,
     *   THREE, FOUR, SXTNTH, NINE, TEN, TWOT7, SQT10, ESXTNT, SXTEEN,
     *   THIRD, FOURTH, FIFTH, SIXTH, SEVNTH, EIGHTH, NINTH, TENTH,
     *   ELEVTH, TWLVTH, THRTTH, FORTTH, TT7TH, TINY2, CBTEP, SQT3,
     *   ODSQT3
C
      INTEGER ITINY2, JTINY2
C
      COMMON /IMATH/   ITINY2, JTINY2
C
C
C  The above common blocks hold mathematical constants which are used in
C  the elementary function routines.
C
C  Variable descriptions
C
C  MAXX        a double precision representation of the largest
C              representable integer
C
C  MXLGM1      an approximation to the logarithm of .25 times the
C              largest representable machine number, minus 1. This
C              should be a rigorous lower bound on the logarithm of the
C              largest representable machine number.  Its default
C              computation in SIMINI is to use the Fortran LOG function.
C              This should be changed if the Fortran LOG function is not
C              sufficiently accurate.
C
C
C  A           an interval enclosure for 1/pi
C
C  PI          an interval enclosure for pi
C
C  PI2         an interval enclosure for pi/2
C
C  PI3         an interval enclosure for pi/3
C
C  PI4         an interval enclosure for pi/4
C
C  PI6         an interval enclosure for pi/6
C
C  PI8         an interval enclosure for pi/8
C
C  E           an interval enclosure for E
C
C  E14         an interval enclosure for e^{1/4}
C
C  ESXTNT      an interval enclosure for e^{1/16}
C
C  SQT10       an interval enclosure for SQRT(10)
C
C  TINY2       the maximum of the smallest machine number and the
C              reciprocal of the largest machine number.  This quantity
C              is checked to avoid overflow in certain places.
C
C  CBTEP       an approximation to the cube root of six times the
C              largest distance between numbers.
C
C  SQT3        an interval enclosure for the square root of 3.
C
C  ODSQT3      the reciprocal of the square root of 3 times
C              1+ 100*MXULP, used in argument reduction in the
C              arctangent routine.
C
C  ITINY2      the logarithm base 10 of TINY2 truncated to an integer.
C
C  JTINY2      sixteen times the logarithm base e of TINY2 truncated to
C              an integer.
C
C  See the statements in the routine SIMINI for the definitions of
C  the other constants.
C
      INTEGER ISIG, IERR, IROUT, IPRTCL, ISEVER, IERPUN
      COMMON /ERFLGS/ ISIG, IERR, IROUT, IPRTCL, ISEVER, IERPUN
C
C The above common block stores signalling information for error
C conditions.  In particular,
C
C ISIG   is set to 0 at the beginning of SIMINI, and is reset
C        to the severity level of the error (1, 2, or 3) when an error
C        condition occurs.  The user may reset the flag to zero after an
C        error condition, depending on the error, if ISEVER is set to
C        allow execution after an error of its severity.
C
C IERR   is the number of the error condition, if an error occurred in
C        the last routine in INTLIB with error checking.  If no errors
C        occurred in the last such routine called, then IERR is zero.
C        The specific error conditions associated with particular error
C        numbers is defined in the routine ERRTST.
C
C IROUT  is the code number of the package routine in which the error
C        occurred
C
C IPRTCL controls the level of error which is printed.  IPRCTL=0 prints
C        all levels, while IPRCTL=1 prints only errors of level 1 or
C        greater.  If IPRTCL=4, then no error information is printed.
C
C ISEVER gives the error level which will stop execution.  ISEVER=0
C        causes any error to stop execution, while ISEVER>=3 causes only
C        errors of severity 3 to stop execution.  Generally, severity 3
C        errors correspond to when it is impossible to assign a result
C        with any meaningful interpretation.
C
C IERPUN is the Fortran unit number to which errors should be printed.
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  ERRTST, MULT, RNDOUT
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Identifying code for this routine --
C
      IROUT = 3
      IERR = 0
C
C  Do usual interval division if zero is not in the denominator --
C
      IF (BB(1).GT.ZERO(2)) THEN
         C(1) = ONE(1)/BB(2)
         C(2) = ONE(2)/BB(1)
         CALL RNDOUT(C,.TRUE.,.TRUE.)
         CALL MULT(AA,C,RESULT)
      ELSE IF (BB(2).LT.ZERO(1)) THEN
         C(1) = ONE(2)/BB(2)
         C(2) = ONE(1)/BB(1)
         CALL RNDOUT(C,.TRUE.,.TRUE.)
         CALL MULT(AA,C,RESULT)
      ELSE
         IERR = 6
         CALL ERRTST(BB)
         RESULT(1) = NEGINF
         RESULT(2) = POSINF
      END IF
C
      RETURN
C
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE MULT (A, B, RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
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   November 17, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  MULT multiplies the interval A and the interval B.  It
C  puts the result into output parameter RESULT.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), B(2), RESULT(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TEMP, A1, B1
      DOUBLE PRECISION AA(2), BB(2)
C
C***********************************************************************
C
C Internal Constant Declarations --
C
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO = 0.0D0)
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Pictures for cases.
C
      IF ((ZERO .LE. A(1)) .AND. (ZERO .LE. B(1))) THEN
C  Case 1   ---------------- 0 -----------------
C        A:                    [==========]
C        B:                         [===========]
          RESULT(1) = A(1) * B(1)
          RESULT(2) = A(2) * B(2)
      ELSE IF ((A(1) .LT. ZERO) .AND. (ZERO .LT. A(2))
     1                          .AND. (ZERO .LE. B(1))) THEN
C  Case 2   ---------------- 0 -----------------
C        A:            [==================]
C        B:                         [===========]
          RESULT(1) = A(1) * B(2)
          RESULT(2) = A(2) * B(2)
      ELSE IF ((A(2) .LE. ZERO) .AND. (ZERO .LE. B(1))) THEN
C  Case 3   ---------------- 0 -----------------
C        A:    [==========]
C        B:                         [===========]
          B1 = B(1)
          RESULT(1) = A(1) * B(2)
          RESULT(2) = A(2) * B1
      ELSE IF ((ZERO .LE. A(1)) .AND. (B(1) .LT. ZERO)
     1                          .AND. (ZERO .LT. B(2))) THEN
C  Case 4   ---------------- 0 -----------------
C        A:                     [==========]
C        B:               [===========]
          RESULT(1) = A(2) * B(1)
          RESULT(2) = A(2) * B(2)
      ELSE IF ((A(2) .LE. ZERO) .AND. (B(1) .LT. ZERO)
     1                          .AND. (ZERO .LT. B(2))) THEN
C  Case 5   ---------------- 0 -----------------
C        A:    [==========]
C        B           [===========]
          A1 = A(1)
          B1 = B(1)
          RESULT(1) = A(1) * B(2)
          RESULT(2) = A1 * B1
      ELSE IF ((ZERO .LE. A(1)) .AND. (B(2) .LE. ZERO)) THEN
C  Case 6   ---------------- 0 -----------------
C        A:                    [==========]
C        B:  [===========]
          A1 = A(1)
          RESULT(1) = A(2) * B(1)
          RESULT(2) = A1 * B(2)
      ELSE IF ((A(1) .LT. ZERO) .AND. (ZERO .LT. A(2))
     1                          .AND. (B(2) .LE. ZERO)) THEN
C  Case 7   ---------------- 0 -----------------
C        A:              [==========]
C        B:  [===========]
          A1 = A(1)
          B1 = B(1)
          RESULT(1) = A(2) * B(1)
          RESULT(2) = A1 * B1
      ELSE IF ((A(2) .LE. ZERO) .AND. (B(2) .LE. ZERO)) THEN
C  Case 8   ---------------- 0 -----------------
C        A:    [==========]
C        B:   [===========]
          A1 = A(1)
          B1 = B(1)
          RESULT(1) = A(2) * B(2)
          RESULT(2) = A1 * B1
      ELSE
C  Case 9   ---------------- 0 -----------------
C        A:              [==========]
C        B:              [===========]
C  Must check two cases.
          AA(1) = A(1)
          AA(2) = A(2)
          BB(1) = B(1)
          BB(2) = B(2)
          RESULT(1) = AA(1) * BB(2)
          TEMP = AA(2) * BB(1)
          IF (TEMP .LT. RESULT(1)) THEN
              RESULT(1) = TEMP
          ELSE
          END IF
              RESULT(2) = AA(1) * BB(1)
          TEMP = AA(2) * BB(2)
          IF (TEMP .GT. RESULT(2)) THEN
              RESULT(2) = TEMP
          ELSE
          END IF
      END IF
C
      CALL RNDOUT(RESULT,.TRUE.,.TRUE.)
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   Modified August, 1991 and April, 1992 by Kaisheng Du and
C   R. Baker Kearfott to be more accurate near the underflow
C   threshold.
C
C Part of the generalized bisection package
C (interval arithmetic subpackage)
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 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
      DOUBLE PRECISION TINY, TEST
      COMMON /MACH2/ TINY, TEST
C
C  See SIMINI for an explanation of the common block MACH2.
C
C***********************************************************************
C
C Beginning of executable statements --
C
      IF (RNDDWN) THEN
         IF (X(1).GE.TEST) THEN
            X(1) = (1.D0 - MXULP ) * X(1)
         ELSE IF (X(1).LE.-TEST) THEN
            X(1) = (1D0 + MXULP ) * X(1)
         ELSE IF (X(1).LE.0.D0) THEN
            X(1) = -TEST
         ELSE
            X(1) = 0.D0
         END IF
      END IF
C
      IF (RNDUP) THEN
         IF (X(2).GE.TEST) THEN
            X(2) = (1.D0 + MXULP )* X(2)
         ELSE IF (X(2).LE.-TEST) THEN
            X(2) = (1.D0 - MXULP ) * X(2)
         ELSE IF(X(2).GE.0D0) THEN
            X(2) = TEST
         ELSE
            X(2) = 0.D0
         ENDIF
      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 Package-supplied functions and subroutines --
C
C  RNDOUT
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 Modified by Manuel Novoa III on March 13, 1990 to clean the code
C   slightly, and to remove the need for MAX and MIN.
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 --  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
      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 IF (R.GT.0D0) THEN
         RESULT(1) = R * T1
         RESULT(2) = R * T2
      ELSE
         RESULT(1) = R * T2
         RESULT(2) = R * T1
      END IF
C
      CALL RNDOUT(RESULT,RNDDWN,RNDUP)
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
