c*** elemfuns.f
      SUBROUTINE IACOS (X,IRCCOS)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C   By Chenyi Hu
C       and
C   Abdulhamid Awad
C
C   Computer and Mathematical Sciences Department
C   University of Houston-Downtown
C   Houston, TX 77002
C
C       and
C
C   R. Baker Kearfott
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   July 12, 1992
C
C Part of the interval elementary function library
C
C
C***********************************************************************
C
C Function --
C
C  This routine computes an interval enclosure for the arccos over
C  the interval X and returns the result in the interval IRCCOS.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(2), IRCCOS(2)
C
C***********************************************************************
C
C Common blocks --
C
        DOUBLE PRECISION MXULP, TTINY2, TOL0
        COMMON /MACH1/MXULP, TTINY2, TOL0
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, IASIN, SUB
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Identifying code for this routine --
C
      IROUT = 8
      IERR = 0
C
      IF (X(2).LT.X(1)) THEN
         IERR = 1
         CALL ERRTST(X)
      END IF
C
      IF (X(2).GT.ONE(1)) THEN
         IERR = 9
         CALL ERRTST(X)
         IRCCOS(1) = NEGINF
         IRCCOS(2) = POSINF
         RETURN
      ELSE IF (X(1).LT.-ONE(1)+MXULP) THEN
         IERR = 10
         CALL ERRTST(X)
         IRCCOS(1) = NEGINF
         IRCCOS(2) = POSINF
         RETURN
      END IF
C
      CALL IASIN(X, IRCCOS)
C
      CALL SUB(PI2,IRCCOS,IRCCOS)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE IACOT(X,IRCCOT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   By Chenyi Hu
C       and
C   Abdulhamid Awad
C
C   Computer and Mathematical Sciences Department
C   University of Houston-Downtown
C   Houston, TX 77002
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   June 16, 1992
C
C Part of the interval elementary function library
C
C***********************************************************************
C
C Function --
C
C   This program finds bounds on the value of ARCCOT(X) as X
C   ranges of the interval XX, and places the resulting interval
C   in IRCCOT.
C
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(2), IRCCOT(2)
C
C***********************************************************************
C
C Common blocks--
C
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
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  IATAN, SUB
C
C***********************************************************************
C
C Beginning of executable statements --
C

      CALL IATAN(X, IRCCOT)
      CALL SUB(PI2,IRCCOT,IRCCOT)

      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE IASIN(XX,IRCSIN)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C   By Chenyi Hu
C       and
C   Abdulhamid Awad
C
C   Computer and Mathematical Sciences Department
C   University of Houston-Downtown
C   Houston, TX 77002
C
C       and
C
C   R. Baker Kearfott
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   July 7, 1992
C
C Part of the interval elementary function library
C
C
C***********************************************************************
C
C Function --
C
C  This routine computes an interval enclosure for the arcsin over
C  the interval XX and returns the result in the interval IRCSIN.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION XX(2), IRCSIN(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION LB, UB
      DOUBLE PRECISION X(2), RSLT(2), TMP(2)
      LOGICAL OVER, NEGATV, FLIP
C
C***********************************************************************
C
C Common blocks --
C
        DOUBLE PRECISION MXULP, TTINY2, TOL0
        COMMON /MACH1/MXULP, TTINY2, TOL0
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 Fortran-supplied functions and subroutines --  none
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  ASNSER, ERRTST, ISQRT, MULT, RNDOUT, SUB
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Identifying code for this routine --
C
      IROUT = 9
      IERR = 0
C
      IF (XX(2).LT.XX(1)) THEN
         IERR = 1
         CALL ERRTST(XX)
      END IF
C
      IF (XX(2).GT.ONE(1)) THEN
         IERR = 9
         CALL ERRTST(XX)
         RETURN
      ELSE IF (XX(1).LT.-ONE(1)+MXULP) THEN
         IERR = 10
         CALL ERRTST(XX)
         RETURN
      END IF
C
      LB = XX(1)
      UB = XX(2)
C
C For the lower end point ---
C----------------------------------------------------------------------
C
C  If LB < 0, convert it to positive --

      IF (LB .LT. ZERO(1)) THEN
         OVER = .TRUE.
         NEGATV = .TRUE.
         TMP(2) = -LB
         CALL RNDOUT(TMP,.FALSE.,.TRUE.)
         LB = TMP(2)
      ELSE
         OVER = .FALSE.
         NEGATV =.FALSE.
      ENDIF
C
C  Transform LB to a small interval X such that 0 <= X <= 0.5 --
C
      X(1) = LB
      X(2) = LB
      IF (NEGATV) CALL RNDOUT(X,.TRUE.,.TRUE.)
C
      IF (LB .LE. OD2F(2) ) THEN
         FLIP = .FALSE.
      ELSE
         FLIP = .TRUE.
C
C   X <-- SQRT( (1 - X)/2 )
C
         CALL SUB(ONE,X,X)
         CALL MULT(OD2F,X,X)
         IF(X(1).LT.0D0) THEN
            X(1)=0D0
         END IF
         CALL ISQRT(X,X)
      ENDIF
C
C  Find a bound on the arcsin at the lower end point --
C
        CALL ASNSER (X, OVER, RSLT)
C
C  Undo the argument transformation, if it was applied --
C
      IF (FLIP) THEN
         CALL MULT(TWO,RSLT,RSLT)
         CALL SUB(PI2,RSLT,RSLT)
      END IF

      IF (NEGATV) THEN
         IRCSIN(1) = -RSLT(2)
         CALL RNDOUT(IRCSIN,.TRUE.,.FALSE.)
      ELSE
         IRCSIN(1) = RSLT(1)
      ENDIF
C
C For the upper end point ---
C----------------------------------------------------------------------
C
      IF (UB .LT. 0D0) THEN
         OVER =.FALSE.
         NEGATV = .TRUE.
         TMP(2) = -UB
         CALL RNDOUT(TMP,.FALSE.,.TRUE.)
         UB = TMP(2)
      ELSE
         OVER = .TRUE.
         NEGATV = .FALSE.
      ENDIF
C
      X(1) = UB
      X(2) = UB
      IF (NEGATV) CALL RNDOUT(X,.TRUE.,.TRUE.)
C
      IF (UB .LE. OD2F(2) ) THEN
         FLIP = .FALSE.
      ELSE
         FLIP = .TRUE.
C
C   X <-- SQRT( (1 - X)/2 )
C
         CALL SUB(ONE,X,X)
         CALL MULT(OD2F,X,X)
         IF(X(1).LT.0D0) THEN
            X(1)=0D0
         END IF
         CALL ISQRT(X,X)
      ENDIF
C
      CALL ASNSER(X, OVER, RSLT)
C
      IF (FLIP) THEN
         CALL MULT(TWO,RSLT,RSLT)
         CALL SUB(PI2,RSLT,RSLT)
      END IF

      IF (NEGATV) THEN
         IRCSIN(2) = -RSLT(1)
         CALL RNDOUT(IRCSIN,.FALSE.,.TRUE.)
      ELSE
         IRCSIN(2) = RSLT(2)
      ENDIF
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
C
      SUBROUTINE ASNSER (X, OVER, RSLT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C   By Chenyi Hu
C       and
C   Abdulhamid Awad
C
C   Computer and Mathematical Sciences Department
C   University of Houston-Downtown
C   Houston, TX 77002
C
C       and
C
C   R. Baker Kearfott
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   July 7, 1992
C
C Part of the interval elementary function library
C
C
C***********************************************************************
C
C Function --
C
C  This routine computes an interval enclosure for the arcsin over
C  the interval X and returns the result in the interval RSLT.  It is
C  assumed that X is nearly a point, between 0 and .5.  This routine is
C  normally called from IASIN.  The argument OVER indicates whether
C  the upper end point or lower end point is important.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(2)
      LOGICAL OVER
      DOUBLE PRECISION RSLT(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TEMP1(2), TEMP2(2), TEMP3(2), SUM(2)
      DOUBLE PRECISION T
      INTEGER I
      DOUBLE PRECISION D2IM1(2), D2I(2), D2IP1(2)
C
C***********************************************************************
C
C Common blocks --
C
      DOUBLE PRECISION MXULP, TTINY2, TOL0
      COMMON /MACH1/MXULP, TTINY2, TOL0
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 Fortran-supplied functions and subroutines --
C
C  ABS, DBLE, MAX
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  ADD, ERRTST, IDIV, MULT, RNDOUT
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Identifying code for this routine --
C
      IROUT = 10
      IERR = 0
C
      SUM(1) = X(1)
      SUM(2) = X(2)
C
      TEMP1(1) = X(1) * X(1)
      TEMP1(2) = X(2) * X(2)
      CALL RNDOUT(TEMP1, .TRUE., .TRUE.)
C
      TEMP2(1) = X(1)
      TEMP2(2) = X(2)
C
      DO 10 I = 1,100
C
         D2IM1(1) = DBLE(2*I-1)
         D2IM1(2) = D2IM1(1)
         CALL RNDOUT(D2IM1,.TRUE.,.TRUE.)
         D2I(1) = DBLE(2*I)
         D2I(2) = D2I(1)
         CALL RNDOUT(D2I,.TRUE.,.TRUE.)
         D2IP1(1) = DBLE(2*I+1)
         D2IP1(2) = D2IP1(1)
         CALL RNDOUT(D2IP1,.TRUE.,.TRUE.)
C
C  TEMP2 <-- (2i-1)TEMP1 TEMP2 / (2i)
C
         CALL MULT(D2IM1,TEMP2,TEMP2)
         CALL MULT(TEMP2,TEMP1,TEMP2)
         CALL IDIV(TEMP2,D2I,TEMP2)
C
         CALL IDIV(TEMP2,D2IP1,TEMP3)
C
         SUM(1) = SUM(1) + TEMP3(1)
         SUM(2) = SUM(2) + TEMP3(2)
         CALL RNDOUT(SUM,.TRUE.,.TRUE.)
         T = MAX( ABS(SUM(1)), TOL0 )
         IF (TEMP3(2)/T .LT. MXULP) GOTO 20
   10 CONTINUE
         IERR=11
         CALL ERRTST(TEMP3)
         RETURN
   20 CONTINUE
C
      IF (OVER) THEN
         CALL MULT(TWO,TEMP3,TEMP3)
         CALL ADD(SUM,TEMP3,RSLT)
      ELSE
         RSLT(1) = SUM(1)
         RSLT(2) = SUM(2)
      ENDIF
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE IATAN(XX,IRCTAN)

C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   By Chenyi Hu
C       and
C   Abdulhamid Awad
C
C   Computer and Mathematical Sciences Department
C   University of Houston-Downtown
C   Houston, TX 77002
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   June 16, 1992
C
C Part of the interval elementary function library
C
C***********************************************************************
C
C Function --
C
C   This program finds bounds on the value of ARCTAN(X) as X
C   ranges over the interval XX, and places the resulting interval
C   in IRCTAN.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION XX(2), IRCTAN(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION LB, UB, X(2), RSLT(2), TMP(2)
      INTEGER CODE
      LOGICAL EVEN, NEGATV, RLB, RUB
C
C***********************************************************************
C
C Common blocks--
C
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
      DOUBLE PRECISION MXULP, TTINY2, TOL0
      COMMON /MACH1/MXULP, TTINY2, TOL0
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  ADD, ATNRED, ATNSER, RNDOUT, SUB
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Identifying code for this routine --
C
      IROUT = 11
      IERR = 0
C
      IF (XX(2).LT.XX(1)) THEN
         IERR = 1
         CALL ERRTST(XX)
      END IF
C
      LB = XX(1)
      UB = XX(2)
C
C  Find the logical value EVEN, used in subroutine ATNSER to determine
C  the number of terms to be computed in the series --
C
      IF (LB .GT. 1.0D0) THEN
         EVEN = .TRUE.
      ELSE IF (LB .GE. 0D0) THEN
         EVEN = .FALSE.
      ELSE IF (LB .GT. -1D0) THEN
         EVEN = .FALSE.
      ELSE
         EVEN = .TRUE.
      ENDIF
C
C  If LB < 0, convert it to positive and update even --

      IF (LB .LT. 0D0) THEN
         EVEN = .NOT. EVEN
         NEGATV = .TRUE.
         TMP(2) = -LB
         CALL RNDOUT(TMP,.FALSE.,.TRUE.)
         LB = TMP(2)
      ELSE
         NEGATV =.FALSE.
      ENDIF
C
C  Transform LB to an interval X such that 0 <= X < 1/sqrt(3).  Here, X
C  is a small interval containing the lower bound of the the interval
C  [LB,UB] after transformation --
C
      CALL ATNRED(LB, CODE, X)
C
C  Find the value of the arctangent for the interval X --
C
      CALL ATNSER(X, EVEN, RSLT)
C
C  Undo the transformation and store the result --
C
      IF (CODE .EQ. 11) THEN
      ELSE IF (CODE .EQ. 12 ) THEN
         CALL ADD(PI6,RSLT,RSLT)
      ELSE IF (CODE .EQ. 21 ) THEN
         CALL SUB(PI2,RSLT,RSLT)
      ELSE IF (CODE .EQ. 22 ) THEN
         CALL SUB(PI3,RSLT,RSLT)
      END IF
C
      IF (NEGATV) THEN
         IRCTAN(1) = -RSLT(2)
         RLB = .TRUE.
      ELSE
         IRCTAN(1) = RSLT(1)
         RLB = .FALSE.
      ENDIF
C
C  Similarly obtain an enclosure for the arctangent at the upper
C  end point --
C
      IF (UB .GT. 1.0D0) THEN
         EVEN = .FALSE.
      ELSE IF (UB .GE. 0D0) THEN
         EVEN = .TRUE.
      ELSE IF (UB .GT. -1D0) THEN
         EVEN = .TRUE.
      ELSE
         EVEN = .FALSE.
      ENDIF
C
      IF (UB .LT. 0D0) THEN
         EVEN = .NOT. EVEN
         NEGATV = .TRUE.
         TMP(2) = -UB
         CALL RNDOUT(TMP,.FALSE.,.TRUE.)
         UB = TMP(2)
      ELSE
         NEGATV = .FALSE.
      ENDIF
C
      CALL ATNRED(UB, CODE, X)
      CALL ATNSER(X, EVEN, RSLT)
C
      IF (CODE .EQ. 11) THEN
      ELSE IF (CODE .EQ. 12 ) THEN
         CALL ADD(PI6,RSLT,RSLT)
      ELSE IF (CODE .EQ. 21 ) THEN
         CALL SUB(PI2,RSLT,RSLT)
      ELSE IF (CODE .EQ. 22 ) THEN
         CALL SUB(PI3,RSLT,RSLT)
      END IF
C
      IF (NEGATV) THEN
         IRCTAN(2) = -RSLT(1)
         RUB = .TRUE.
      ELSE
         IRCTAN(2) = RSLT(2)
         RUB = .FALSE.
      ENDIF
C
      CALL RNDOUT(IRCTAN,RLB,RUB)
C
      RETURN
      END
C**********************************************************************
C**********************************************************************
      SUBROUTINE ATNRED(PP,CODE,X)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   By Chenyi Hu
C       and
C   Abdulhamid Awad
C
C   Computer and Mathematical Sciences Department
C   University of Houston-Downtown
C   Houston, TX 77002
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   June 16, 1992
C
C Part of the interval elementary function library
C
C***********************************************************************
C
C Function --
C
C  This is an auxiliary routine for IATAN.  It performs an argument
C  reduction in preparation for computing a series approximation to the
C  arctangent function.
C
C  The subroutine transforms a given nonnegative number P to a
C  nonnegative number less than the reciprocal of the square root of 3,
C  then converts it to an interval.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION PP, X(2)
      INTEGER CODE
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 PP     A floating point representation of an endpoint, to be reduced.
C        (INPUT)
C
C X      A small interval enclosure for the reduced quantity
C        corresponding to PP.
C        (OUTPUT)
C
C CODE   A flag indicating how the argument was reduced, as follows:
C
C          CODE = 11      no transformation
C          CODE = 12      1 / sqrt(3) < P < = 1
C          CODE = 21      P > 1 and 0 < 1/P < 1 / sqrt(3)
C          CODE = 22      P > 1 and 1/P >= 1 / sqrt(3)
C
C        (OUTPUT)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION P(2), TMP(2)
C
C***********************************************************************
C
C Common blocks --
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
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  ADD, IDIV, MULT, SUB
C
C***********************************************************************
C
C Beginning of executable statements --
C
      P(1) = PP
      P(2) = PP
C
      IF (P(1) .LE. ODSQT3(1)) THEN
         CODE = 11
      ELSE IF (P(1) .LE. ONE(2)) THEN
C
C  P <-- (sqrt(3)P - 1)/(sqrt(3)+P):
C
         CALL MULT(P,SQT3,TMP)
         CALL SUB(TMP,ONE,TMP)
         CALL ADD(SQT3,P,P)
         CALL IDIV(TMP,P,P)
         CODE = 12
      ELSE
C
C  P <-- 1/P:
C
         CALL IDIV(ONE,P,P)
         IF (P(1) .LE. ODSQT3(1)) THEN
            CODE = 21
         ELSE
            CALL MULT(P,SQT3,TMP)
            CALL SUB(TMP,ONE,TMP)
            CALL ADD(SQT3,P,P)
            CALL IDIV(TMP,P,P)
            CODE = 22
         ENDIF
      ENDIF
C
      X(1) = P(1)
      X(2) = P(2)
C
      RETURN
      END
C**********************************************************************
C**********************************************************************
      SUBROUTINE ATNSER (X,EVEN,RSLT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   By Chenyi Hu
C       and
C   Abdulhamid Awad
C
C   Computer and Mathematical Sciences Department
C   University of Houston-Downtown
C   Houston, TX 77002
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   June 16, 1992
C
C Part of the interval elementary function library
C
C***********************************************************************
C
C Function --
C
C  This subroutine computes an enclosure for arctan X and puts the
C  result in RSLT.  It is assumed that X is an interval of small width
C  between 0 and 1 / sqrt(3).  The argument EVEN is used to determine
C  whether an even or odd number of terms should be taken.  This routine
C  is normally called from IATAN.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(2), RSLT(2)
      LOGICAL EVEN
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TEMP1(2), TEMP2(2), TEMP3(2), SUM(2), XS(2)
      DOUBLE PRECISION DI(2), T
      INTEGER I, SIGN
C
C***********************************************************************
C
C Common blocks --
C
      DOUBLE PRECISION MXULP, TTINY2, TOL0
      COMMON /MACH1/ MXULP, TTINY2, TOL0
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 Fortran-supplied functions and subroutines --
C
C  DBLE
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  ADD, IDIV, MULT, POWER, RNDOUT
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Identifying code for this routine --
C
      IROUT = 12
      IERR = 0
C
      IF (X(1).LT.0D0) THEN
         IERR=12
         CALL ERRTST(X)
         RETURN
      END IF
C
      CALL POWER(X,2,XS)
C
      IF ((X(1) .EQ. 0D0) .AND. (X(2) .EQ. 0D0)) THEN
         RSLT(1) = ZERO(1)
         RSLT(2) = ZERO(2)
         RETURN
      ENDIF
C
      SUM(1) = X(1)
      SUM(2) = X(2)
C
      SIGN = -1
      CALL POWER(X,3,TEMP1)
C
      DO 10 I = 1,10000
         SIGN = -SIGN
         CALL MULT(TWO,TEMP1,TEMP2)
         DI(1) = DBLE (4*I*I - 1)
         DI(2) = DI(1)
         CALL RNDOUT(DI,.TRUE.,.TRUE.)
         CALL IDIV(TEMP2,DI,TEMP2)
C
         IF (SIGN .GT. 0) THEN
            SUM(1) = SUM(1) + TEMP2(1)
            SUM(2) = SUM(2) + TEMP2(2)
         ELSE
            SUM(1) = SUM(1) - TEMP2(2)
            SUM(2) = SUM(2) - TEMP2(1)
         ENDIF
         CALL RNDOUT(SUM,.TRUE.,.TRUE.)
C
         TEMP1(1) = TEMP1(1) * XS(1)
         TEMP1(2) = TEMP1(2) * XS(2)
         CALL RNDOUT(TEMP1, .TRUE., .TRUE.)
C
         T = MAX( ABS(SUM(1)),2D0*TOL0 )
         IF (TEMP1(2)/T .GT. MXULP) GOTO 10
         IF (EVEN .AND. (SIGN .LT. 0)) GOTO 10
         IF (EVEN .OR. (SIGN .LE. 0)) GOTO 20
   10 CONTINUE
         IERR=11
         CALL ERRTST(TEMP2)
         RETURN
   20 CONTINUE
C
      CALL ADD(XS,ONE,TEMP3)
      CALL IDIV(SUM,TEMP3,RSLT)
C
        RETURN
        END
C***********************************************************************
C***********************************************************************
      SUBROUTINE ICOS(X,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Kaisheng Du
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C
C
C***********************************************************************
C
C Function --
C
C  This routine returns the interval value of the cosine function
C  (evaluated over the interval X) in the interval variable
C  RESULT.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(2), RESULT(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION V1(2), V2(2)
      DOUBLE PRECISION VAL1(2),VAL2(2)
      DOUBLE PRECISION TMP(2), T1, T2
C
      INTEGER N1, N2
C
      LOGICAL EVEN
C
C***********************************************************************
C
C Common blocks --
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 Fortran-supplied functions and subroutines --
C
C     ABS, ERRTST, IDINT, MAX, MIN, MOD
C
C***********************************************************************
C
C Package-supplied functions and subroutines -- none
C
C     EXTERNAL ERRTST, MULT, RCOS, RNDOUT
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Identifying code for this routine --
C
      IROUT = 2
      IERR = 0
C
      IF (X(2).LT.X(1)) THEN
         IERR = 1
         CALL ERRTST(X)
      END IF
C
C  Compute the value at the left endpoint using interval arithmetic --
C
      V1(1) = X(1)
      V1(2) = X(1)
C
C  If X1/PI is out of range, return the interval [-1,1] --
C
      CALL MULT(V1,A,TMP)
      T1 =ABS(TMP(1))
      IF(T1.GE.MAXX) THEN
         IERR = 5
         CALL ERRTST(X)
         RESULT(1) = -1D0
         RESULT(2) =  1D0
         RETURN
      ENDIF
C
C  Compute the value at the right endpoint using interval arithmetic --
C
      V2(1) = X(2)
      V2(2) = X(2)
      CALL RNDOUT(V2,.TRUE.,.TRUE.)
C
C  If X2/pi is out of range, return the interval [-1,1] --
C
      CALL MULT(V2,A,TMP)
      T2 = ABS(TMP(2))
      IF(T2.GE.MAXX) THEN
         IERR = 5
         CALL ERRTST(X)
         RESULT(1) = -1D0
         RESULT(2) =  1D0
         RETURN
      ENDIF
C
C  If both X1/pi and X2/pi are within range, calculate the function
C  values at the left and right endpoints using interval arithmetic --
C
      CALL RCOS(V1,VAL1)
C
      CALL RCOS(V2,VAL2)
C
C  Calculate the number of half-periods from zero for the left
C  and right endpoints in order to normalize to [-pi/2,pi/2] --
C
      N1 = IDINT(T1)
      N2 = IDINT(T2)
      IF(V1(1).LT.0D0) N1 = -N1-1
      IF(V2(2).LT.0D0) N2 = -N2-1
C
C   In even half periods, the function is decreasing --
C
      EVEN = MOD(N1,2).EQ.0
C
C   If X1 and X2 are in the same half period, then the lower bound
C   and upper bound on the range occur at the endpoints --
C
      IF(N1.EQ.N2) THEN
         RESULT(1) = MIN(VAL1(1),VAL2(1))
         RESULT(2) = MAX(VAL1(2),VAL2(2))
         RETURN
      ENDIF
C
C   Consider the case X2 is in the half period adjacent to X1 --
C
      IF(N2.EQ.N1+1) THEN
C
C    If X1 is in the increasing half period, then the upper bound should
C    equal 1;  otherwise, the lower bound should be -1 --
C
         IF(EVEN) THEN
            RESULT(1) = -1.D0
            RESULT(2) = MAX(VAL1(2),VAL2(2))
         ELSE
           RESULT(1) = MIN(VAL1(1),VAL2(1))
           RESULT(2) = 1.D0
         ENDIF
C
C   If X1 and X2 are not in adjecent half periods, then the lower
C   bound should be -1 and the upper bound should be 1 --
C
      ELSE
         RESULT(1) = -1.D0
         RESULT(2) =  1.D0
      ENDIF
C
      CALL RNDOUT(RESULT,.TRUE.,.TRUE.)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE IEXP(X,RESULT)
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Kaisheng Du
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C
C***********************************************************************
C
C Function --
C
C  This routine returns the interval exp(X) in RESULT.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(2), RESULT(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TEMP(2), TVAL(2), T
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 Fortran-supplied functions --
C
C     ABS
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
      EXTERNAL ERRTST, REXP
      DOUBLE PRECISION D1MACH
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Identifying code for this routine --
C
      IROUT = 1
      IERR = 0
C
      IF (X(2).LT.X(1)) THEN
         IERR = 1
         CALL ERRTST(X)
      END IF
C
      IF (X(1).EQ.ZERO(1)) THEN
         RESULT(1) =ONE(1)
      ELSE
         T = ABS(X(1))
         IF (T.GT.MXLGM1) THEN
            IF(X(1).GT.ZERO(1))THEN
               IERR = 3
            ELSE
               IERR = 2
            END IF
            CALL ERRTST(X)
         END IF
         IF (IERR.EQ.2) THEN
            RESULT(1) = ZERO(1)
         ELSE IF (IERR.EQ.3) THEN
            RESULT(1) = POSINF
            RESULT(2) = POSINF
            RETURN
         ELSE
            TEMP(1) = X(1)
            TEMP(2) = X(1)
            CALL REXP(TEMP,TVAL)
            RESULT(1) = TVAL(1)
         END IF
      END IF
C
      IF (X(2).EQ.ZERO(2)) THEN
         RESULT(2) =ONE(2)
      ELSE
         IF (X(2).GT.MXLGM1) THEN
            IERR = 4
            CALL ERRTST(X)
         END IF
         IF (IERR.EQ.4) THEN
            RESULT(2) = D1MACH(2)
            RETURN
         ELSE
            TEMP(1) = X(2)
            TEMP(2) = X(2)
            CALL REXP(TEMP,TVAL)
            RESULT(2) = TVAL(2)
         END IF
      END IF
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE IIPOWR(XX, YY, RESULT)
C
C  Kaisheng Du and R. B. Kearfott, Summer, 1991
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Kaisheng Du
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C
C**********************************************************************
C
C Function --
C
C  This routine returns the interval X**Y, provided that X is a
C  positive interval.
C  RESULT.
C
C**********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION XX(2), YY(2), RESULT(2)
C
C**********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION X(2), Y(2)
      LOGICAL LBZERO

C**********************************************************************
C
C Common blocks --
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
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 Fortran-supplied functions and subroutines --
C
C   ABS, MAX
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C   ERRTST, IEXP, ILOG, MULT
C
C**********************************************************************
C
C Beginning of executable statements --
C
C  Identifying code for this routine --
C
      IROUT = 14
      IERR = 0
C
      X(1) = XX(1)
      X(2) = XX(2)
      Y(1) = YY(1)
      Y(2) = YY(2)
C
      IF (X(2).LT.X(1)) THEN
         IERR = 1
         CALL ERRTST(X)
      END IF
C
      LBZERO = .FALSE.
      IF (X(1).EQ.0D0) THEN
         IF (Y(1).LT.0D0) THEN
            IF(Y(2).GE.0D0) THEN
               IERR = 14
               RESULT(1) = NEGINF
               RESULT(2) = POSINF
               RETURN
            END IF
         END IF
         LBZERO = .TRUE.
         X(1) = X(2)
      END IF

      IF (X(1).LE.ZERO(2)) THEN
         IERR = 7
         CALL ERRTST(X)
         RETURN
      END IF
C
      CALL ILOG(X,RESULT)
C
      CALL MULT(RESULT,Y,RESULT)
C
      IF (ABS(RESULT(1)).GT.MXLGM1) THEN
         IF(RESULT(1).GT.ZERO(1)) THEN
            IERR = 3
         ELSE
            IERR = 2
         END IF
         CALL ERRTST(RESULT)
      END IF
      IF (RESULT(2).GT.MXLGM1) THEN
         IERR = 4
         CALL ERRTST(RESULT)
      END IF
C
      CALL IEXP(RESULT,RESULT)
      IF(LBZERO) THEN
         IF (Y(1).LT.0D0) THEN
            RESULT(2) = POSINF
            IERR = 15
         ELSE
            RESULT(1) = 0D0
         END IF
      END IF
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE ILOG(X,RESULT)
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Kaisheng Du
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C
C***********************************************************************
C
C Function --
C
C  This routine returns the interval log(X) in RESULT.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(2), RESULT(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TEMP(2), TVAL(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, RLOG
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Identifying code for this routine --
C
      IROUT = 4
      IERR = 0
C
      IF (X(2).LT.X(1)) THEN
         IERR = 1
         CALL ERRTST(X)
      END IF
C
      IF (X(1).LE.ZERO(2)) THEN
         IERR = 7
         CALL ERRTST(X)
         RETURN
      END IF
C
      TEMP(1) = X(1)
      TEMP(2) = X(1)
      CALL RLOG(TEMP,TVAL)
      RESULT(1) = TVAL(1)
C
      TEMP(1) = X(2)
      TEMP(2) = X(2)
      CALL RLOG(TEMP,TVAL)
      RESULT(2) = TVAL(2)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE ISIN(X,RESULT)
C
C  Kaisheng Du and R. B. Kearfott, Summer, 1991
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Kaisheng Du
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C
C**********************************************************************
C
C Function --
C
C  This routine returns the interval value of the sine function
C  (evaluated over the interval X) in the interval variable
C  RESULT.
C
C**********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(2), RESULT(2)
C
C**********************************************************************
C
C Common blocks --
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
C***********************************************************************
C
C Fortran supplied functions and subroutines --
C
C      ABS, MAX
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C      ICOS, MULT, POWER, SUB
C
C**********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TMP1(2), TMP2(2), TMP3(2)
      LOGICAL SMALL1, SMALL2
C
C**********************************************************************
C
C
C Beginning of executable statements --
C
C  Assure good relative accuracy near zero by the case separately;
C  otherwise, shift the argument and use the interval cosine function --
C
      SMALL1=.FALSE.
      SMALL2=.FALSE.
      IF ( MAX( ABS(X(1)), ABS(X(2)) ) .LT. PI2(1) ) THEN
         IF(ABS(X(1)).LT.CBTEP) THEN
            SMALL1 = .TRUE.
            TMP1(1) = X(1)
            TMP1(2) = X(1)
            IF (X(1).GT.ZERO(2)) THEN
               CALL POWER(TMP1,3,TMP2)
               CALL MULT(TMP2,OD3F,TMP2)
               CALL SUB(TMP1,TMP2,TMP1)
            END IF
         END IF
         IF(ABS(X(2)).LT.CBTEP) THEN
            SMALL2 = .TRUE.
            TMP3(1) = X(2)
            TMP3(2) = X(2)
            IF (X(2).LT.ZERO(1)) THEN
               CALL POWER(TMP3,3,TMP2)
               CALL MULT(TMP2,OD3F,TMP2)
               CALL SUB(TMP3,TMP2,TMP3)
            END IF
         END IF
      END IF
C
      IF (SMALL1 .AND. SMALL2) THEN
         RESULT(1) = TMP1(1)
         RESULT(2) = TMP3(2)
         RETURN
      ELSE
         CALL SUB(X,PI2,RESULT)
         CALL ICOS(RESULT,RESULT)
      END IF
C
      IF (SMALL1) RESULT(1) = MAX(TMP1(1),RESULT(1))
      IF (SMALL2) RESULT(2) = MIN(TMP3(2),RESULT(2))
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE ISINH(XX,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   Chenyi Hu
C
C   Computer and Mathematical Sciences Department
C   University of Houston-Downtown
C   Houston, TX 77002
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   October 20, 1992
C
C Part of the interval elementary function library
C
C***********************************************************************
C
C Function --
C
C  This routine returns an interval inclusion for the hyperbolic sine
C  over the interval XX in RESULT.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION XX(2), RESULT(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION AA, BB, X(2), RSLT(2)
      DOUBLE PRECISION TMP(2), TMP2(2)
      INTEGER CODE
      LOGICAL OVER, NEGATV
C
C***********************************************************************
C
C Common blocks --
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
      DOUBLE PRECISION MXULP, TTINY2, TOL0
      COMMON /MACH1/MXULP, TTINY2, TOL0
C
C***********************************************************************
C
C Package-supplied functions and subroutines --

C  ERRTST, IDIV, IEXP, ISHSER, ISNRED, ISNVAL, MULT, RNDOUT, SUB
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Identifying code for this routine --
C
      IROUT = 6
      IERR = 0
C
      IF (XX(2).LT.XX(1)) THEN
         IERR = 1
         CALL ERRTST(XX)
      END IF
C
      AA = XX(1)
      BB = XX(2)
C
C-------------------------------------------------------------------
C  Compute the value of the hyperbolic sine at the lower endpoint A.
C-------------------------------------------------------------------
C
      IF (AA .LT. ZERO(1)) THEN
         OVER = .TRUE.
         NEGATV = .TRUE.
         TMP(2) = -AA
         CALL RNDOUT(TMP,.FALSE.,.TRUE.)
         AA = TMP(2)
      ELSE
         OVER = .FALSE.
         NEGATV =.FALSE.
      ENDIF
C
C  Transform AA to an interval X such that 0 <= X <= 1 --
C
      CALL ISNRED(AA, CODE, X)
C
      IF (CODE .EQ. 20) THEN
C
C  Use the exponential function for arguments larger than 27 --
C
         CALL IEXP(X, TMP)
         CALL IDIV(ONE,TMP,TMP2)
         CALL SUB(TMP,TMP2,TMP)
         CALL MULT(OD2F,TMP,RSLT)
      ELSE
C
C  Evaluate the series X + X^3/3! + X^5/5! + ... for 0 <= X <= 1 --
C
         CALL ISHSER(X, OVER, RSLT)
         CALL ISNVAL(CODE, RSLT)
       END IF
C
C  Undo the possible argument reduction from ISNRED --
C
      IF (NEGATV) THEN
         RESULT(1) = - RSLT(2)
         CALL RNDOUT(RESULT,.TRUE.,.FALSE.)
      ELSE
         RESULT(1) = RSLT(1)
      ENDIF
C
C-------------------------------------------------------------------
C  Compute the value of the hyperbolic sine at the upper endpoint B
C  (completely analogous to computation for the lower endpoint A).
C-------------------------------------------------------------------
C
      IF (BB .LT. ZERO(1)) THEN
         OVER =.FALSE.
         NEGATV = .TRUE.
         TMP(2) = -BB
         CALL RNDOUT(TMP,.FALSE.,.TRUE.)
         BB = TMP(2)
      ELSE
         OVER = .TRUE.
         NEGATV = .FALSE.
      ENDIF
C
      CALL ISNRED(BB, CODE, X)
C
      IF (CODE .EQ. 20) THEN
         CALL IEXP(X, TMP)
         CALL IDIV(ONE,TMP,TMP2)
         CALL SUB(TMP,TMP2,TMP)
         CALL MULT(OD2F,TMP,RSLT)
      ELSE
         CALL ISHSER(X, OVER, RSLT)
         CALL ISNVAL(CODE, RSLT)
      END IF
C
      IF (NEGATV) THEN
         RESULT(2) = - RSLT(1)
         CALL RNDOUT(RESULT,.FALSE.,.TRUE.)
      ELSE
         RESULT(2) = RSLT(2)
      ENDIF
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE ISNRED(P,CODE,X)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   Chenyi Hu
C
C   Computer and Mathematical Sciences Department
C   University of Houston-Downtown
C   Houston, TX 77002
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   October 20, 1992
C
C Part of the interval elementary function library
C
C***********************************************************************
C
C Function --
C
C  This subroutine transforms the nonnegative number P to a nonnegative
C  number which less than 1, and then converts the result to an
C  interval.  The following code numbers are set to record how the
C  transformation was done:
C
C     CODE = 11      0 < P < =  1
C     CODE = 12      1 < P < =  3
C     CODE = 13      3 < P < =  9
C     CODE = 14      9 < P < = 27
C     CODE = 20          P >   27
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION P
      INTEGER CODE
      DOUBLE PRECISION X(2)
C
C***********************************************************************
C
C Common blocks --
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
C***********************************************************************
C
C Beginning of executable statements --
C
      X(1) = P
      X(2) = P
C
      IF (P .LE. ONE(1) ) THEN
         CODE = 11
         GO TO 10
      ELSE IF (P .LE. THREE(1)) THEN
         CALL MULT(X,THIRD,X)
         CODE = 12
         GO TO 10
      ELSE IF (P .LE. NINE(1)) THEN
         CALL MULT(X,NINTH,X)
         CODE = 13
         GO TO 10
      ELSE IF (P .LE. TWOT7(1)) THEN
         CALL MULT(X,TT7TH,X)
         CODE = 14
         GO TO 10
      ELSE
         CODE = 20
         GO TO 10
      ENDIF
C
   10 CALL RNDOUT(X, .TRUE., .TRUE.)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
        SUBROUTINE ISHSER (X, OVER, RSLT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   Chenyi Hu
C
C   Computer and Mathematical Sciences Department
C   University of Houston-Downtown
C   Houston, TX 77002
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   October 20, 1992
C
C Part of the interval elementary function library
C
C***********************************************************************
C
C Function --
C
C   This subroutine uses a series to bound the value of the hyperbolic
C   sine function over the nonnegative interval X of small width.  The
C   result is returned in RSLT.  The logical variable OVER indicates
C   whether negation was applied during the argument reduction.  This
C   routine is not meant to be called by the user.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(2)
      LOGICAL OVER
      DOUBLE PRECISION RSLT(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TEMP1(2), TEMP2(2), SUM(2), ERR(2), T
      INTEGER I
C
C***********************************************************************
C
C Common blocks --
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
      DOUBLE PRECISION MXULP, TTINY2, TOL0
      COMMON /MACH1/ MXULP, TTINY2, TOL0
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines --
C
C  ABS, DBLE, FLOAT, MAX
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  MULT, POWER, RNDOUT
C
C***********************************************************************
C
C Beginning of executable statements --
C
      SUM(1) = X(1)
      SUM(2) = X(2)

      I = 1
C
C  Use temp1 to hold X^2 --
C
      CALL POWER(X,2,TEMP1)

      TEMP2(1) = X(1)
      TEMP2(2) = X(2)
C
10    CONTINUE
      I = I + 2
         TEMP2(1) = TEMP1(1) * TEMP2(1)
         TEMP2(2) = TEMP1(2) * TEMP2(2)
         CALL RNDOUT(TEMP2, .TRUE., .TRUE.)
C
         TEMP2(1) = TEMP2(1)/DBLE(FLOAT(I))
         TEMP2(2) = TEMP2(2)/DBLE(FLOAT(I))
         CALL RNDOUT(TEMP2, .TRUE., .TRUE.)
C
         TEMP2(1) = TEMP2(1)/DBLE(FLOAT(I-1))
         TEMP2(2) = TEMP2(2)/DBLE(FLOAT(I-1))
         CALL RNDOUT(TEMP2, .TRUE., .TRUE.)
C
         SUM(1) = SUM(1) + TEMP2(1)
         SUM(2) = SUM(2) + TEMP2(2)
         CALL RNDOUT(SUM,.TRUE.,.TRUE.)
C
         CALL MULT(TWO,TEMP2,ERR)
         T = MAX(ABS(SUM(1)),ABS(SUM(2)),2D0*TOL0)
         IF ((ERR(1)/T) .GT. MXULP) GOTO 10
C
      IF (OVER) THEN
         RSLT(1) = SUM(1) + ERR(1)
         RSLT(2) = SUM(2) + ERR(2)
         CALL RNDOUT(RSLT, .TRUE., .TRUE.)
      ELSE
         RSLT(1) = SUM(1)
         RSLT(2) = SUM(2)
      ENDIF
C
      RETURN
      END

C **********************************************************
      SUBROUTINE ISNVAL (CODE,RSLT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Chenyi Hu
C
C   Computer and Mathematical Sciences Department
C   University of Houston-Downtown
C   Houston, TX 77002
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   October 20, 1992
C
C Part of the interval elementary function library
C
C***********************************************************************
C
C Function --
C
C   This routine transforms the result from the reduced argument
C   computation for the arcsine function back to the value corresponding
C   to the original argument.  The type of argument reduction which was
C   applied is given in the integer variable CODE.  Both the reduced
C   value input to this routine and the final value output from this
C   routine are stored in RSLT.
C
C***********************************************************************
C
C Argument declarations --
C
      INTEGER CODE
      DOUBLE PRECISION RSLT(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TEMP(2)
C
C***********************************************************************
C
C Common blocks --
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
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  ADD, MULT, RNDOUT
C
C
C***********************************************************************
C
C Beginning of executable statements --
C
      IF (CODE .EQ. 11) RETURN
C
   20 CONTINUE
      TEMP(1) = RSLT(1) * RSLT(1)
      TEMP(2) = RSLT(2) * RSLT(2)
      CALL RNDOUT(TEMP, .TRUE., .TRUE.)
C
C  rslt <-- rslt * (3 + 4 rslt^2)
C
      CALL MULT(FOUR,TEMP,TEMP)
      CALL ADD(THREE,TEMP,TEMP)
      CALL MULT(RSLT,TEMP,RSLT)
C
      IF (CODE .EQ. 12) THEN
         RETURN
      ELSE IF (CODE .EQ. 13) THEN
         CODE = 12
         GO TO 20
      ELSE IF (CODE .EQ. 14) THEN
         CODE = 13
         GO TO 20
      ENDIF
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE ISQRT(X,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Kaisheng Du
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C
C***********************************************************************
C
C Function --
C
C  This subroutine returns the interval value of the square root of X
C  in RESULT.  It uses the routine RSQRT to get rigorous values at the
C  endpoints.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(2), RESULT(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION V(2), RL(2), RR(2)
C
C***********************************************************************
C
C  Common block declarations --
C
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     ERRTST, RSQRT
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Identifying code for this routine --
C
      IROUT = 5
      IERR = 0
C
      IF (X(2).LT.X(1)) THEN
         IERR = 1
         CALL ERRTST(X)
      END IF
C
      IF (X(1).LT.0D0) THEN
         IERR = 7
         CALL ERRTST(X)
         RETURN
      ELSE IF (X(1).EQ.0D0) THEN
         RL(1) = 0D0
         RL(2) = 0D0
      ELSE
         V(1) = X(1)
         V(2) = X(1)
         CALL RSQRT(V,RL)
      END IF
C
      IF (X(2).EQ.0D0) THEN
         RR(1) = 0D0
         RR(2) = 0D0
      ELSE
         V(1) = X(2)
         V(2) = X(2)
         CALL RSQRT(V,RR)
      END IF
C
      RESULT(1) = RL(1)
      RESULT(2) = RR(2)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE POWER(AA,NDUM,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 fix a bug which
C   caused overestimation for N odd and zero in the interval's interior.
C
C   Modified August, 1991 and April, 1992 by Kaisheng Du and
C   R. Baker Kearfott to allow computation of negative integer powers.
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 NDUM-th power, of the interval A, where
C  NDUM is an integer.  It simulates directed roundings with the routine
C  RNDOUT; the interval result should contain the interval which would
C  have been obtained with exact interval arithmetic.  However, in
C  general it will not be the smallest possible machine-representable
C  such containing interval.  See the documentation in subroutine RNDOUT
C  for more detailed information.
C
C  This routine can clearly be made more efficient and more readable on
C  Fortran systems for which the intrinsic represented by A**N is
C  optimally accurate.  See the documentation in the subsidiary routine
C  RRPOWR for more information.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION AA(2), RESULT(2)
      INTEGER NDUM
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION B(2), TMP(2)
      LOGICAL EVEN ,L
      DOUBLE PRECISION TEMP
      INTEGER N
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 Fortran-supplied functions and subroutines --
C
C    MAX, MOD
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C    ERRTST, MULT, RNDOUT, RRPOWR
C
C***********************************************************************
C
C Beginning of executable statements --
C
C
C  Identifying code for this routine --
C
      IROUT = 7
      IERR = 0
C
      IF (AA(2).LT.AA(1)) THEN
         IERR = 1
         CALL ERRTST(AA)
      END IF
C
C  If N is equal to 0 then return the point interval [1,1] --
C
      N = NDUM
C
      IF (N.EQ.0) THEN
         RESULT(1) = ONE(1)
         RESULT(2) = ONE(2)
         RETURN
      END IF
C
C     If N = 1 then return AA.
C
      IF (N.EQ.1) THEN
         RESULT(1) = AA(1)
         RESULT(2) = AA(2)
         RETURN
      END IF
C
C  If N < 0, let N = -N.
C
      L = (N.LT.0)
      IF(L) THEN
         IF (AA(1).LT.ZERO(1) .AND. AA(2).GT.ZERO(2)) THEN
            IERR=8
            CALL ERRTST(AA)
            RESULT(1) = NEGINF
            RESULT(2) = POSINF
            RETURN
         END IF
         N = -N
      END IF
C
C  If N > 1, check cases for 0 in or not in the interval --
C
      EVEN = (MOD(N,2).EQ.0)
C
      IF (AA(1).GT.ZERO(2)) THEN
         CALL RRPOWR(AA(1),N,TMP)
         B(1) = TMP(1)
         CALL RRPOWR(AA(2),N,TMP)
         B(2) = TMP(2)
      ELSE IF (AA(2).LT.ZERO(1)) THEN
         IF (EVEN) THEN
            CALL RRPOWR(AA(2),N,TMP)
            B(1) = TMP(1)
            CALL RRPOWR(AA(1),N,TMP)
            B(2) = TMP(2)
         ELSE
            CALL RRPOWR(AA(1),N,TMP)
            B(1) = TMP(1)
            CALL RRPOWR(AA(2),N,TMP)
            B(2) = TMP(2)
         END IF
      ELSE IF (EVEN) THEN
         B(1) = ZERO(1)
         TMP(2) = -AA(1)
         CALL RNDOUT(TMP,.FALSE.,.TRUE.)
         TEMP = MAX(TMP(2),AA(2))
         CALL RRPOWR(TEMP,N,TMP)
         B(2) = TMP(2)
      ELSE
         CALL RRPOWR(AA(1),N,TMP)
         B(1) = TMP(1)
         CALL RRPOWR(AA(2),N,TMP)
         B(2) = TMP(2)
      END IF
C
      IF (L) THEN
         IF(B(1).GT.0D0)  THEN
            TEMP = B(1)
            B(1) = ONE(1)/B(2)
            B(2) = ONE(2)/TEMP
            CALL RNDOUT(B,.TRUE.,.TRUE.)
         ELSE
            TEMP = B(1)
            B(1) = ONE(2)/B(2)
            B(2) = ONE(1)/TEMP
            CALL RNDOUT(B,.TRUE.,.TRUE.)
         END IF
      END IF
C
      RESULT(1) = B(1)
      RESULT(2) = B(2)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE RRPOWR(AA,NDUM,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   November 14, 1992
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C***********************************************************************
C
C Called by --
C
C  POWER
C
C***********************************************************************
C
C Function --
C
C  This routine computes the NDUM-th power, of the interval point AA,
C  using simulated directed roundings to rigorously bound the roundoff
C  error.  The results is placed in RESULT.  The reason for this routine
C  is to allow rigor on machines on which exponentiation by an integer
C  is not accurate to within one digit in the last place.  The POWER
C  routine can certainly be made more efficient without this routine
C  if the Fortran intrinsic A**NDUM is optimally accurate.  This routine
C  can also be made more efficient by rewriting it in Fortran 90, and
C  making good use of recursion.  (See ``Programmer's Guide to Fortran
C  90, W. S. Brainerd, C. H. Goldberg, J. C. Adams, McGraw Hill, 1990
C  p. 222 ff.)
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION AA, RESULT(2)
      INTEGER NDUM
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION X(2)
      INTEGER I
C
C***********************************************************************
C
C Beginning of executable statements --
C
      X(1) = AA
      X(2) = AA
C
      RESULT(1) = X(1)
      RESULT(2) = X(2)
      IF (NDUM.EQ.1) RETURN
C
      DO 10 I = 2,NDUM,1
         CALL MULT(RESULT,X,RESULT)
   10 CONTINUE
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE RCOS(XA,VAL)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Kaisheng Du
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C
C***********************************************************************
C
C Function --
C
C  This routine is used to compute one of the endpoints of the
C  cosine of the interval X.  It is assumed that X is an interval of
C  very small width.  In particular, it is assumed that the length of
C  the interval is at most pi/16.  It is also assumed that the
C  absolute value of the argument is no greater than MAXX / pi,
C  where MAXX is the largest representable integer in the machine.
C  Since this routine doesn't check these conditions, it is assumed
C  that the calling routine has done so.
C
C  This routine is usually called by ISIN, which returns interval
C  values of the sin for general interval arguments.
C
C  This routine translates the interval to [0,pi/8], and possibly
C  uses a double angle formula twice.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION XA(2), VAL(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  XA    is the interval argument
C       (INPUT)
C
C  VAL   is the resulting interval cosine value
C        (OUTPUT)
C
C***********************************************************************
C
C Internal variable declarations --
C
      LOGICAL L2
      INTEGER N
      DOUBLE PRECISION  DN, ERTRM(2), TEM(2), TMP, V(2), VALT(2),
     *                  VSQ(2), X(2)
C
C***********************************************************************
C
C Common blocks --
C
      DOUBLE PRECISION MXULP, TOL0, TTINY2
      COMMON /MACH1/ MXULP, TTINY2, TOL0
C
C  The above 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
      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
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines
C
C     DBLE, MOD, NINT
C
C***********************************************************************
C
C Package-supplied functions and subroutines
C
C     ADD, MULT, POWER, RNDOUT, SCLMLT, SUB
C
C***********************************************************************
C
C Beginning of executable statements --
C
      X(1) = XA(1)
      X(2) = XA(2)
      IF (X(2).LT.0D0) THEN
         TMP = X(2)
         X(2) = -X(1)
         X(1) = -TMP
      END IF
C
C  Reduce the angle to between 0 and pi using translation by a
C  multiple of pi/2, making sure that the resulting interval is either
C  entirely positive or entirely negative.
C
      CALL MULT (X,A,TEM)
      N = NINT((TEM(1)+TEM(2))/2D0)
      DN = DBLE(N)
C
      CALL SCLMLT (DN,PI,V)
      CALL SUB (X,V,V)
C
      IF (V(2).LT.0D0) THEN
         TMP = V(2)
         V(2) = -V(1)
         V(1) = -TMP
         CALL RNDOUT(V,.TRUE.,.TRUE.)
      END IF
C
C  Further reduce the angle to between 0 and pi/8, using angle halving.
C
      L2 = V(2) .GE. PI8(1)
      IF (L2) THEN
         V(1) = V(1)/4D0
         V(2) = V(2)/4D0
         CALL RNDOUT(V,.TRUE.,.TRUE.)
      END IF
C
C  Compute value of the cosine of the reduced angle.
C
      CALL POWER(V,2,VSQ)
C
C  val  = one + vsq * (-od2f + vsq * (od4f + vsq * (-od6f
C        + vsq * (od8f + vsq * (-od10f + vsq * od12f)))))
C
      CALL MULT (VSQ, OD12F, VAL)
      CALL SUB (VAL, OD10F, VAL)
      CALL MULT (VSQ, VAL, VAL)
      CALL ADD (OD8F, VAL, VAL)
      CALL MULT (VSQ, VAL, VAL)
      CALL SUB (VAL, OD6F, VAL)
      CALL MULT (VSQ,VAL, VAL)
      CALL ADD (OD4F, VAL, VAL)
      CALL MULT (VSQ, VAL, VAL)
      CALL SUB (VAL, OD2F, VAL)
      CALL MULT (VSQ, VAL, VAL)
      CALL ADD (ONE, VAL, VAL)
C
C  valt = - v **14 * od14f
C
      CALL POWER (V, 14, VALT)
      CALL MULT (VALT, OD14F, VALT)
      TMP = VALT(2)
      VALT(2) = -VALT(1)
      VALT(1) = -TMP
      CALL RNDOUT(VALT,.TRUE.,.FALSE.)

      ERTRM(1) = VALT(1)
      ERTRM(2) = 0D0
C
      CALL ADD (VAL,ERTRM,VAL)
C
C  Use double angle formulas to get back to original angle in
C  [0, pi/2]. ( val <-- 8val^4 - 8val^2 + 1 )
C
      IF (L2) THEN
         CALL POWER(VAL, 2, TEM)
         CALL SUB  (TEM, ONE, VAL)
         CALL MULT (VAL, TEM, VAL)
         CALL MULT (VAL, EIGHT, VAL)
         CALL ADD  (ONE, VAL, VAL)
      END IF
C
C  Change the sign of the result if the shift was by an odd multiple of
C  pi --
C
      IF (MOD(N,2).NE.0) THEN
         TMP = VAL(1)
         VAL(1) = -VAL(2)
         VAL(2) = -TMP
      END IF
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE REXP(X,VAL)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Kaisheng Du
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C
C***********************************************************************
C
C Function --
C
C  This routine is used to compute the interval value exp(X).
C  It is assumed that the width of the interval X is small.
C
C  This routine is usually called by the routine IEXP,
C  to compute interval values of endpoints.
C
C  This routine translates the midpoint to [0,1/16], by first making
C  it positive,  then subtracting the integer part,  and finally
C  possibly dividing by 2 a number of times (To translate
C  back to the original interval, it then squares the result a number of
C  times, multiplies by e**N for some N, and possibly takes a
C  reciprocal.)
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(2), VAL(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TMP(2)
      INTEGER N
      DOUBLE PRECISION ERTRM(2), G(2), VALT(2)
C
C***********************************************************************
C
C Common block declarations --
C
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 MXULP, TOL0, TTINY2
      COMMON /MACH1/ MXULP, TTINY2, TOL0
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
C***********************************************************************
C
C Fortran-supplied functions and subroutines --
C
C     DBLE, DMAX1, DMIN1, NINT
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C     ADD, MULT, POWER, RNDOUT, SCLMLT
C
C***********************************************************************
C
C Beginning of executable statements --
C
C   Reduce the arguments with the formula X = G + N /4, where G is
C   between -.25 and .25 --
C
      CALL MULT (FOUR,X,TMP(1))
      N=NINT(TMP(1))
      TMP(1) = DBLE(N)
      TMP(2) = DBLE(N)
      CALL RNDOUT(TMP,.TRUE.,.TRUE.)
      CALL IDIV(TMP,FOUR,TMP)
      CALL SUB(X, TMP, G)
C
C  Compute Taylor series approximation --
C
c     val = one + g*(one + g*(od2f + g*(od3f + g*(od4f + g*(od5f +
c    *     g*(od6f + g*(od7f + g*(od8f +g*(od9f +g*(od10f +g*od11f
c    *     )))))))))) --
C
      CALL MULT (G,OD11F,TMP)
      CALL ADD (OD10F,TMP,TMP)
      CALL MULT (G,TMP,TMP)
      CALL ADD (OD9F,TMP,TMP)
      CALL MULT (G,TMP,TMP)
      CALL ADD (OD8F,TMP,TMP)
      CALL MULT (G,TMP,TMP)
      CALL ADD (OD7F,TMP,TMP)
      CALL MULT (G,TMP,TMP)
      CALL ADD (OD6F,TMP,TMP)
      CALL MULT (G,TMP,TMP)
      CALL ADD (OD5F,TMP,TMP)
      CALL MULT (G,TMP,TMP)
      CALL ADD (OD4F,TMP,TMP)
      CALL MULT (G,TMP,TMP)
      CALL ADD (OD3F,TMP,TMP)
      CALL MULT (G,TMP,TMP)
      CALL ADD (OD2F,TMP,TMP)
      CALL MULT (G,TMP,TMP)
      CALL ADD (ONE,TMP,TMP)
      CALL MULT (G,TMP,TMP)
      CALL ADD (ONE,TMP,VAL)
C
C  Compute the error term and add to the approximation --
C
C     VALT = G**12 * OD12F * A3
C
      CALL POWER (G,12,TMP)
      CALL MULT (TMP,OD12F,TMP)
      CALL MULT (TMP,A3,VALT)
C
      ERTRM(1) = DMIN1(0D0,VALT(1))
      ERTRM(2) = DMAX1(0D0,VALT(2))
C
      CALL ADD (VAL,ERTRM,VAL)
C
C  Translate back to original interval -- val=val*e14**n --
C
      CALL POWER (E14,N,TMP)
      CALL MULT (VAL,TMP,VAL)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE RLOG(XX,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Rebecca Yun
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C
C
C***********************************************************************
C
C Function --
C
C  This subroutine returns the interval value of the logarithm over the
C  interval XX in RESULT.  It is assumed that XX is relatively small,
C  and that it is possible.  The usual use of this routine is from the
C  routine ILOG, which checks the argument values and calls this routine
C  to get rigorous bounds on values at the end points.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION XX(2), RESULT(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      INTEGER K
      DOUBLE PRECISION ERRTRM(2), VAL(2), X(2), XI(2), Y(2)
C
C
C***********************************************************************
C
C Common blocks --
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
C***********************************************************************
C
C Fortran-supplied functions and subroutines --
C
C     DBLE, INT, LOG
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C   ADD, IDIV, MULT, POWER, RNDOUT, SUB
C
C***********************************************************************
C
C Beginning of executable statements --
C
      X(1) = XX(1)
      X(2) = XX(2)
C
C  Argument reduction: rewrite X <-- R^K * X
C               where  1/R <= X <= R = EXP(1/16)  ---
C
C   It is O.K. to use the Fortran LOG function here, since it it not
C   critical that X be exactly between 1/R and R ---
C
      IF (X(1).GE.ONE(1)) THEN
         K=INT(LOG((X(1)+X(2))/TWO(1))*SXTEEN(1))
      ELSE
         K = INT(LOG((X(1)+X(2))/TWO(1))*SXTEEN(1))-1
      END IF
C
      IF (K.LT.JTINY2) K=JTINY2
C
C    X <-- X/(R^K) ---
C
      CALL POWER(ESXTNT,K,VAL)
      CALL IDIV(X,VAL,X)
C
C  Initialize the values ---
C
      CALL SUB (X,ONE,Y)
C
      RESULT(1) = ZERO(1)
      RESULT(2) = ZERO(2)
C
C  Use Horner's scheme to evaluate the thirteenth degree Taylor
C  polynomial centered at one --
C
      CALL MULT(Y,THRTTH,RESULT)
      CALL SUB(RESULT,TWLVTH,RESULT)
      CALL MULT(Y,RESULT,RESULT)
      CALL ADD(ELEVTH,RESULT,RESULT)
      CALL MULT(Y,RESULT,RESULT)
      CALL SUB(RESULT,TENTH,RESULT)
      CALL MULT(Y,RESULT,RESULT)
      CALL ADD(NINTH,RESULT,RESULT)
      CALL MULT(Y,RESULT,RESULT)
      CALL SUB(RESULT,EIGHTH,RESULT)
      CALL MULT(Y,RESULT,RESULT)
      CALL ADD(SEVNTH,RESULT,RESULT)
      CALL MULT(Y,RESULT,RESULT)
      CALL SUB(RESULT,SIXTH,RESULT)
      CALL MULT(Y,RESULT,RESULT)
      CALL ADD(FIFTH,RESULT,RESULT)
      CALL MULT(Y,RESULT,RESULT)
      CALL SUB(RESULT,FOURTH,RESULT)
      CALL MULT(Y,RESULT,RESULT)
      CALL ADD(THIRD,RESULT,RESULT)
      CALL MULT(Y,RESULT,RESULT)
      CALL SUB(RESULT,OD2F,RESULT)
      CALL MULT(Y,RESULT,RESULT)
      CALL ADD(ONE,RESULT,RESULT)
      CALL MULT(Y,RESULT,RESULT)
C
C  Compute the error term ---
C
      IF (X(1).LT.ONE(1)) THEN
         XI(1) = X(1)
      ELSE
         XI(1) = ONE(1)
      END IF
C
C
      IF(X(2).GT.ONE(2)) THEN
         XI(2) = X(2)
      ELSE
         XI(2) = ONE(2)
      END IF
C
C    ERRTRM <-- Y^14 / [14*XI^14]
C
      CALL IDIV(Y,XI,ERRTRM)
      CALL POWER(ERRTRM,14,ERRTRM)
      CALL MULT(ERRTRM,FORTTH,ERRTRM)
C
C    RESULT <-- RESULT - ERRTRM
C
      CALL SUB (RESULT,ERRTRM,RESULT)
C
C  Undo the argument reduction: RESULT <-- RESULT+ K/16 --
C
      Y(1) = DBLE(K)
      Y(2) = Y(1)
      CALL RNDOUT(Y,.TRUE.,.TRUE.)
      CALL MULT(Y,SXTNTH,Y)
      CALL ADD (RESULT,Y,RESULT)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE RSQRT(X,VAL)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Kaisheng Du
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C
C
C***********************************************************************
C
C Function --
C
C  This routine returns the square root of the interval X in VAL.  It is
C  assumed that the interval X has relatively small width, and that its
C  lower bound is positive.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(2), VAL(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION AM1(2), TMP1(2), TMP2(2), ERRTRM(2)
      DOUBLE PRECISION XTMP(2), XMID(2), XOLD(2)
      DOUBLE PRECISION CMIN, CMAX, T
      LOGICAL L, LSMALL
      INTEGER ITER, N
C
C***********************************************************************
C
C Common blocks --
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
C***********************************************************************
C
C Fortran-supplied functions and subroutines --
C
C   LOG10, MAX, MIN, NINT
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C   ADD, IDIV, MULT, POWER, RNDOUT, SUB
C
C***********************************************************************
C
C Beginning of executable statements --
C
      XTMP(1) = X(1)
      XTMP(2) = X(2)
C
C  Reduce to between 1/5 and 5 by dividing by the appropriate power
C  of 10.  Use the Fortran logarithm function to get the power.  This is
C  O.K., provided the argument to this function is within range,
C  since it is not crucial that the argument be exactly in [1/5,5] --
C
      N = NINT( LOG10( (XTMP(1)+XTMP(2))/2D0 ) )
C
C  Treat values between the smallest number which can be reciprocated
C  and the smallest representable number (if these two values are
C  different) as a special case --
C
      IF (N.LT.ITINY2) THEN
         N=ITINY2
         XTMP(1) = TINY2
         XTMP(2) = TINY2
         LSMALL = .TRUE.
      ELSE
         LSMALL = .FALSE.
      END IF
      CALL POWER(TEN,-N,TMP1)
      CALL MULT(XTMP,TMP1,XTMP)
C
      L = .FALSE.
C
      IF (XTMP(2).LT.ONE(1)) L = .TRUE.
C
C  Take the reciprocal if initial value is greater than one, in
C  order to have stable behavior in the interval Newton method --
C
      IF (L .AND. .NOT.LSMALL) THEN
         T = ONE(1)/XTMP(2)
         XTMP(2) = ONE(2)/XTMP(1)
         XTMP(1) = T
         CALL RNDOUT(XTMP,.TRUE.,.TRUE.)
      END IF
C
C  The reduced argument can be less than zero, depending on the
C  rounding properties and size of the original argument --
C
      IF (XTMP(1).LE.ZERO(2)) THEN
         WRITE(6,*) 'IN RSQRT, ARGUMENT WAS DETERMINED TO BE TOO'
         WRITE(6,*) 'LARGE OR TOO SMALL.'
         STOP
      END IF
C
C  Make an initial estimate for the result using a degree two
C  Taylor series with remainder term, provided the interval is
C  sufficiently bounded away from zero
C  (i.e. V <-- 1 + 1/2 (c-1) - 1/8 (c-1)^2 + 1/32 I (c-1)^3 )
C
         VAL(1) = ONE(1)
         VAL(2) = ONE(2)
         CALL SUB(XTMP,ONE,AM1)
         CALL MULT(AM1,OD2F,TMP1)
         CALL ADD(VAL,TMP1,VAL)
         CALL POWER(AM1,2,TMP1)
         CALL MULT(TMP1,EIGHTH,TMP1)
         CALL SUB(VAL,TMP1,VAL)
C
         CMIN = MIN(XTMP(1),ONE(1))
         CMAX = MAX(XTMP(2),ONE(2))
         TMP1(1) = ONE(1)/CMAX
         TMP1(2) = ONE(2)/CMIN
         CALL RNDOUT(TMP1,.TRUE.,.TRUE.)
         CALL POWER(TMP1,2,TMP2)
         CALL POWER(TMP1,3,TMP1)
         ERRTRM(1) = MIN(TMP1(1),TMP2(1))
         ERRTRM(2) = MAX(TMP1(2),TMP2(2))
         CALL MULT(ERRTRM,SXTNTH,ERRTRM)
         CALL POWER(AM1,3,TMP1)
         CALL MULT(ERRTRM,TMP1,ERRTRM)
C
         CALL ADD(ERRTRM,VAL,VAL)
         VAL(1) = MAX(VAL(1),TINY2)
C
C  Do an interval Newton iteration until it becomes stationary --
C
      DO 10    ITER = 1,100
C
         XOLD(1) = VAL(1)
         XOLD(2) = VAL(2)
         XMID(1) = (VAL(1)+VAL(2))/2D0
         XMID(2) = XMID(1)
C
C       VAL <-- MID(VAL) - (MID(VAL)^2-A) / (TWO*VAL) --
C
         CALL POWER (XMID,2,TMP1)
         CALL SUB (TMP1,XTMP,TMP1)
         CALL MULT (TWO,VAL,TMP2)
         CALL IDIV(TMP1,TMP2,TMP2)
         CALL SUB (XMID,TMP2,VAL)
C
         VAL(1) = MAX(VAL(1),XOLD(1))
         VAL(2) = MIN(VAL(2),XOLD(2))
C
         IF (VAL(1).EQ.XOLD(1) .AND. VAL(2).EQ.XOLD(2)) GO TO 20
C
   10 CONTINUE
   20 CONTINUE
C
C  Take the reciprocal of the result, if the argument was
C  reciprocated --
C
      IF (L.AND. .NOT. LSMALL) THEN
         T = VAL(1)
         VAL(1) = ONE(1)/VAL(2)
         VAL(2) = ONE(2)/T
         CALL RNDOUT(VAL,.TRUE.,.TRUE.)
      END IF
C
C  Transform back to the original position by multiplying by the
C  appropriate power of the square root of 10 --
C
      CALL POWER(SQT10,N,TMP1)
      CALL MULT(TMP1,VAL,VAL)
C
      IF (LSMALL) VAL(1)=ZERO(1)
C
      RETURN
      END