c*** miscmach.f
      SUBROUTINE ERRTST(X)
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 prints error conditions which have been signalled
C  in other routines.  It halts execution, changes X, or returns,
C  depending on the particular error and on the way the error handling
C  control flags are set.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(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  X  is an interval which depends on the error set (I/O)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION 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 Beginning of executable statements --
C
      ISIG = 0
C
      IF (IERR.EQ.1) THEN
         ISIG = 2
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*)
     *    'SEVERITY 2 ERROR:  LOWER BOUND ON INTERVAL IS GREATER',
     *    ' THAN UPPER BOUND'
            WRITE(IERPUN,100) 'X:',X
         END IF
         T = X(1)
         X(1) = X(2)
         X(2) = T
      ELSE IF (IERR.EQ.2) THEN
         ISIG = 1
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*)
     *    'SEVERITY 1 ERROR: LOWER BOUND ON X FOR IEXP WOULD UNDERFLOW'
         END IF
      ELSE IF (IERR.EQ.3) THEN
         ISIG = 2
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*)
     *    'SEVERITY 2 ERROR: LOWER BOUND ON X FOR IEXP WOULD OVERFLOW'
         END IF
      ELSE IF (IERR.EQ.4) THEN
         ISIG = 2
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*)
     *    'SEVERITY 2 ERROR: UPPER BOUND ON X FOR IEXP WOULD OVERFLOW'
         END IF
      ELSE IF (IERR.EQ.5) THEN
         ISIG = 0
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*)
     *    'WARNING: Loss of accuracy in trig function due to the',
     *    ' argument range.'
            WRITE(IERPUN,*) 'X:',X
         END IF
      ELSE IF (IERR.EQ.6) THEN
         ISIG = 2
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*)
     *    'SEVERITY 2 ERROR: ZERO IN DENOMINATOR IN ORDINARY INTERVAL',
     *    ' DIVISION.'
            WRITE(IERPUN,*) 'DENOMINATOR:',X
         END IF
      ELSE IF (IERR.EQ.7) THEN
         ISIG = 3
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*)
     *    'SEVERITY 3 ERROR: ARGUMENT CONTAINS ZERO IN ELEMENTARY',
     *    ' FUNCTION.'
            WRITE(IERPUN,*) 'ARGUMENT:',X
         END IF
      ELSE IF (IERR.EQ.8) THEN
         ISIG = 2
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*)
     *    'SEVERITY 2 ERROR: NEGATIVE POWER OF A ZERO-CONTAINING'
            WRITE(IERPUN,*)
     *    'INTERVAL IS UNDEFINED IN ORDINARY INTERVAL ARITHMETIC.'
            WRITE(IERPUN,*) 'ARGUMENT:',X
         END IF
      ELSE IF (IERR.EQ.9) THEN
         ISIG = 3
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*)
     *    'SEVERITY 3 ERROR: ARGUMENT TO INVERSE TRIG FUNCTION'
            WRITE(IERPUN,*)
     *    'MAY CONTAIN NUMBERS WHICH ARE GREATER THAN 1.'
            WRITE(IERPUN,*) 'ARGUMENT:',X
         END IF
      ELSE IF (IERR.EQ.10) THEN
         ISIG = 3
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*)
     *    'SEVERITY 3 ERROR: ARGUMENT TO INVERSE TRIG FUNCTION'
            WRITE(IERPUN,*)
     *    'MAY CONTAIN NUMBERS WHICH ARE LESS THAN -1.'
            WRITE(IERPUN,*) 'ARGUMENT:',X
         END IF
      ELSE IF (IERR.EQ.11) THEN
         ISIG = 3
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*)
     *    'SEVERITY 3 ERROR: APPROXIMATING SERIES FOR REDUCED ARGUMENT',
     *    ' DID NOT CONVERGE.'
            WRITE(IERPUN,*) 'CURRENT TERM IN SERIES:',X
         END IF
      ELSE IF (IERR.EQ.12) THEN
         ISIG = 3
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*) 'INTERNAL ERROR;  NEGATIVE ARGUMENT:'
            WRITE(IERPUN,*) X
            WRITE(IERPUN,*) 'POSSIBLY, AN ARRAY WAS DIMENSIONED',
     *                   ' INCORRECTLY OR MEMORY WAS NOT PROPERLY'
            WRITE(IERPUN,*) 'ALLOCATED.  IF THE ERROR PERSISTS, THEN:'
            WRITE(IERPUN,*) 'Contact R. Baker Kearfott, Dept. Math.,'
            WRITE(IERPUN,*) 'USL Box 4-1010, Lafayette, LA 70504-1010.'
            WRITE(IERPUN,*) 'email: rbk@usl.edu'
         END IF
      ELSE IF (IERR.EQ.13) THEN
         ISIG = 0
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*)
     *        'Warning:  Disjoint intervals in intersection.'
            X(1) = POSINF
            X(2) = NEGINF
         END IF
      ELSE IF (IERR.EQ.14) THEN
         ISIG = 1
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*)
     *        'Power function with operands containing 0^0 occured.'
            WRITE(IERPUN,*)
     *        'Result is sent to [NEGINF,POSINF].'
         END IF
      ELSE IF (IERR.EQ.15) THEN
         ISIG = 1
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*)
     *        'Power function with operands of the form'
            WRITE(IERPUN,*)
     *        '[0,pos]^[pos,pos] occurred. Result is set to'
            WRITE(IERPUN,*)
     *        '[pos,POSINF].'
         END IF
      ELSE
         ISIG=3
         IF (IPRTCL.LE.ISIG) THEN
            WRITE(IERPUN,*) 'INTERNAL ERROR;  UNKNOWN ERROR TYPE,',IERR,
     *                      ' IN ERRTST.'
            WRITE(IERPUN,*) 'POSSIBLY, AN ARRAY WAS DIMENSIONED',
     *                   ' INCORRECTLY OR MEMORY WAS NOT PROPERLY'
            WRITE(IERPUN,*) 'ALLOCATED.  IF THE ERROR PERSISTS, THEN:'
            WRITE(IERPUN,*) 'Contact R. Baker Kearfott, Dept. Math.,'
            WRITE(IERPUN,*) 'USL Box 4-1010, Lafayette, LA 70504-1010.'
            WRITE(IERPUN,*) 'email: rbk@usl.edu'
         END IF
      END IF

C
      IF (IPRTCL.LE.ISIG) THEN
         IF (IROUT.EQ.1) THEN
            WRITE(IERPUN,*) 'Error occurred in routine IEXP.'
         ELSE IF (IROUT.EQ.2) THEN
            WRITE(IERPUN,*) 'Error occurred in routine ICOS.'
         ELSE IF (IROUT.EQ.3) THEN
            WRITE(IERPUN,*) 'Error occurred in routine IDIV.'
         ELSE IF (IROUT.EQ.4) THEN
            WRITE(IERPUN,*) 'Error occurred in routine ILOG.'
         ELSE IF (IROUT.EQ.5) THEN
            WRITE(IERPUN,*) 'Error occurred in routine ISQRT.'
         ELSE IF (IROUT.EQ.6) THEN
            WRITE(IERPUN,*) 'Error occurred in routine ISINH.'
         ELSE IF (IROUT.EQ.7) THEN
            WRITE(IERPUN,*) 'Error occurred in routine POWER.'
         ELSE IF (IROUT.EQ.8) THEN
            WRITE(IERPUN,*) 'Error occurred in routine IACOS.'
         ELSE IF (IROUT.EQ.9) THEN
            WRITE(IERPUN,*) 'Error occurred in routine IASIN.'
         ELSE IF (IROUT.EQ.10) THEN
            WRITE(IERPUN,*) 'Error occurred in routine ASNSER, ',
     *                      'an auxiliary routine called by IASIN.'
         ELSE IF (IROUT.EQ.11) THEN
            WRITE(IERPUN,*) 'Error occurred in routine IATAN.'
         ELSE IF (IROUT.EQ.12) THEN
            WRITE(IERPUN,*) 'Error occurred in routine ATNSER, ',
     *                      'an auxiliary routine called by IATAN.'
         ELSE IF (IROUT.EQ.13) THEN
            WRITE(IERPUN,*) 'Error occurred in routine ICAP, ',
     *                      'a utility function routine.'
         ELSE IF (IROUT.EQ.14) THEN
            WRITE(IERPUN,*) 'Error occurred in routine IIPOWR.'
         ELSE IF (IROUT.EQ.15) THEN
            WRITE(IERPUN,*) 'Error occurred in routine IVL2, ',
     *                      'a utility function routine.'
         ELSE
            ISIG = 3
            WRITE(IERPUN,*) 'INTERNAL ERROR;  UNKNOWN ROUTINE NUMBER,',
     *                      IROUT, ' IN ERRTST.'
            WRITE(IERPUN,*) 'POSSIBLY, AN ARRAY WAS DIMENSIONED',
     *                   ' INCORRECTLY OR MEMORY WAS NOT PROPERLY'
            WRITE(IERPUN,*) 'ALLOCATED.  IF THE ERROR PERSISTS, THEN:'
            WRITE(IERPUN,*) 'Contact R. Baker Kearfott, Dept. Math.,'
            WRITE(IERPUN,*) 'USL Box 4-1010, Lafayette, LA 70504-1010.'
            WRITE(IERPUN,*) 'email: rbk@usl.edu'
         END IF
      END IF
C
      IF (ISIG.GE.ISEVER) STOP
C
      RETURN
C
  100 FORMAT(1X,A,3X,D23.15,2X,D23.15)
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SIMINI
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   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 Called by -- GENBIS
C
C***********************************************************************
C
C Function --
C
C  This routine sets certain machine parameters used to simulate
C  directed roundings in a reasonably transportable way.  In
C  particular, it sets the amount by which to decrease the left endpoint
C  and increase the right endpoint of an interval computed using usual
C  floating point arithmetic to guarantee that the resulting interval
C  will contain the result which would have been obtained with true
C  interval arithmetic.
C
C
C  In addition to setting certain parameters, certain constants, such as
C  pi and the natural logarithm base e, are set here.  The data
C  for the representation for these constants which appears here
C  should be accurate at least to the number of significant figures
C  present in the machine.  If it is possible to input a binary form
C  for the double precision representation of this data, then that
C  form should be given to be correct to all representable digits.
C
C  This routine assumes that the four elementary floating point
C  operations, and unary negation, will give results with a maximum
C  error of one ULP (unit in the last place).  If this is not so, change
C  the value of MAXERR in the data statement below to the maximum number
C  of ULP's by which a floating point result can differ from the true
C  result (for '+', '-', '*', '/' and conversion from integer
C  to double precision).
C
C  An additional assumption is that the standard routines MIN and MAX
C  return exact values corresponding to one of their arguments.
C
C  When determining the maximum error of the result A op B, where
C  A and B are floating point numbers, we assume that A and B are
C  represented exactly.  For example, if A and B are almost equal,
C  then it is not unreasonable to assume that A - B, where the
C  subtraction is a floating point subtraction, is within a few
C  units of the last place of the true result.
C
C  Throughout the elementary function routines, it is assumed that
C  storing the double precision expressions 0D0 leads to an exact
C  floating point representation of 0. It is also assumed that a
C  floating point assignment statement (such as A=B) causes exactly the
C  same value to be in A as in B.
C
C  Note:  On some machines, an underflow naturally occurs when this
C         routine is executed.  (It would happen in the computation of
C         TINY2.)  This is not a problem.
C
C  If SIMINI is installed correctly, then the conclusions this
C  package prints out will have mathematical rigor.
C
C***********************************************************************
C
C Common block declarations --
C
      DOUBLE PRECISION MXULP, TTINY2, TOL0
      COMMON /MACH1/ MXULP, TTINY2, TOL0
C
C  This common block holds machine parameters which are set here
C  and used in RNDOUT.
C
C  Variable descriptions
C
C  MXULP       (machine epsilon)
C                * (maximum error in ULP's of the floating pt. op's)
C
C  TTINY2      2 * (smallest representable positive machine number)
C                * (maximum error in ULP's of the floating pt. op's)
C
C  TOL0        TTINY2 / MXULP
C
      DOUBLE PRECISION TINY, TEST
      COMMON /MACH2/ TINY, TEST
C
C  Common block MACH2 stores machine constants used in RNDOUT.  TINY
C  is the smallest representable machine number, while TEST is the
C  smallest number which can be safely rounded to something other than
C  zero.
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 -- DBLE, INT, MAX
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  D1MACH, I1MACH (the SLATEC routines for machine constants)
C
C  RNDOUT
C
      DOUBLE PRECISION D1MACH
      INTEGER I1MACH
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C I/O functions -- none
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TMP(2)
C
C***********************************************************************
C
C Internal constant declarations --
C
      INTEGER MAXERR
      DATA MAXERR/1/
C
C***********************************************************************
C
C Internal constant descriptions --
C
C  MAXERR      is the maximum number of ULP's (units in the last
C              place) by which a result of one of the floating point
C              operations (+, -, *, /, ** N) can differ from the
C              true result.  (See explanation above.)
C
C********** WARNING: The value of MAXERR is machine dependent and
C                    must be manually set.
C
C***********************************************************************
C
C Beginning of executable statements --
C
      TINY = D1MACH(1)
      MXULP = DBLE(MAXERR) * D1MACH(4)
      TEST = TINY/(1D0-3D0*MXULP)
      TTINY2 = 2D0 * DBLE(MAXERR) * TINY
      TOL0 = TTINY2 / MXULP
C
      TMP(1) = 1D0/D1MACH(2)
      TMP(2) = TMP(1)
      CALL RNDOUT (TMP,.TRUE.,.TRUE.)
      TINY2 = MAX(TINY,TMP(2))
      CBTEP = (6D0*D1MACH(4))**(1D0/3D0)
      ITINY2 = INT(LOG10(TINY2))
      JTINY2 = 16*INT(LOG(TINY2))
C
      MAXX = DBLE(I1MACH(9))
C
      MXLGM1 = LOG(D1MACH(2)) -LOG(4D0) - 1D0
      NEGINF = -D1MACH(2)
      POSINF =  D1MACH(2)
C
      A(1) = 0.31830988618379067153776752674502864D+00
      A(2) = 0.31830988618379067153776752674502864D+00
      CALL RNDOUT(A,.TRUE.,.TRUE.)
C
      PI(1) = 0.31415926535897932384626433832795028D+01
      PI(2) = 0.31415926535897932384626433832795028D+01
      CALL RNDOUT(PI,.TRUE.,.TRUE.)

      E(1) =  0.27182818284590452353602874713526625D+01
      E(2) =  0.27182818284590452353602874713526625D+01
      CALL RNDOUT(E,.TRUE.,.TRUE.)
C
      PI2(1) = 0.15707963267948966192313216916397499D+01
      PI2(2) = 0.15707963267948966192313216916397499D+01
      CALL RNDOUT(PI2,.TRUE.,.TRUE.)
C
      PI3(1) = 1.0471975511965977461542144610931676D0
      PI3(2) = 1.0471975511965977461542144610931676D0
      CALL RNDOUT(PI3,.TRUE.,.TRUE.)
C
      PI4(1) = 0.7853981633974483096156608458198757D0
      PI4(2) = 0.7853981633974483096156608458198757D0
      CALL RNDOUT(PI4,.TRUE.,.TRUE.)
C
      PI6(1) = 0.5235987755982988730771072305465838D0
      PI6(2) = 0.5235987755982988730771072305465838D0
      CALL RNDOUT(PI6,.TRUE.,.TRUE.)
C
      PI8(1) = 0.39269908169872415480783042290993785D+00
      PI8(2) = 0.39269908169872415480783042290993785D+00
      CALL RNDOUT(PI8,.TRUE.,.TRUE.)
C
      E14(1) = 0.12840254166877414840734205680624368D+01
      E14(2) = 0.12840254166877414840734205680624368D+01
      CALL RNDOUT(E14,.TRUE.,.TRUE.)
C
      A3(1) = 0.11331484530668263168290072278117947D+01
      A3(2) = 0.11331484530668263168290072278117947D+01
      CALL RNDOUT(A3,.TRUE.,.TRUE.)
C
      ONE  (1) = 1D0
      OD2F (1) = 1D0 /           2D0
      OD3F (1) = 1D0 /           6D0
      OD4F (1) = 1D0 /          24D0
      OD5F (1) = 1D0 /         120D0
      OD6F (1) = 1D0 /         720D0
      OD7F (1) = 1D0 /        5040D0
      OD8F (1) = 1D0 /       40320D0
      OD9F (1) = 1D0 /      362880D0
      OD10F(1) = 1D0 /     3628800D0
      OD11F(1) = 1D0 /    39916800D0
      OD12F(1) = 1D0 /   479001600D0
      OD14F(1) = 1D0 / 87178291200D0
C
      ONE  (2) = ONE  (1)
      OD2F (2) = OD2F (1)
      OD3F (2) = OD3F (1)
      OD4F (2) = OD4F (1)
      OD5F (2) = OD5F (1)
      OD6F (2) = OD6F (1)
      OD7F (2) = OD7F (1)
      OD8F (2) = OD8F (1)
      OD9F (2) = OD9F (1)
      OD10F(2) = OD10F(1)
      OD11F(2) = OD11F(1)
      OD12F(2) = OD12F(1)
      OD14F(2) = OD14F(1)
C
      CALL RNDOUT(  ONE,.TRUE.,.TRUE.)
      CALL RNDOUT( OD2F,.TRUE.,.TRUE.)
      CALL RNDOUT( OD4F,.TRUE.,.TRUE.)
      CALL RNDOUT( OD6F,.TRUE.,.TRUE.)
      CALL RNDOUT( OD8F,.TRUE.,.TRUE.)
      CALL RNDOUT(OD10F,.TRUE.,.TRUE.)
      CALL RNDOUT(OD12F,.TRUE.,.TRUE.)
      CALL RNDOUT(OD14F,.TRUE.,.TRUE.)
C
      EIGHT(1) = 8D0
      EIGHT(2) = 8D0
      CALL RNDOUT(EIGHT,.TRUE.,.TRUE.)
C
      ZERO(1) = 0D0
      ZERO(2) = 0D0
C
      TWO(1) = 2D0
      TWO(2) = 2D0
      CALL RNDOUT(TWO,.TRUE.,.TRUE.)
C
      THREE(1) = 3D0
      THREE(2) = 3D0
      CALL RNDOUT(THREE,.TRUE.,.TRUE.)
C
      FOUR(1) = 4D0
      FOUR(2) = 4D0
      CALL RNDOUT(FOUR,.TRUE.,.TRUE.)
C
      NINE(1) = 9D0
      NINE(2) = 9D0
      CALL RNDOUT(NINE,.TRUE.,.TRUE.)
C
      SXTNTH(1) = 1D0/16D0
      SXTNTH(2) = SXTNTH(1)
      CALL RNDOUT(SXTNTH,.TRUE.,.TRUE.)
C
      TEN(1) = 10D0
      TEN(2) = 10D0
      CALL RNDOUT(TEN,.TRUE.,.TRUE.)
C
      TWOT7(1) = 27D0
      TWOT7(2) = 27D0
      CALL RNDOUT(TWOT7,.TRUE.,.TRUE.)
C
      SQT10(1) = 3.1622776601683793319988935444327185D0
      SQT10(2) = 3.1622776601683793319988935444327185D0
      CALL RNDOUT(SQT10,.TRUE.,.TRUE.)
C
      ESXTNT(1) = 1.0644944589178594295633905946428909D0
      ESXTNT(2) = 1.0644944589178594295633905946428909D0
      CALL RNDOUT(ESXTNT,.TRUE.,.TRUE.)
C
      SXTEEN(1) = 16D0
      SXTEEN(2) = 16D0
      CALL RNDOUT(SXTEEN, .TRUE.,.TRUE.)
C
      THIRD(1) = 1D0/3D0
      THIRD(2) = 1D0/3D0
      CALL RNDOUT(THIRD,.TRUE.,.TRUE.)
C
      FOURTH(1) = .25D0
      FOURTH(2) = .25D0
      CALL RNDOUT(FOURTH,.TRUE.,.TRUE.)
C
      FIFTH(1) = .2D0
      FIFTH(2) = .2D0
      CALL RNDOUT(FIFTH,.TRUE.,.TRUE.)
C
      SIXTH(1) = 1D0/6D0
      SIXTH(2) = 1D0/6D0
      CALL RNDOUT(SIXTH,.TRUE.,.TRUE.)
C
      SEVNTH(1) = 1D0/7D0
      SEVNTH(2) = 1D0/7D0
      CALL RNDOUT(SEVNTH,.TRUE.,.TRUE.)
C
      EIGHTH(1) = .125D0
      EIGHTH(2) = .125D0
      CALL RNDOUT(EIGHTH,.TRUE.,.TRUE.)
C
      NINTH(1) = 1D0/9D0
      NINTH(2) = 1D0/9D0
      CALL RNDOUT(NINTH,.TRUE.,.TRUE.)
C
      TENTH(1) = .1D0
      TENTH(2) = .1D0
      CALL RNDOUT(TENTH,.TRUE.,.TRUE.)
C
      ELEVTH(1) = 1D0/11D0
      ELEVTH(2) = 1D0/11D0
      CALL RNDOUT(ELEVTH,.TRUE.,.TRUE.)
C
      TWLVTH(1) = 1D0/12D0
      TWLVTH(2) = 1D0/12D0
      CALL RNDOUT(TWLVTH,.TRUE.,.TRUE.)
C
      THRTTH(1) = 1D0/13D0
      THRTTH(2) = 1D0/13D0
      CALL RNDOUT(THRTTH,.TRUE.,.TRUE.)
C
      FORTTH(1) = 1D0/14D0
      FORTTH(2) = 1D0/14D0
      CALL RNDOUT(FORTTH,.TRUE.,.TRUE.)
C
      TT7TH(1) = 1D0/27D0
      TT7TH(2) = 1D0/27D0
      CALL RNDOUT(TT7TH,.TRUE.,.TRUE.)
C
      SQT3(1) = 1.7320508075688772935274463415058724D0
      SQT3(2) = 1.7320508075688772935274463415058724D0
      CALL RNDOUT(SQT3,.TRUE.,.TRUE.)
C
      ODSQT3(1) = 0.57735026918962576450914878050195746D0
      ODSQT3(2) = 0.57735026918962576450914878050195746D0
      CALL RNDOUT(ODSQT3,.TRUE.,.TRUE.)
      CALL SCLMLT(1D0+100D0*MXULP,ODSQT3,ODSQT3)
C
C  Set default values for the error checking routine --
C
      ISIG = 0
      IERR = 0
      IPRTCL = 0
      ISEVER = 3
      IERPUN = 6
C
      RETURN
      END
