C      ALGORITHM 737, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 20, NO. 4, DECEMBER, 1994, P. 447-459.

This package consists of 5 files

1) D1I1MACH.F which contains the routines D1MACH and I1MACH from the 
   SLATEC package.

2) BASICOPS.F which contains the routines for interval inclusions for
   the elementary arithmetic operations.  These are ADD, CANCEL, IDIV,
   MULT, RNDOUT, SCLADD, SCLMLT and SUB.  Documentation for these 
   routines appears in their prologues.

3) UTILFUNS.F which contains basic utility functions for operations
   such as intersection of intervals, checking inclusion of one interval
   in another, and so on.  The routines are ICAP, IDISJ, IHULL, IILEI,
   IILTI, IINF, IMID, IMIG, INEG, INTABS, IRLEI, IRLTI, ISUP, IVL1, IVL2
   and IWID.
   
4) ELEMFUNS.F which contains routines for the elementary functions.
   The routines are IACOS, IACOT, IASIN, IATAN, ICOS, IEXP, IIPOWR, 
   ILOG, ISIN, ISINH, ISQRT, AND POWER.

5) MISCMACH.F which contains the routine ERRTST, used to print
   information about error conditions.  It also contains the routine
   SIMINI, in which global variables specifying interval inclusions for
   mathematical constants are set.  Some assumptions on the arithmetic 
   are given in the prologue to SIMINI.

C*** d1i1mach.f
Caveat receptor.  (Jack) dongarra@anl-mcs, (Eric Grosse) research!ehg
Compliments of netlib   Wed Mar 11 09:10:56 CST 1987
      DOUBLE PRECISION FUNCTION D1MACH(I)
C
C  DOUBLE-PRECISION MACHINE CONSTANTS
C
C  D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
C
C  D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
C
C  D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
C
C  D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
C
C  D1MACH( 5) = LOG10(B)
C
C  TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT,
C  THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY
C  REMOVING THE C FROM COLUMN 1.
C  ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED.
C  (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.)
C
C  FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), ONE OF THE FIRST
C  TWO SETS OF CONSTANTS BELOW SHOULD BE APPROPRIATE.
C
C  WHERE POSSIBLE, OCTAL OR HEXADECIMAL CONSTANTS HAVE BEEN USED
C  TO SPECIFY THE CONSTANTS EXACTLY, WHICH HAS IN SOME CASES
C  REQUIRED THE USE OF EQUIVALENT INTEGER ARRAYS.
C
      INTEGER SMALL(4)
      INTEGER LARGE(4)
      INTEGER RIGHT(4)
      INTEGER DIVER(4)
      INTEGER LOG10(4)
C
      DOUBLE PRECISION DMACH(5)
C
      EQUIVALENCE (DMACH(1),SMALL(1))
      EQUIVALENCE (DMACH(2),LARGE(1))
      EQUIVALENCE (DMACH(3),RIGHT(1))
      EQUIVALENCE (DMACH(4),DIVER(1))
      EQUIVALENCE (DMACH(5),LOG10(1))
C
C     MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T
C     3B SERIES AND MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T
C     PC 7300), IN WHICH THE MOST SIGNIFICANT BYTE IS STORED FIRST.
C
C      DATA SMALL(1),SMALL(2) /    1048576,          0 /
C      DATA LARGE(1),LARGE(2) / 2146435071,         -1 /
C      DATA RIGHT(1),RIGHT(2) / 1017118720,          0 /
C      DATA DIVER(1),DIVER(2) / 1018167296,          0 /
C      DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 /
C
C     MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES AND 8087-BASED
C     MICROS, SUCH AS THE IBM PC AND AT&T 6300, IN WHICH THE LEAST
C     SIGNIFICANT BYTE IS STORED FIRST.
C
       DATA SMALL(1),SMALL(2) /          0,    1048576 /
       DATA LARGE(1),LARGE(2) /         -1, 2146435071 /
       DATA RIGHT(1),RIGHT(2) /          0, 1017118720 /
       DATA DIVER(1),DIVER(2) /          0, 1018167296 /
       DATA LOG10(1),LOG10(2) / 1352628735, 1070810131 /
C
C     MACHINE CONSTANTS FOR AMDAHL MACHINES.
C
C      DATA SMALL(1),SMALL(2) /    1048576,          0 /
C      DATA LARGE(1),LARGE(2) / 2147483647,         -1 /
C      DATA RIGHT(1),RIGHT(2) /  856686592,          0 /
C      DATA DIVER(1),DIVER(2) /  873463808,          0 /
C      DATA LOG10(1),LOG10(2) / 1091781651, 1352628735 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C      DATA SMALL(1) / ZC00800000 /
C      DATA SMALL(2) / Z000000000 /
C
C      DATA LARGE(1) / ZDFFFFFFFF /
C      DATA LARGE(2) / ZFFFFFFFFF /
C
C      DATA RIGHT(1) / ZCC5800000 /
C      DATA RIGHT(2) / Z000000000 /
C
C      DATA DIVER(1) / ZCC6800000 /
C      DATA DIVER(2) / Z000000000 /
C
C      DATA LOG10(1) / ZD00E730E7 /
C      DATA LOG10(2) / ZC77800DC0 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
C
C      DATA SMALL(1) / O1771000000000000 /
C      DATA SMALL(2) / O0000000000000000 /
C
C      DATA LARGE(1) / O0777777777777777 /
C      DATA LARGE(2) / O0007777777777777 /
C
C      DATA RIGHT(1) / O1461000000000000 /
C      DATA RIGHT(2) / O0000000000000000 /
C
C      DATA DIVER(1) / O1451000000000000 /
C      DATA DIVER(2) / O0000000000000000 /
C
C      DATA LOG10(1) / O1157163034761674 /
C      DATA LOG10(2) / O0006677466732724 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
C
C      DATA SMALL(1) / O1771000000000000 /
C      DATA SMALL(2) / O7770000000000000 /
C
C      DATA LARGE(1) / O0777777777777777 /
C      DATA LARGE(2) / O7777777777777777 /
C
C      DATA RIGHT(1) / O1461000000000000 /
C      DATA RIGHT(2) / O0000000000000000 /
C
C      DATA DIVER(1) / O1451000000000000 /
C      DATA DIVER(2) / O0000000000000000 /
C
C      DATA LOG10(1) / O1157163034761674 /
C      DATA LOG10(2) / O0006677466732724 /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C
C      DATA SMALL(1) / 00604000000000000000B /
C      DATA SMALL(2) / 00000000000000000000B /
C
C      DATA LARGE(1) / 37767777777777777777B /
C      DATA LARGE(2) / 37167777777777777777B /
C
C      DATA RIGHT(1) / 15604000000000000000B /
C      DATA RIGHT(2) / 15000000000000000000B /
C
C      DATA DIVER(1) / 15614000000000000000B /
C      DATA DIVER(2) / 15010000000000000000B /
C
C      DATA LOG10(1) / 17164642023241175717B /
C      DATA LOG10(2) / 16367571421742254654B /
C
C     MACHINE CONSTANTS FOR CONVEX C-1
C
C      DATA SMALL(1),SMALL(2) / '00100000'X, '00000000'X /
C      DATA LARGE(1),LARGE(2) / '7FFFFFFF'X, 'FFFFFFFF'X /
C      DATA RIGHT(1),RIGHT(2) / '3CC00000'X, '00000000'X /
C      DATA DIVER(1),DIVER(2) / '3CD00000'X, '00000000'X /
C      DATA LOG10(1),LOG10(2) / '3FF34413'X, '509F79FF'X /
C
C     MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3.
C
C      DATA SMALL(1) / 201354000000000000000B /
C      DATA SMALL(2) / 000000000000000000000B /
C
C      DATA LARGE(1) / 577767777777777777777B /
C      DATA LARGE(2) / 000007777777777777776B /
C
C      DATA RIGHT(1) / 376434000000000000000B /
C      DATA RIGHT(2) / 000000000000000000000B /
C
C      DATA DIVER(1) / 376444000000000000000B /
C      DATA DIVER(2) / 000000000000000000000B /
C
C      DATA LOG10(1) / 377774642023241175717B /
C      DATA LOG10(2) / 000007571421742254654B /
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE -
C     STATIC DMACH(5)
C
C      DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/
C      DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/
C      DATA LOG10/40423K,42023K,50237K,74776K/
C
C     MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7
C
C      DATA SMALL(1),SMALL(2) / '20000000, '00000201 /
C      DATA LARGE(1),LARGE(2) / '37777777, '37777577 /
C      DATA RIGHT(1),RIGHT(2) / '20000000, '00000333 /
C      DATA DIVER(1),DIVER(2) / '20000000, '00000334 /
C      DATA LOG10(1),LOG10(2) / '23210115, '10237777 /
C
C     MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES.
C
C      DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 /
C      DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 /
C      DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 /
C      DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 /
C      DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /
C
C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C     THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86.
C
C      DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 /
C      DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF /
C      DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 /
C      DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 /
C      DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF /
C
C     MACHINE CONSTANTS FOR THE INTERDATA 8/32
C     WITH THE UNIX SYSTEM FORTRAN 77 COMPILER.
C
C     FOR THE INTERDATA FORTRAN VII COMPILER REPLACE
C     THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S.
C
C      DATA SMALL(1),SMALL(2) / Z'00100000', Z'00000000' /
C      DATA LARGE(1),LARGE(2) / Z'7EFFFFFF', Z'FFFFFFFF' /
C      DATA RIGHT(1),RIGHT(2) / Z'33100000', Z'00000000' /
C      DATA DIVER(1),DIVER(2) / Z'34100000', Z'00000000' /
C      DATA LOG10(1),LOG10(2) / Z'41134413', Z'509F79FF' /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
C
C      DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 /
C      DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 /
C      DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 /
C      DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 /
C      DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
C
C      DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 /
C      DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 /
C      DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 /
C      DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 /
C      DATA LOG10(1),LOG10(2) / "177464202324, "047674776746 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C      DATA SMALL(1),SMALL(2) /    8388608,           0 /
C      DATA LARGE(1),LARGE(2) / 2147483647,          -1 /
C      DATA RIGHT(1),RIGHT(2) /  612368384,           0 /
C      DATA DIVER(1),DIVER(2) /  620756992,           0 /
C      DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /
C
C      DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 /
C      DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 /
C      DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 /
C      DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 /
C      DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
C     16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C      DATA SMALL(1),SMALL(2) /    128,      0 /
C      DATA SMALL(3),SMALL(4) /      0,      0 /
C
C      DATA LARGE(1),LARGE(2) /  32767,     -1 /
C      DATA LARGE(3),LARGE(4) /     -1,     -1 /
C
C      DATA RIGHT(1),RIGHT(2) /   9344,      0 /
C      DATA RIGHT(3),RIGHT(4) /      0,      0 /
C
C      DATA DIVER(1),DIVER(2) /   9472,      0 /
C      DATA DIVER(3),DIVER(4) /      0,      0 /
C
C      DATA LOG10(1),LOG10(2) /  16282,   8346 /
C      DATA LOG10(3),LOG10(4) / -31493, -12296 /
C
C      DATA SMALL(1),SMALL(2) / O000200, O000000 /
C      DATA SMALL(3),SMALL(4) / O000000, O000000 /
C
C      DATA LARGE(1),LARGE(2) / O077777, O177777 /
C      DATA LARGE(3),LARGE(4) / O177777, O177777 /
C
C      DATA RIGHT(1),RIGHT(2) / O022200, O000000 /
C      DATA RIGHT(3),RIGHT(4) / O000000, O000000 /
C
C      DATA DIVER(1),DIVER(2) / O022400, O000000 /
C      DATA DIVER(3),DIVER(4) / O000000, O000000 /
C
C      DATA LOG10(1),LOG10(2) / O037632, O020232 /
C      DATA LOG10(3),LOG10(4) / O102373, O147770 /
C
C     MACHINE CONSTANTS FOR THE PRIME 50 SERIES SYSTEMS
C     WTIH 32-BIT INTEGERS AND 64V MODE INSTRUCTIONS,
C     SUPPLIED BY IGOR BRAY.
C
C      DATA SMALL(1),SMALL(2) / :10000000000, :00000100001 /
C      DATA LARGE(1),LARGE(2) / :17777777777, :37777677775 /
C      DATA RIGHT(1),RIGHT(2) / :10000000000, :00000000122 /
C      DATA DIVER(1),DIVER(2) / :10000000000, :00000000123 /
C      DATA LOG10(1),LOG10(2) / :11504046501, :07674600177 /
C
C     MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000
C
C      DATA SMALL(1),SMALL(2) / $00000000,  $00100000 /
C      DATA LARGE(1),LARGE(2) / $FFFFFFFF,  $7FEFFFFF /
C      DATA RIGHT(1),RIGHT(2) / $00000000,  $3CA00000 /
C      DATA DIVER(1),DIVER(2) / $00000000,  $3CB00000 /
C      DATA LOG10(1),LOG10(2) / $509F79FF,  $3FD34413 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
C
C      DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 /
C      DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 /
C      DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 /
C      DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 /
C      DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /
C
C     MACHINE CONSTANTS FOR THE VAX UNIX F77 COMPILER
C
C      DATA SMALL(1),SMALL(2) /        128,           0 /
C      DATA LARGE(1),LARGE(2) /     -32769,          -1 /
C      DATA RIGHT(1),RIGHT(2) /       9344,           0 /
C      DATA DIVER(1),DIVER(2) /       9472,           0 /
C      DATA LOG10(1),LOG10(2) /  546979738,  -805796613 /
C
C     MACHINE CONSTANTS FOR THE VAX-11 WITH
C     FORTRAN IV-PLUS COMPILER
C
C      DATA SMALL(1),SMALL(2) / Z00000080, Z00000000 /
C      DATA LARGE(1),LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
C      DATA RIGHT(1),RIGHT(2) / Z00002480, Z00000000 /
C      DATA DIVER(1),DIVER(2) / Z00002500, Z00000000 /
C      DATA LOG10(1),LOG10(2) / Z209A3F9A, ZCFF884FB /
C
C     MACHINE CONSTANTS FOR VAX/VMS VERSION 2.2
C
C      DATA SMALL(1),SMALL(2) /       '80'X,        '0'X /
C      DATA LARGE(1),LARGE(2) / 'FFFF7FFF'X, 'FFFFFFFF'X /
C      DATA RIGHT(1),RIGHT(2) /     '2480'X,        '0'X /
C      DATA DIVER(1),DIVER(2) /     '2500'X,        '0'X /
C      DATA LOG10(1),LOG10(2) / '209A3F9A'X, 'CFF884FB'X /
C
      IF (I .LT. 1  .OR.  I .GT. 5) GOTO 999
      D1MACH = DMACH(I)
      RETURN
  999 WRITE(I1MACH(2),1999) I
 1999 FORMAT(' D1MACH - I OUT OF BOUNDS',I10)
      STOP
      END
C***********************************************************************
C***********************************************************************
      INTEGER FUNCTION I1MACH(I)
C
C  I/O UNIT NUMBERS.
C
C    I1MACH( 1) = THE STANDARD INPUT UNIT.
C
C    I1MACH( 2) = THE STANDARD OUTPUT UNIT.
C
C    I1MACH( 3) = THE STANDARD PUNCH UNIT.
C
C    I1MACH( 4) = THE STANDARD ERROR MESSAGE UNIT.
C
C  WORDS.
C
C    I1MACH( 5) = THE NUMBER OF BITS PER INTEGER STORAGE UNIT.
C
C    I1MACH( 6) = THE NUMBER OF CHARACTERS PER INTEGER STORAGE UNIT.
C
C  INTEGERS.
C
C    ASSUME INTEGERS ARE REPRESENTED IN THE S-DIGIT, BASE-A FORM
C
C               SIGN ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
C
C               WHERE 0 .LE. X(I) .LT. A FOR I=0,...,S-1.
C
C    I1MACH( 7) = A, THE BASE.
C
C    I1MACH( 8) = S, THE NUMBER OF BASE-A DIGITS.
C
C    I1MACH( 9) = A**S - 1, THE LARGEST MAGNITUDE.
C
C  FLOATING-POINT NUMBERS.
C
C    ASSUME FLOATING-POINT NUMBERS ARE REPRESENTED IN THE T-DIGIT,
C    BASE-B FORM
C
C               SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
C
C               WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T,
C               0 .LT. X(1), AND EMIN .LE. E .LE. EMAX.
C
C    I1MACH(10) = B, THE BASE.
C
C  SINGLE-PRECISION
C
C    I1MACH(11) = T, THE NUMBER OF BASE-B DIGITS.
C
C    I1MACH(12) = EMIN, THE SMALLEST EXPONENT E.
C
C    I1MACH(13) = EMAX, THE LARGEST EXPONENT E.
C
C  DOUBLE-PRECISION
C
C    I1MACH(14) = T, THE NUMBER OF BASE-B DIGITS.
C
C    I1MACH(15) = EMIN, THE SMALLEST EXPONENT E.
C
C    I1MACH(16) = EMAX, THE LARGEST EXPONENT E.
C
C  TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT,
C  THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY
C  REMOVING THE C FROM COLUMN 1.  ALSO, THE VALUES OF
C  I1MACH(1) - I1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY
C  WITH THE LOCAL OPERATING SYSTEM.  FOR FORTRAN 77, YOU MAY WISH
C  TO ADJUST THE DATA STATEMENT SO IMACH(6) IS SET TO 1, AND
C  THEN TO COMMENT OUT THE EXECUTABLE TEST ON I .EQ. 6 BELOW.
C  ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED.
C  (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.)
C
C  FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), THE FIRST
C  SET OF CONSTANTS BELOW SHOULD BE APPROPRIATE, EXCEPT PERHAPS
C  FOR IMACH(1) - IMACH(4).
C
      INTEGER IMACH(16),OUTPUT
C
      EQUIVALENCE (IMACH(4),OUTPUT)
C
C     MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T
C     3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T
C     PC 7300), AND 8087 BASED MICROS (E.G. IBM PC AND AT&T 6300).
C
       DATA IMACH( 1) /    5 /
       DATA IMACH( 2) /    6 /
       DATA IMACH( 3) /    7 /
       DATA IMACH( 4) /    6 /
       DATA IMACH( 5) /   32 /
       DATA IMACH( 6) /    4 /
       DATA IMACH( 7) /    2 /
       DATA IMACH( 8) /   31 /
       DATA IMACH( 9) / 2147483647 /
       DATA IMACH(10) /    2 /
       DATA IMACH(11) /   24 /
       DATA IMACH(12) / -125 /
       DATA IMACH(13) /  128 /
       DATA IMACH(14) /   53 /
       DATA IMACH(15) / -1021 /
       DATA IMACH(16) /  1024 /
C
C     MACHINE CONSTANTS FOR AMDAHL MACHINES.
C
C      DATA IMACH( 1) /   5 /
C      DATA IMACH( 2) /   6 /
C      DATA IMACH( 3) /   7 /
C      DATA IMACH( 4) /   6 /
C      DATA IMACH( 5) /  32 /
C      DATA IMACH( 6) /   4 /
C      DATA IMACH( 7) /   2 /
C      DATA IMACH( 8) /  31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /  16 /
C      DATA IMACH(11) /   6 /
C      DATA IMACH(12) / -64 /
C      DATA IMACH(13) /  63 /
C      DATA IMACH(14) /  14 /
C      DATA IMACH(15) / -64 /
C      DATA IMACH(16) /  63 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C      DATA IMACH( 1) /    7 /
C      DATA IMACH( 2) /    2 /
C      DATA IMACH( 3) /    2 /
C      DATA IMACH( 4) /    2 /
C      DATA IMACH( 5) /   36 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   33 /
C      DATA IMACH( 9) / Z1FFFFFFFF /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   24 /
C      DATA IMACH(12) / -256 /
C      DATA IMACH(13) /  255 /
C      DATA IMACH(14) /   60 /
C      DATA IMACH(15) / -256 /
C      DATA IMACH(16) /  255 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
C
C      DATA IMACH( 1) /   5 /
C      DATA IMACH( 2) /   6 /
C      DATA IMACH( 3) /   7 /
C      DATA IMACH( 4) /   6 /
C      DATA IMACH( 5) /  48 /
C      DATA IMACH( 6) /   6 /
C      DATA IMACH( 7) /   2 /
C      DATA IMACH( 8) /  39 /
C      DATA IMACH( 9) / O0007777777777777 /
C      DATA IMACH(10) /   8 /
C      DATA IMACH(11) /  13 /
C      DATA IMACH(12) / -50 /
C      DATA IMACH(13) /  76 /
C      DATA IMACH(14) /  26 /
C      DATA IMACH(15) / -50 /
C      DATA IMACH(16) /  76 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
C
C      DATA IMACH( 1) /   5 /
C      DATA IMACH( 2) /   6 /
C      DATA IMACH( 3) /   7 /
C      DATA IMACH( 4) /   6 /
C      DATA IMACH( 5) /  48 /
C      DATA IMACH( 6) /   6 /
C      DATA IMACH( 7) /   2 /
C      DATA IMACH( 8) /  39 /
C      DATA IMACH( 9) / O0007777777777777 /
C      DATA IMACH(10) /   8 /
C      DATA IMACH(11) /  13 /
C      DATA IMACH(12) / -50 /
C      DATA IMACH(13) /  76 /
C      DATA IMACH(14) /  26 /
C      DATA IMACH(15) / -32754 /
C      DATA IMACH(16) /  32780 /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   60 /
C      DATA IMACH( 6) /   10 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   48 /
C      DATA IMACH( 9) / 00007777777777777777B /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   48 /
C      DATA IMACH(12) / -974 /
C      DATA IMACH(13) / 1070 /
C      DATA IMACH(14) /   96 /
C      DATA IMACH(15) / -927 /
C      DATA IMACH(16) / 1070 /
C
C     MACHINE CONSTANTS FOR CONVEX C-1.
C
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   32 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   24 /
C      DATA IMACH(12) / -128 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   53 /
C      DATA IMACH(15) /-1024 /
C      DATA IMACH(16) / 1023 /
C
C     MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3.
C
C      DATA IMACH( 1) /     5 /
C      DATA IMACH( 2) /     6 /
C      DATA IMACH( 3) /   102 /
C      DATA IMACH( 4) /     6 /
C      DATA IMACH( 5) /    64 /
C      DATA IMACH( 6) /     8 /
C      DATA IMACH( 7) /     2 /
C      DATA IMACH( 8) /    63 /
C      DATA IMACH( 9) /  777777777777777777777B /
C      DATA IMACH(10) /     2 /
C      DATA IMACH(11) /    47 /
C      DATA IMACH(12) / -8189 /
C      DATA IMACH(13) /  8190 /
C      DATA IMACH(14) /    94 /
C      DATA IMACH(15) / -8099 /
C      DATA IMACH(16) /  8190 /
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200.
C
C      DATA IMACH( 1) /   11 /
C      DATA IMACH( 2) /   12 /
C      DATA IMACH( 3) /    8 /
C      DATA IMACH( 4) /   10 /
C      DATA IMACH( 5) /   16 /
C      DATA IMACH( 6) /    2 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   15 /
C      DATA IMACH( 9) /32767 /
C      DATA IMACH(10) /   16 /
C      DATA IMACH(11) /    6 /
C      DATA IMACH(12) /  -64 /
C      DATA IMACH(13) /   63 /
C      DATA IMACH(14) /   14 /
C      DATA IMACH(15) /  -64 /
C      DATA IMACH(16) /   63 /
C
C     MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7.
C
C      DATA IMACH( 1) /       5 /
C      DATA IMACH( 2) /       6 /
C      DATA IMACH( 3) /       0 /
C      DATA IMACH( 4) /       6 /
C      DATA IMACH( 5) /      24 /
C      DATA IMACH( 6) /       3 /
C      DATA IMACH( 7) /       2 /
C      DATA IMACH( 8) /      23 /
C      DATA IMACH( 9) / 8388607 /
C      DATA IMACH(10) /       2 /
C      DATA IMACH(11) /      23 /
C      DATA IMACH(12) /    -127 /
C      DATA IMACH(13) /     127 /
C      DATA IMACH(14) /      38 /
C      DATA IMACH(15) /    -127 /
C      DATA IMACH(16) /     127 /
C
C     MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES.
C
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /   43 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   36 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   35 /
C      DATA IMACH( 9) / O377777777777 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   27 /
C      DATA IMACH(12) / -127 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   63 /
C      DATA IMACH(15) / -127 /
C      DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C     THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86.
C
C      DATA IMACH( 1) /   5 /
C      DATA IMACH( 2) /   6 /
C      DATA IMACH( 3) /   7 /
C      DATA IMACH( 4) /   6 /
C      DATA IMACH( 5) /  32 /
C      DATA IMACH( 6) /   4 /
C      DATA IMACH( 7) /   2 /
C      DATA IMACH( 8) /  31 /
C      DATA IMACH( 9) / Z7FFFFFFF /
C      DATA IMACH(10) /  16 /
C      DATA IMACH(11) /   6 /
C      DATA IMACH(12) / -64 /
C      DATA IMACH(13) /  63 /
C      DATA IMACH(14) /  14 /
C      DATA IMACH(15) / -64 /
C      DATA IMACH(16) /  63 /
C
C     MACHINE CONSTANTS FOR THE INTERDATA 8/32
C     WITH THE UNIX SYSTEM FORTRAN 77 COMPILER.
C
C     FOR THE INTERDATA FORTRAN VII COMPILER REPLACE
C     THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S.
C
C      DATA IMACH( 1) /   5 /
C      DATA IMACH( 2) /   6 /
C      DATA IMACH( 3) /   6 /
C      DATA IMACH( 4) /   6 /
C      DATA IMACH( 5) /  32 /
C      DATA IMACH( 6) /   4 /
C      DATA IMACH( 7) /   2 /
C      DATA IMACH( 8) /  31 /
C      DATA IMACH( 9) / Z'7FFFFFFF' /
C      DATA IMACH(10) /  16 /
C      DATA IMACH(11) /   6 /
C      DATA IMACH(12) / -64 /
C      DATA IMACH(13) /  62 /
C      DATA IMACH(14) /  14 /
C      DATA IMACH(15) / -64 /
C      DATA IMACH(16) /  62 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
C
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   36 /
C      DATA IMACH( 6) /    5 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   35 /
C      DATA IMACH( 9) / "377777777777 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   27 /
C      DATA IMACH(12) / -128 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   54 /
C      DATA IMACH(15) / -101 /
C      DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
C
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   36 /
C      DATA IMACH( 6) /    5 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   35 /
C      DATA IMACH( 9) / "377777777777 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   27 /
C      DATA IMACH(12) / -128 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   62 /
C      DATA IMACH(15) / -128 /
C      DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
C     32-BIT INTEGER ARITHMETIC.
C
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   32 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   24 /
C      DATA IMACH(12) / -127 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   56 /
C      DATA IMACH(15) / -127 /
C      DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
C     16-BIT INTEGER ARITHMETIC.
C
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   16 /
C      DATA IMACH( 6) /    2 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   15 /
C      DATA IMACH( 9) / 32767 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   24 /
C      DATA IMACH(12) / -127 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   56 /
C      DATA IMACH(15) / -127 /
C      DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR THE PRIME 50 SERIES SYSTEMS
C     WTIH 32-BIT INTEGERS AND 64V MODE INSTRUCTIONS,
C     SUPPLIED BY IGOR BRAY.
C
C      DATA IMACH( 1) /            1 /
C      DATA IMACH( 2) /            1 /
C      DATA IMACH( 3) /            2 /
C      DATA IMACH( 4) /            1 /
C      DATA IMACH( 5) /           32 /
C      DATA IMACH( 6) /            4 /
C      DATA IMACH( 7) /            2 /
C      DATA IMACH( 8) /           31 /
C      DATA IMACH( 9) / :17777777777 /
C      DATA IMACH(10) /            2 /
C      DATA IMACH(11) /           23 /
C      DATA IMACH(12) /         -127 /
C      DATA IMACH(13) /         +127 /
C      DATA IMACH(14) /           47 /
C      DATA IMACH(15) /       -32895 /
C      DATA IMACH(16) /       +32637 /
C
C     MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000.
C
C      DATA IMACH( 1) /     0 /
C      DATA IMACH( 2) /     0 /
C      DATA IMACH( 3) /     7 /
C      DATA IMACH( 4) /     0 /
C      DATA IMACH( 5) /    32 /
C      DATA IMACH( 6) /     1 /
C      DATA IMACH( 7) /     2 /
C      DATA IMACH( 8) /    31 /
C      DATA IMACH( 9) /  2147483647 /
C      DATA IMACH(10) /     2 /
C      DATA IMACH(11) /    24 /
C      DATA IMACH(12) /  -125 /
C      DATA IMACH(13) /   128 /
C      DATA IMACH(14) /    53 /
C      DATA IMACH(15) / -1021 /
C      DATA IMACH(16) /  1024 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
C
C     NOTE THAT THE PUNCH UNIT, I1MACH(3), HAS BEEN SET TO 7
C     WHICH IS APPROPRIATE FOR THE UNIVAC-FOR SYSTEM.
C     IF YOU HAVE THE UNIVAC-FTN SYSTEM, SET IT TO 1.
C
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   36 /
C      DATA IMACH( 6) /    6 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   35 /
C      DATA IMACH( 9) / O377777777777 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   27 /
C      DATA IMACH(12) / -128 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   60 /
C      DATA IMACH(15) /-1024 /
C      DATA IMACH(16) / 1023 /
C
C     MACHINE CONSTANTS FOR VAX.
C
C      DATA IMACH( 1) /    5 /
C      DATA IMACH( 2) /    6 /
C      DATA IMACH( 3) /    7 /
C      DATA IMACH( 4) /    6 /
C      DATA IMACH( 5) /   32 /
C      DATA IMACH( 6) /    4 /
C      DATA IMACH( 7) /    2 /
C      DATA IMACH( 8) /   31 /
C      DATA IMACH( 9) / 2147483647 /
C      DATA IMACH(10) /    2 /
C      DATA IMACH(11) /   24 /
C      DATA IMACH(12) / -127 /
C      DATA IMACH(13) /  127 /
C      DATA IMACH(14) /   56 /
C      DATA IMACH(15) / -127 /
C      DATA IMACH(16) /  127 /
      IF (I .LT. 1  .OR.  I .GT. 16) GO TO 999
      I1MACH=IMACH(I)
C/6S
C/7S
      IF(I.EQ.6) I1MACH=1
C/
      RETURN
  999 WRITE(OUTPUT,1999) I
 1999 FORMAT(' I1MACH - I OUT OF BOUNDS',I10)
      STOP
      END

C*** basicops.f
      SUBROUTINE ADD (A,B,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C***********************************************************************
C
C Called by --
C
C  any routine requiring interval addition.
C
C***********************************************************************
C
C Function --
C
C  This routine adds the interval A and the interval B.  It
C  simulates directed roundings with the routine RNDOUT;  the interval
C  result should contain the interval which would have been obtained
C  with exact interval arithmetic.  However, in general it will not
C  be the smallest possible machine-representable such containing
C  interval.  See the documentation in subroutine RNDOUT for more
C  detailed information.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), B(2), RESULT(2)
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C  A        is the first operand to the addition.
C           (INPUT)
C
C  B        is the second operand to the addition.
C           (INPUT)
C
C  RESULT   is the interval-arithmetic sum of A and B.
C           (OUTPUT)
C
C***********************************************************************
C
C Internal variable declarations --
C
      LOGICAL RNDDWN, RNDUP
C
C***********************************************************************
C
C Internal variable descriptions --
C
C RNDDWN    is set to .TRUE. if RNDOUT has to round down, and is set
C           to .FALSE. otherwise.
C
C RNDUP     is set to .TRUE. if RNDOUT has to round up, and is set
C           to .FALSE. otherwise.
C
C***********************************************************************
C
C Common block declarations -- none
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines -- none
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C RNDOUT
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C I/O functions --  none
C
C***********************************************************************
C
C Internal constant declarations -- none
C
C***********************************************************************
C
C Beginning of executable statements --
C
      RNDDWN = (A(1).NE.0D0).AND.(B(1).NE.0D0)
      RNDUP = (A(2).NE.0D0).AND.(B(2).NE.0D0)
C
      RESULT(1) = A(1) + B(1)
      RESULT(2) = A(2) + B(2)
C
      CALL RNDOUT(RESULT,RNDDWN,RNDUP)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE CANCEL (A,B,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   March 13, 1990
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C***********************************************************************
C
C Called by --
C
C  Any routine requiring an interval addend to be canceled (removed)
C  from a previously accumulated interval sum.
C
C***********************************************************************
C
C Function --
C
C  Given an interval B, and a previously accumulated interval sum A
C  for which B was an addend, this routine returns an interval which
C  contains, and is hopefully close to, the sum of the other addends
C  for A.  Directed roundings are simulated with the routine RNDOUT;
C  the interval result should contain the interval which would have
C  been obtained with exact interval arithmetic.  However, in general
C  it will not be the smallest possible machine-representable such
C  containing interval.  See the documentation in subroutine RNDOUT for
C  more detailed information.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), B(2), RESULT(2)
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C  A        is the first operand to the subtraction.
C           (INPUT)
C
C  B        is the second operand to the subtraction
C           (INPUT)
C
C  RESULT   is the interval-arithmetic value of A - B.
C           (OUTPUT)
C
C***********************************************************************
C
C Internal variable declarations --
C
      LOGICAL RNDDWN, RNDUP
C
C***********************************************************************
C
C Internal variable descriptions --
C
C
C RNDDWN    is set to .TRUE. if RNDOUT has to round down, and is set
C           to .FALSE. otherwise.
C
C RNDUP     is set to .TRUE. if RNDOUT has to round up, and is set
C           to .FALSE. otherwise.
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C RNDOUT
C
C***********************************************************************
C
C Beginning of executable statements --
C
      RNDDWN = (B(1).NE.0D0)
      RNDUP = (B(2).NE.0D0)
C
      RESULT(1) = A(1) - B(1)
      RESULT(2) = A(2) - B(2)
C
      CALL RNDOUT(RESULT,RNDDWN,RNDUP)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE IDIV (AA,BB,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   May 21, 1992
C
C Part of the interval arithmetic elementary function package.
C
C***********************************************************************
C
C Function --
C
C  This routine performs interval division.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION AA(2), BB(2), RESULT(2)
C
C***********************************************************************
C
C Internal variable declarations --

      DOUBLE PRECISION C(2)
C
C***********************************************************************
C
C  Common block declarations --
C
      DOUBLE PRECISION MAXX, MXLGM1, NEGINF, POSINF, A(2), PI(2), E(2),
     *   PI2(2), PI3(2), PI4(2), PI6(2), PI8(2), E14(2), A3(2), ONE(2),
     *   OD2F(2), OD3F(2), OD4F(2), OD5F(2), OD6F(2), OD7F(2), OD8F(2),
     *   OD9F(2), OD10F(2), OD11F(2), OD12F(2), OD14F(2), EIGHT(2),
     *   ZERO(2), TWO(2), THREE(2), FOUR(2), SXTNTH(2), NINE(2), TEN(2),
     *   TWOT7(2), SQT10(2), ESXTNT(2), SXTEEN(2), THIRD(2), FOURTH(2),
     *   FIFTH(2), SIXTH(2), SEVNTH(2), EIGHTH(2), NINTH(2), TENTH(2),
     *   ELEVTH(2), TWLVTH(2), THRTTH(2), FORTTH(2), TT7TH(2), TINY2,
     *   CBTEP, SQT3(2), ODSQT3(2)
C
      COMMON /MTHCNS/  MAXX, MXLGM1, NEGINF, POSINF, A, PI, E, PI2, PI3,
     *   PI4, PI6, PI8, ONE, OD2F, OD3F, E14, A3, OD4F, OD5F, OD6F,
     *   OD7F, OD8F, OD9F, OD10F, OD11F, OD12F, OD14F, EIGHT, ZERO, TWO,
     *   THREE, FOUR, SXTNTH, NINE, TEN, TWOT7, SQT10, ESXTNT, SXTEEN,
     *   THIRD, FOURTH, FIFTH, SIXTH, SEVNTH, EIGHTH, NINTH, TENTH,
     *   ELEVTH, TWLVTH, THRTTH, FORTTH, TT7TH, TINY2, CBTEP, SQT3,
     *   ODSQT3
C
      INTEGER ITINY2, JTINY2
C
      COMMON /IMATH/   ITINY2, JTINY2
C
C
C  The above common blocks hold mathematical constants which are used in
C  the elementary function routines.
C
C  Variable descriptions
C
C  MAXX        a double precision representation of the largest
C              representable integer
C
C  MXLGM1      an approximation to the logarithm of .25 times the
C              largest representable machine number, minus 1. This
C              should be a rigorous lower bound on the logarithm of the
C              largest representable machine number.  Its default
C              computation in SIMINI is to use the Fortran LOG function.
C              This should be changed if the Fortran LOG function is not
C              sufficiently accurate.
C
C
C  A           an interval enclosure for 1/pi
C
C  PI          an interval enclosure for pi
C
C  PI2         an interval enclosure for pi/2
C
C  PI3         an interval enclosure for pi/3
C
C  PI4         an interval enclosure for pi/4
C
C  PI6         an interval enclosure for pi/6
C
C  PI8         an interval enclosure for pi/8
C
C  E           an interval enclosure for E
C
C  E14         an interval enclosure for e^{1/4}
C
C  ESXTNT      an interval enclosure for e^{1/16}
C
C  SQT10       an interval enclosure for SQRT(10)
C
C  TINY2       the maximum of the smallest machine number and the
C              reciprocal of the largest machine number.  This quantity
C              is checked to avoid overflow in certain places.
C
C  CBTEP       an approximation to the cube root of six times the
C              largest distance between numbers.
C
C  SQT3        an interval enclosure for the square root of 3.
C
C  ODSQT3      the reciprocal of the square root of 3 times
C              1+ 100*MXULP, used in argument reduction in the
C              arctangent routine.
C
C  ITINY2      the logarithm base 10 of TINY2 truncated to an integer.
C
C  JTINY2      sixteen times the logarithm base e of TINY2 truncated to
C              an integer.
C
C  See the statements in the routine SIMINI for the definitions of
C  the other constants.
C
      INTEGER ISIG, IERR, IROUT, IPRTCL, ISEVER, IERPUN
      COMMON /ERFLGS/ ISIG, IERR, IROUT, IPRTCL, ISEVER, IERPUN
C
C The above common block stores signalling information for error
C conditions.  In particular,
C
C ISIG   is set to 0 at the beginning of SIMINI, and is reset
C        to the severity level of the error (1, 2, or 3) when an error
C        condition occurs.  The user may reset the flag to zero after an
C        error condition, depending on the error, if ISEVER is set to
C        allow execution after an error of its severity.
C
C IERR   is the number of the error condition, if an error occurred in
C        the last routine in INTLIB with error checking.  If no errors
C        occurred in the last such routine called, then IERR is zero.
C        The specific error conditions associated with particular error
C        numbers is defined in the routine ERRTST.
C
C IROUT  is the code number of the package routine in which the error
C        occurred
C
C IPRTCL controls the level of error which is printed.  IPRCTL=0 prints
C        all levels, while IPRCTL=1 prints only errors of level 1 or
C        greater.  If IPRTCL=4, then no error information is printed.
C
C ISEVER gives the error level which will stop execution.  ISEVER=0
C        causes any error to stop execution, while ISEVER>=3 causes only
C        errors of severity 3 to stop execution.  Generally, severity 3
C        errors correspond to when it is impossible to assign a result
C        with any meaningful interpretation.
C
C IERPUN is the Fortran unit number to which errors should be printed.
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  ERRTST, MULT, RNDOUT
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Identifying code for this routine --
C
      IROUT = 3
      IERR = 0
C
C  Do usual interval division if zero is not in the denominator --
C
      IF (BB(1).GT.ZERO(2)) THEN
         C(1) = ONE(1)/BB(2)
         C(2) = ONE(2)/BB(1)
         CALL RNDOUT(C,.TRUE.,.TRUE.)
         CALL MULT(AA,C,RESULT)
      ELSE IF (BB(2).LT.ZERO(1)) THEN
         C(1) = ONE(2)/BB(2)
         C(2) = ONE(1)/BB(1)
         CALL RNDOUT(C,.TRUE.,.TRUE.)
         CALL MULT(AA,C,RESULT)
      ELSE
         IERR = 6
         CALL ERRTST(BB)
         RESULT(1) = NEGINF
         RESULT(2) = POSINF
      END IF
C
      RETURN
C
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE MULT (A, B, RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   November 17, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  MULT multiplies the interval A and the interval B.  It
C  puts the result into output parameter RESULT.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), B(2), RESULT(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TEMP, A1, B1
      DOUBLE PRECISION AA(2), BB(2)
C
C***********************************************************************
C
C Internal Constant Declarations --
C
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO = 0.0D0)
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Pictures for cases.
C
      IF ((ZERO .LE. A(1)) .AND. (ZERO .LE. B(1))) THEN
C  Case 1   ---------------- 0 -----------------
C        A:                    [==========]
C        B:                         [===========]
          RESULT(1) = A(1) * B(1)
          RESULT(2) = A(2) * B(2)
      ELSE IF ((A(1) .LT. ZERO) .AND. (ZERO .LT. A(2))
     1                          .AND. (ZERO .LE. B(1))) THEN
C  Case 2   ---------------- 0 -----------------
C        A:            [==================]
C        B:                         [===========]
          RESULT(1) = A(1) * B(2)
          RESULT(2) = A(2) * B(2)
      ELSE IF ((A(2) .LE. ZERO) .AND. (ZERO .LE. B(1))) THEN
C  Case 3   ---------------- 0 -----------------
C        A:    [==========]
C        B:                         [===========]
          B1 = B(1)
          RESULT(1) = A(1) * B(2)
          RESULT(2) = A(2) * B1
      ELSE IF ((ZERO .LE. A(1)) .AND. (B(1) .LT. ZERO)
     1                          .AND. (ZERO .LT. B(2))) THEN
C  Case 4   ---------------- 0 -----------------
C        A:                     [==========]
C        B:               [===========]
          RESULT(1) = A(2) * B(1)
          RESULT(2) = A(2) * B(2)
      ELSE IF ((A(2) .LE. ZERO) .AND. (B(1) .LT. ZERO)
     1                          .AND. (ZERO .LT. B(2))) THEN
C  Case 5   ---------------- 0 -----------------
C        A:    [==========]
C        B           [===========]
          A1 = A(1)
          B1 = B(1)
          RESULT(1) = A(1) * B(2)
          RESULT(2) = A1 * B1
      ELSE IF ((ZERO .LE. A(1)) .AND. (B(2) .LE. ZERO)) THEN
C  Case 6   ---------------- 0 -----------------
C        A:                    [==========]
C        B:  [===========]
          A1 = A(1)
          RESULT(1) = A(2) * B(1)
          RESULT(2) = A1 * B(2)
      ELSE IF ((A(1) .LT. ZERO) .AND. (ZERO .LT. A(2))
     1                          .AND. (B(2) .LE. ZERO)) THEN
C  Case 7   ---------------- 0 -----------------
C        A:              [==========]
C        B:  [===========]
          A1 = A(1)
          B1 = B(1)
          RESULT(1) = A(2) * B(1)
          RESULT(2) = A1 * B1
      ELSE IF ((A(2) .LE. ZERO) .AND. (B(2) .LE. ZERO)) THEN
C  Case 8   ---------------- 0 -----------------
C        A:    [==========]
C        B:   [===========]
          A1 = A(1)
          B1 = B(1)
          RESULT(1) = A(2) * B(2)
          RESULT(2) = A1 * B1
      ELSE
C  Case 9   ---------------- 0 -----------------
C        A:              [==========]
C        B:              [===========]
C  Must check two cases.
          AA(1) = A(1)
          AA(2) = A(2)
          BB(1) = B(1)
          BB(2) = B(2)
          RESULT(1) = AA(1) * BB(2)
          TEMP = AA(2) * BB(1)
          IF (TEMP .LT. RESULT(1)) THEN
              RESULT(1) = TEMP
          ELSE
          END IF
              RESULT(2) = AA(1) * BB(1)
          TEMP = AA(2) * BB(2)
          IF (TEMP .GT. RESULT(2)) THEN
              RESULT(2) = TEMP
          ELSE
          END IF
      END IF
C
      CALL RNDOUT(RESULT,.TRUE.,.TRUE.)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE RNDOUT (X,RNDDWN,RNDUP)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C   Modified August, 1991 and April, 1992 by Kaisheng Du and
C   R. Baker Kearfott to be more accurate near the underflow
C   threshold.
C
C Part of the generalized bisection package
C (interval arithmetic subpackage)
C
C***********************************************************************
C
C Function --
C
C  This routine is intended to simulate directed roundings in a
C  reasonably transportable way.  It is called for each elementary
C  operation involving intervals.  The endpoints of the result interval
C  are first computed with the machine's usual floating point
C  arithmetic.
C
C  If RNDDWN = .TRUE., then this routine decreases the left
C  endpoint of that approximate result by the absolute value of
C  that endpoint times a rigorous estimate for the maximum relative
C  error in an elementary operation.
C
C  If RNDUP = .TRUE., then this routine increases the right
C  endpoint of that approximate result by the absolute value of
C  that endpoint times a rigorous estimate for the maximum relative
C  error in an elementary operation.
C
C  For this routine to work properly, a machine-dependent parameter
C  must be installed in the routine SIMINI.  See the documentation in
C  that routine for details.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(2)
      LOGICAL RNDDWN, RNDUP
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C X         is the interval to be adjusted.
C           (I/O)
C
C RNDDWN    is set to .TRUE. if the left endpoint is to be adjusted, and
C           is set to .FALSE. otherwise.
C           (INPUT)
C
C RNDUP     is set to .TRUE. if the right endpoint is to be adjusted,
C           and is set to .FALSE. otherwise.
C           (INPUT)
C
C***********************************************************************
C
C Common block declarations --
C
      DOUBLE PRECISION MXULP, TTINY2, TOL0
      COMMON /MACH1/ MXULP, TTINY2, TOL0
C
C  This common block holds machine parameters which are set in
C  SIMINI and used here.
C
C  Variable descriptions
C
C  MXULP       (machine epsilon)
C                * (maximum error in ULP's of the floating pt. op's)
C
C  TTINY2      2 * (smallest representable positive machine number)
C                * (maximum error in ULP's of the floating pt. op's)
C
C  TOL0        TTINY2 / MXULP
C
      DOUBLE PRECISION TINY, TEST
      COMMON /MACH2/ TINY, TEST
C
C  See SIMINI for an explanation of the common block MACH2.
C
C***********************************************************************
C
C Beginning of executable statements --
C
      IF (RNDDWN) THEN
         IF (X(1).GE.TEST) THEN
            X(1) = (1.D0 - MXULP ) * X(1)
         ELSE IF (X(1).LE.-TEST) THEN
            X(1) = (1D0 + MXULP ) * X(1)
         ELSE IF (X(1).LE.0.D0) THEN
            X(1) = -TEST
         ELSE
            X(1) = 0.D0
         END IF
      END IF
C
      IF (RNDUP) THEN
         IF (X(2).GE.TEST) THEN
            X(2) = (1.D0 + MXULP )* X(2)
         ELSE IF (X(2).LE.-TEST) THEN
            X(2) = (1.D0 - MXULP ) * X(2)
         ELSE IF(X(2).GE.0D0) THEN
            X(2) = TEST
         ELSE
            X(2) = 0.D0
         ENDIF
      END IF
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SCLADD (R,A,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C***********************************************************************
C
C Called by --
C
C  Any routine requiring addition of a point value to an interval.
C
C***********************************************************************
C
C FUNCTION --
C
C  This routine adds the interval A to the point R.  It simulates
C  directed roundings with the routine RNDOUT;  the interval
C  result should contain the interval which would have been obtained
C  with exact interval arithmetic.  However, in general it will not
C  be the smallest possible machine-representable such containing
C  interval.  See the documentation in subroutine RNDOUT for more
C  detailed information.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION R, A(2), RESULT(2)
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C  R        is the point to be added to the interval.
C           (INPUT)
C
C  A        is the interval to be added to the point.
C           (INPUT)
C
C  RESULT   is the interval-arithmetic sum of R and B.
C           (OUTPUT)
C
C***********************************************************************
C
C Internal variable declarations --
C
      LOGICAL RNDDWN, RNDUP
C
C***********************************************************************
C
C Internal variable descriptions --
C
C RNDDWN    is set to .TRUE. if RNDOUT has to round down, and is set
C           to .FALSE. otherwise.
C
C RNDUP     is set to .TRUE. if RNDOUT has to round up, and is set
C           to .FALSE. otherwise.
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  RNDOUT
C
C***********************************************************************
C
C Beginning of executable statements --
C
      RNDDWN = (R.NE.0D0).AND.(A(1).NE.0D0)
      RNDUP = (R.NE.0D0).AND.(A(2).NE.0D0)
C
      RESULT(1) = R + A(1)
      RESULT(2) = R + A(2)
C
      CALL RNDOUT(RESULT,RNDDWN,RNDUP)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SCLMLT (R,A,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C Modified by Manuel Novoa III on March 13, 1990 to clean the code
C   slightly, and to remove the need for MAX and MIN.
C
C***********************************************************************
C
C Called by --
C
C  Any routine requiring multiplication of an interval and a point
C  value.
C
C***********************************************************************
C
C Function --
C
C  This routine multiplies the interval A and the point R.  It
C  simulates directed roundings with the routine RNDOUT;  the interval
C  result should contain the interval which would have been obtained
C  with exact interval arithmetic.  However, in general it will not
C  be the smallest possible machine-representable such containing
C  interval.  See the documentation in subroutine RNDOUT for more
C  detailed information.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION R, A(2), RESULT(2)
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C  R        is the point to be multiplied to the interval.
C           (INPUT)
C
C  A        is the interval to be multiplied to the point.
C           (INPUT)
C
C  RESULT   is the interval-arithmetic product R * B.
C           (OUTPUT)
C
C***********************************************************************
C
C Internal variable declarations --
C
      LOGICAL RNDDWN, RNDUP
      DOUBLE PRECISION T1, T2
C
C***********************************************************************
C
C Internal variable descriptions --
C
C RNDDWN    is set to .TRUE. if RNDOUT is to round down, and is set
C           to .FALSE. otherwise.
C
C RNDUP     is set to .TRUE. if RNDOUT is to round up, and is set
C           to .FALSE. otherwise.
C
C T1 and T2 are temporary variables.
C
C***********************************************************************
C
C Common block declarations -- none
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines --  none
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  RNDOUT
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C I/O functions --  none
C
C***********************************************************************
C
C Internal constant declarations -- none
C
C***********************************************************************
C
C Beginning of executable statements --
C
      IF ((R.EQ.0D0).OR.((A(1).EQ.0D0).AND.(A(2).EQ.0D0))) THEN
         RESULT(1) = 0D0
         RESULT(2) = 0D0
         RETURN
      END IF
C
      T1 = A(1)
      T2 = A(2)
      RNDDWN = .TRUE.
      RNDUP = .TRUE.
C
      IF (T1.EQ.0D0) THEN
         IF (R.LT.0D0) THEN
            RESULT(1) = R * T2
            RESULT(2) = 0D0
            RNDUP = .FALSE.
         ELSE
            RESULT(1) = 0D0
            RESULT(2) = R * T2
            RNDDWN = .FALSE.
         END IF
      ELSE IF (T2.EQ.0D0) THEN
         IF (R.LT.0D0) THEN
            RESULT(1) = 0D0
            RESULT(2) = R * T1
            RNDDWN = .FALSE.
         ELSE
            RESULT(1) = R * T1
            RESULT(2) = 0D0
            RNDUP = .FALSE.
         END IF
      ELSE IF (R.GT.0D0) THEN
         RESULT(1) = R * T1
         RESULT(2) = R * T2
      ELSE
         RESULT(1) = R * T2
         RESULT(2) = R * T1
      END IF
C
      CALL RNDOUT(RESULT,RNDDWN,RNDUP)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SUB (A,B,RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C   R. Baker Kearfott
C
C         and
C
C   Manuel Novoa III
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   September 29, 1987
C
C Part of the generalized bisection package
C (interval arithmetic subpackage).
C
C***********************************************************************
C
C Called by --
C
C  Any routine requiring interval subtraction.
C
C***********************************************************************
C
C Function --
C
C  This routine subtracts the interval B from the interval A.  It
C  simulates directed roundings with the routine RNDOUT;  the interval
C  result should contain the interval which would have been obtained
C  with exact interval arithmetic.  However, in general it will not
C  be the smallest possible machine-representable such containing
C  interval.  See the documentation in subroutine RNDOUT for more
C  detailed information.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), B(2), RESULT(2)
C
C***********************************************************************
C
C Argument descriptions -- (INPUT  = set on entry and not alterable)
C                          (OUTPUT = to be set by the routine)
C                          (I/O    = set on entry but alterable)
C
C  A        is the first operand to the subtraction.
C           (INPUT)
C
C  B        is the second operand to the subtraction
C           (INPUT)
C
C  RESULT   is the interval-arithmetic value of A - B.
C           (OUTPUT)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TA1, TA2, TB1, TB2
      LOGICAL RNDDWN, RNDUP
C
C***********************************************************************
C
C Internal variable descriptions --
C
C TA1, TA2, TB1, and TB2 are temporaries.
C
C RNDDWN    is set to .TRUE. if RNDOUT has to round down, and is set
C           to .FALSE. otherwise.
C
C RNDUP     is set to .TRUE. if RNDOUT has to round up, and is set
C           to .FALSE. otherwise.
C
C***********************************************************************
C
C Common block declarations -- none
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines -- none
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C RNDOUT
C
C***********************************************************************
C
C User-supplied functions and subroutines -- none
C
C***********************************************************************
C
C I/O functions --  none
C
C***********************************************************************
C
C Internal constant declarations -- none
C
C***********************************************************************
C
C Beginning of executable statements --
C
      TA1 = A(1)
      TA2 = A(2)
      TB1 = B(1)
      TB2 = B(2)
C
      RNDDWN = (TB2.NE.0D0)
      RNDUP = (TB1.NE.0D0)
C
      RESULT(1) = TA1 - TB2
      RESULT(2) = TA2 - TB1
C
      CALL RNDOUT(RESULT,RNDDWN,RNDUP)
C
      RETURN
      END

c*** utilfuns.f
      SUBROUTINE ICAP (A, B, RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  DICAP computes the intersection of the two intervals A and B, and
C  places  the result in RESULT.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), B(2), RESULT(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
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 Fortran-supplied functions and subroutines --
C
C   MIN, MAX
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C   ERRTST
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION T(2)
C
C***********************************************************************
C
C Beginning of executable statements --
C
C  Identifying code for this routine --
C
      IROUT = 13
      IERR = 0
C
      T(1) = MAX (A(1), B(1))
      T(2) = MIN (A(2), B(2))
      RESULT(1) = T(1)
      RESULT(2) = T(2)
      IF (RESULT(1).GT.RESULT(2)) THEN
         IERR=13
         CALL ERRTST(RESULT)
      END IF
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      LOGICAL FUNCTION IDISJ (A, B)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  This function returns .TRUE. if the intervals A and B are disjoint,
C  and .FALSE. otherwise.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), B(2)
C
C***********************************************************************
C
C Beginning of executable statements --
C
      IDISJ = (A(2) .LT. B(1)) .OR. (B(2) .LT. A(1))
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE IHULL (A, B, RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  IHULL returns the convex-hull of the interval A and the interval B
C  in RESULT.  If one of the intervals is empty (lower bound is greater
C  than upper bound), then just the upper interval is returned.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), B(2), RESULT(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION T(2)
C
C***********************************************************************
C
C Beginning of executable statements --
C
      IF ( A(1).GT.A(2) ) THEN
         IF (B(1).GT.B(2)) THEN
            T(1) = MAX(A(2),B(2))
            T(2) = MIN(A(1),B(2))
         ELSE
            T(1) = B(1)
            T(2) = B(2)
         END IF
      ELSE IF ( B(1).GT.B(2) ) THEN
         T(1) = A(1)
         T(2) = A(2)
      ELSE
         T(1) = MIN (A(1), B(1))
         T(2) = MAX (A(2), B(2))
         RESULT(1) = T(1)
         RESULT(2) = T(2)
      END IF
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      LOGICAL FUNCTION IILEI (A, B)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  IILEI sets its value to .TRUE. if interval A is in the closure of
C  interval B.  The value of IILEI is .FALSE. otherwise.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), B(2)
C
C***********************************************************************
C
C Beginning of executable statements --
C
      IILEI = (A(1) .GE. B(1)) .AND. (A(2) .LE. B(2))
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      LOGICAL FUNCTION IILTI (A, B)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  IILE sets its value to .TRUE. if interval A is in the interior of
C  interval B.  The value of IILE is .FALSE. otherwise.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), B(2)
C
C***********************************************************************
C
C Beginning of executable statements --
C
      IILTI = (A(1) .GT. B(1)) .AND. (A(2) .LT. B(2))
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      DOUBLE PRECISION FUNCTION IINF (A)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  IINF returns the lower endpoint of the interval A.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2)
C
C***********************************************************************
C
C Beginning of executable statements --
C
      IINF = A(1)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      DOUBLE PRECISION FUNCTION IMID (X)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  IMID returns a floating point approximation to the midpoint
C  of the interval A, using available floating point arithmetic.  The
C  value returned by this routine can be considered to DEFINE the
C  midpoint.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION X(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
C***********************************************************************
C
C Beginning of executable statements --
C
      IMID = X(1) + (X(2) - X(1)) / TWO(1)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      DOUBLE PRECISION FUNCTION IMIG (A)
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   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  IMIG returns the mignitude of the interval A. Since ABS is not
C  assumed to give an exact result, the result is rounded down.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TMP(2)
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines --
C
C  ABS, MIN
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  RNDOUT
C
C***********************************************************************
C
C Beginning of executable statements --
C
      IF ( ((A(1).GT.0D0) .AND. (A(2).GT.0D0))
     * .OR.((A(1).LT.0D0) .AND. (A(2).LT.0D0)) ) THEN
         TMP(1) = ABS(A(2))
         CALL RNDOUT(TMP,.TRUE.,.FALSE.)
         TMP(2)=TMP(1)
         TMP(1) = ABS(A(1))
         CALL RNDOUT(TMP,.TRUE.,.FALSE.)
         IMIG = MIN (TMP(1), TMP(2))
      ELSE
         IMIG = 0D0
      END IF
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE INEG (A, RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  INEG performs interval unary negation.  The result is rounded out in
C  case the negatives of the endpoints are not representable.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), RESULT(2)
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  RNDOUT
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION T(2)
C
C***********************************************************************
C
C Beginning of executable statements --
C
      T(1) = A(1)
      T(2) = A(2)
      RESULT(1) = -T(2)
      RESULT(2) = -T(1)
      CALL RNDOUT (RESULT, .TRUE.,.TRUE.)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      DOUBLE PRECISION FUNCTION INTABS (A)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  INTABS returns the maximum absolute value of the interval A.
C  Since ABS is not assumed to give an exact result, the result is
C  rounded up.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TMP(2)
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines --
C
C  ABS, MAX
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  RNDOUT
C
C***********************************************************************
C
C Beginning of executable statements --
C
      TMP(2) = ABS(A(1))
      CALL RNDOUT(TMP,.FALSE.,.TRUE.)
      TMP(1)=TMP(2)
      TMP(2) = ABS(A(2))
      CALL RNDOUT(TMP,.FALSE.,.TRUE.)
C
      INTABS = MAX (TMP(1), TMP(2))
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      LOGICAL FUNCTION IRLEI (A, B)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  IRLEI sets its value to .TRUE. if double precision value A is in the
C  closure of the interval B. The value of IRLEI is set to .FALSE.
C  otherwise.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A, B(2)
C
C***********************************************************************
C
C Beginning of executable statements --
C
      IRLEI = (A .GE. B(1)) .AND. (A .LE. B(2))
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      LOGICAL FUNCTION IRLTI (A, B)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  IRLEI sets its value to .TRUE. if double precision value A is in the
C  interior of the interval B. The value of IRLEI is set to .FALSE.
C  otherwise.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A, B(2)
C
C***********************************************************************
C
C Beginning of executable statements --
C
      IRLTI = (A .GT. B(1)) .AND. (A .LT. B(2))
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      DOUBLE PRECISION FUNCTION ISUP (A)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  IINF returns the upper endpoint of the interval A.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2)
C
C***********************************************************************
C
C Beginning of executable statements --
C
      ISUP = A(2)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE IVL1 (R, RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  IVL1 constructs an interval RESULT from the double precision variable
C  R.  This is done by using simulated directed roundings with RNDOUT.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION R, RESULT(2)
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  RNDOUT
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION T
C
C***********************************************************************
C
C Beginning of executable statements --
C
      T = R
      RESULT(1) = T
      RESULT(2) = T
      IF (T.NE.0D0) CALL RNDOUT (RESULT, .TRUE.,.TRUE.)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE IVL2 (R1, R2, RESULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C Written by:
C
C
c   Hong Jiang
c   Dept. of Math., Stat. and Computer Science
c   Marquette University
C   Milwaukee, WI 53233
C
C         and
C
C   R. Baker Kearfott
C
C   Department of Mathematics
C   U.S.L. Box 4-1010
C   Lafayette, LA 70504
C
C   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  IVL2 constructs an interval RESULT, roughly equal to [R1,R2], from
C  the double precision variables R1 and R2.  This is done by using
C  simulated directed roundings with RNDOUT.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION R1, R2, RESULT(2)
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  RNDOUT
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION T(2)
      LOGICAL RNDUP, RNDDWN
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 Beginning of executable statements --
C
C  Identifying code for this routine --
C
      IROUT = 15
      IERR = 0
C
      T(1) = R1
      T(2) = R2
      RESULT(1) = T(1)
      RESULT(2) = T(2)
C
      IF (RESULT(2).LT.RESULT(1)) THEN
         IERR = 1
         CALL ERRTST(RESULT)
      END IF
C
      RNDDWN = T(1).NE.0D0
      RNDUP = T(2).NE.0D0
      CALL RNDOUT (RESULT, RNDDWN, RNDUP)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE IVLABS (A, 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   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  IVLABS returns a rigorous bound on the range of the absolute value
C  of the interval A.  The intrinsic function ABS is assumed to be
C  accurate to the same accuracy as the elementary operations.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), RESULT(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TMP(2), TA(2)
C
C***********************************************************************
C
C Fortran-supplied functions and subroutines --
C
C  ABS
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  RNDOUT
C
C***********************************************************************
C
C Beginning of executable statements --
C
      TA(1) = A(1)
      TA(2) = A(2)

      TMP(1) = ABS(TA(1))
      TMP(2) = ABS(TA(2))

      IF (TMP(1).LE.TMP(2)) THEN
         RESULT(1) = TMP(1)
         RESULT(2) = TMP(2)
      ELSE
         RESULT(1) = TMP(2)
         RESULT(2) = TMP(1)
      END IF

      CALL RNDOUT(RESULT,.TRUE.,.TRUE.)

      IF ( ( TA(1).LE.0 .AND. TA(2).GE.0 ) .OR. RESULT(1).LT. 0)
     1   RESULT(1) = 0D0

      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE IVLI (A, 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   January 8, 1993
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  IVLI places the contents of interval A in RESULT.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2), RESULT(2)
C
C***********************************************************************
C
C Beginning of executable statements --
C
      RESULT(1) = A(1)
      RESULT(2) = A(2)
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      DOUBLE PRECISION FUNCTION IWID (A)
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   April 11, 1992
C
C Part of the interval elementary function library
C (Utility function)
C
C***********************************************************************
C
C Function --
C
C  IWID returns the width of the interval A, rounded up.
C
C***********************************************************************
C
C Argument declarations --
C
      DOUBLE PRECISION A(2)
C
C***********************************************************************
C
C Internal variable declarations --
C
      DOUBLE PRECISION TMP(2)
C
C***********************************************************************
C
C Package-supplied functions and subroutines --
C
C  RNDOUT
C
C***********************************************************************
C
C Beginning of executable statements --
C
      TMP(2) = A(2)-A(1)
      CALL RNDOUT(TMP,.FALSE.,.TRUE.)
C
      IWID = TMP(2)
C
      RETURN
      END

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