C      ALGORITHM 763, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 22, NO. 4, December, 1996, P.  385--392
C

#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	Doc
#	Libs
#	Src
#	Drivers
# This archive created: Wed Feb 19 19:33:13 1997
export PATH; PATH=/bin:$PATH
if test ! -d 'Doc'
then
	mkdir 'Doc'
fi
cd 'Doc'
if test -f 'makefile'
then
	echo shar: will not over-write existing file "'makefile'"
else
cat << \SHAR_EOF > 'makefile'
f90= frt -Am -Hasu -fi -v9
f90= f90
f90= xlf90 
f90= epcf90 -d20 -C 

INT=f90d1mach.o ivl_def.o intaritb.o testsys.o intlib.o
TEST=f90d1mach.o ivl_def.o intaritb.o intlib.o test_f90_intarith.o
 
HSL12= /rutherford/num-arcu/hsl12/source/
HSL12t= /rutherford/num-arcu/hsl12/tests/

all:	test int
 
test:	$(TEST)
	$(f90)  $(TEST)
	a.out 
	diff TEST_F90_INTARITH.OUT tinysmpl.out

int:	$(INT)
	$(f90)  $(INT)
	a.out 
	diff INTARITH.OUT sample.out

f90d1mach.o:	f90d1mach.f90
	$(f90) -c f90d1mach.f90 

ivl_def.o:	ivl_def.f90
	$(f90) -c ivl_def.f90 

intaritb.o:	intaritb.f90
	$(f90) -c intaritb.f90 

testsys.o:	testsys.f90
	$(f90) -c -o testsys.o testsys.f90 

intlib.o:	intlib.f90
	$(f90) -c intlib.f90 

test_f90_intarith.o:	test_f90_intarith.f90
	$(f90) -c test_f90_intarith.f90 

SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Libs'
then
	mkdir 'Libs'
fi
cd 'Libs'
if test -f 'i1mach.f90'
then
	echo shar: will not over-write existing file "'i1mach.f90'"
else
cat << \SHAR_EOF > 'i1mach.f90'
!***********************************************************************
!***********************************************************************
      INTEGER FUNCTION I1MACH(I)
      INTEGER I
!
!  I/O UNIT NUMBERS.
!
!    I1MACH( 1) = THE STANDARD INPUT UNIT.
!
!    I1MACH( 2) = THE STANDARD OUTPUT UNIT.
!
!    I1MACH( 3) = THE STANDARD PUNCH UNIT.
!
!    I1MACH( 4) = THE STANDARD ERROR MESSAGE UNIT.
!
!  WORDS.
!
!    I1MACH( 5) = THE NUMBER OF BITS PER INTEGER STORAGE UNIT.
!
!    I1MACH( 6) = THE NUMBER OF CHARACTERS PER INTEGER STORAGE UNIT.
!
!  INTEGERS.
!
!    ASSUME INTEGERS ARE REPRESENTED IN THE S-DIGIT, BASE-A FORM
!
!               SIGN ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
!
!               WHERE 0 .LE. X(I) .LT. A FOR I=0,...,S-1.
!
!    I1MACH( 7) = A, THE BASE.
!
!    I1MACH( 8) = S, THE NUMBER OF BASE-A DIGITS.
!
!    I1MACH( 9) = A**S - 1, THE LARGEST MAGNITUDE.
!
!  FLOATING-POINT NUMBERS.
!
!    ASSUME FLOATING-POINT NUMBERS ARE REPRESENTED IN THE T-DIGIT,
!    BASE-B FORM
!
!               SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
!
!               WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T,
!               0 .LT. X(1), AND EMIN .LE. E .LE. EMAX.
!
!    I1MACH(10) = B, THE BASE.
!
!  SINGLE-PRECISION
!
!    I1MACH(11) = T, THE NUMBER OF BASE-B DIGITS.
!
!    I1MACH(12) = EMIN, THE SMALLEST EXPONENT E.
!
!    I1MACH(13) = EMAX, THE LARGEST EXPONENT E.
!
!  DOUBLE-PRECISION
!
!    I1MACH(14) = T, THE NUMBER OF BASE-B DIGITS.
!
!    I1MACH(15) = EMIN, THE SMALLEST EXPONENT E.
!
!    I1MACH(16) = EMAX, THE LARGEST EXPONENT E.
!
!  TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT,
!  THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY
!  REMOVING THE C FROM COLUMN 1.  ALSO, THE VALUES OF
!  I1MACH(1) - I1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY
!  WITH THE LOCAL OPERATING SYSTEM.  FOR FORTRAN 77, YOU MAY WISH
!  TO ADJUST THE DATA STATEMENT SO IMACH(6) IS SET TO 1, AND
!  THEN TO COMMENT OUT THE EXECUTABLE TEST ON I .EQ. 6 BELOW.
!  ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED.
!  (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.)
!
!  FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), THE FIRST
!  SET OF CONSTANTS BELOW SHOULD BE APPROPRIATE, EXCEPT PERHAPS
!  FOR IMACH(1) - IMACH(4).
!
      INTEGER IMACH(16),OUTPUT
!
      EQUIVALENCE (IMACH(4),OUTPUT)
!
!     MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T
!     3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T
!     PC 7300), AND 8087 BASED MICROS (E.G. IBM PC AND AT&T 6300).
!
       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 /
!
!     MACHINE CONSTANTS FOR AMDAHL MACHINES.
!
!      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) /  16 /
!      DATA IMACH(11) /   6 /
!      DATA IMACH(12) / -64 /
!      DATA IMACH(13) /  63 /
!      DATA IMACH(14) /  14 /
!      DATA IMACH(15) / -64 /
!      DATA IMACH(16) /  63 /
!
!     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
!
!      DATA IMACH( 1) /    7 /
!      DATA IMACH( 2) /    2 /
!      DATA IMACH( 3) /    2 /
!      DATA IMACH( 4) /    2 /
!      DATA IMACH( 5) /   36 /
!      DATA IMACH( 6) /    4 /
!      DATA IMACH( 7) /    2 /
!      DATA IMACH( 8) /   33 /
!      DATA IMACH( 9) / Z1FFFFFFFF /
!      DATA IMACH(10) /    2 /
!      DATA IMACH(11) /   24 /
!      DATA IMACH(12) / -256 /
!      DATA IMACH(13) /  255 /
!      DATA IMACH(14) /   60 /
!      DATA IMACH(15) / -256 /
!      DATA IMACH(16) /  255 /
!
!     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
!
!      DATA IMACH( 1) /   5 /
!      DATA IMACH( 2) /   6 /
!      DATA IMACH( 3) /   7 /
!      DATA IMACH( 4) /   6 /
!      DATA IMACH( 5) /  48 /
!      DATA IMACH( 6) /   6 /
!      DATA IMACH( 7) /   2 /
!      DATA IMACH( 8) /  39 /
!      DATA IMACH( 9) / O0007777777777777 /
!      DATA IMACH(10) /   8 /
!      DATA IMACH(11) /  13 /
!      DATA IMACH(12) / -50 /
!      DATA IMACH(13) /  76 /
!      DATA IMACH(14) /  26 /
!      DATA IMACH(15) / -50 /
!      DATA IMACH(16) /  76 /
!
!     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
!
!      DATA IMACH( 1) /   5 /
!      DATA IMACH( 2) /   6 /
!      DATA IMACH( 3) /   7 /
!      DATA IMACH( 4) /   6 /
!      DATA IMACH( 5) /  48 /
!      DATA IMACH( 6) /   6 /
!      DATA IMACH( 7) /   2 /
!      DATA IMACH( 8) /  39 /
!      DATA IMACH( 9) / O0007777777777777 /
!      DATA IMACH(10) /   8 /
!      DATA IMACH(11) /  13 /
!      DATA IMACH(12) / -50 /
!      DATA IMACH(13) /  76 /
!      DATA IMACH(14) /  26 /
!      DATA IMACH(15) / -32754 /
!      DATA IMACH(16) /  32780 /
!
!     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
!
!      DATA IMACH( 1) /    5 /
!      DATA IMACH( 2) /    6 /
!      DATA IMACH( 3) /    7 /
!      DATA IMACH( 4) /    6 /
!      DATA IMACH( 5) /   60 /
!      DATA IMACH( 6) /   10 /
!      DATA IMACH( 7) /    2 /
!      DATA IMACH( 8) /   48 /
!      DATA IMACH( 9) / 00007777777777777777B /
!      DATA IMACH(10) /    2 /
!      DATA IMACH(11) /   48 /
!      DATA IMACH(12) / -974 /
!      DATA IMACH(13) / 1070 /
!      DATA IMACH(14) /   96 /
!      DATA IMACH(15) / -927 /
!      DATA IMACH(16) / 1070 /
!
!     MACHINE CONSTANTS FOR CONVEX C-1.
!
!      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) / -128 /
!      DATA IMACH(13) /  127 /
!      DATA IMACH(14) /   53 /
!      DATA IMACH(15) /-1024 /
!      DATA IMACH(16) / 1023 /
!
!     MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3.
!
!      DATA IMACH( 1) /     5 /
!      DATA IMACH( 2) /     6 /
!      DATA IMACH( 3) /   102 /
!      DATA IMACH( 4) /     6 /
!      DATA IMACH( 5) /    64 /
!      DATA IMACH( 6) /     8 /
!      DATA IMACH( 7) /     2 /
!      DATA IMACH( 8) /    63 /
!      DATA IMACH( 9) /  777777777777777777777B /
!      DATA IMACH(10) /     2 /
!      DATA IMACH(11) /    47 /
!      DATA IMACH(12) / -8189 /
!      DATA IMACH(13) /  8190 /
!      DATA IMACH(14) /    94 /
!      DATA IMACH(15) / -8099 /
!      DATA IMACH(16) /  8190 /
!
!     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200.
!
!      DATA IMACH( 1) /   11 /
!      DATA IMACH( 2) /   12 /
!      DATA IMACH( 3) /    8 /
!      DATA IMACH( 4) /   10 /
!      DATA IMACH( 5) /   16 /
!      DATA IMACH( 6) /    2 /
!      DATA IMACH( 7) /    2 /
!      DATA IMACH( 8) /   15 /
!      DATA IMACH( 9) /32767 /
!      DATA IMACH(10) /   16 /
!      DATA IMACH(11) /    6 /
!      DATA IMACH(12) /  -64 /
!      DATA IMACH(13) /   63 /
!      DATA IMACH(14) /   14 /
!      DATA IMACH(15) /  -64 /
!      DATA IMACH(16) /   63 /
!
!     MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7.
!
!      DATA IMACH( 1) /       5 /
!      DATA IMACH( 2) /       6 /
!      DATA IMACH( 3) /       0 /
!      DATA IMACH( 4) /       6 /
!      DATA IMACH( 5) /      24 /
!      DATA IMACH( 6) /       3 /
!      DATA IMACH( 7) /       2 /
!      DATA IMACH( 8) /      23 /
!      DATA IMACH( 9) / 8388607 /
!      DATA IMACH(10) /       2 /
!      DATA IMACH(11) /      23 /
!      DATA IMACH(12) /    -127 /
!      DATA IMACH(13) /     127 /
!      DATA IMACH(14) /      38 /
!      DATA IMACH(15) /    -127 /
!      DATA IMACH(16) /     127 /
!
!     MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES.
!
!      DATA IMACH( 1) /    5 /
!      DATA IMACH( 2) /    6 /
!      DATA IMACH( 3) /   43 /
!      DATA IMACH( 4) /    6 /
!      DATA IMACH( 5) /   36 /
!      DATA IMACH( 6) /    4 /
!      DATA IMACH( 7) /    2 /
!      DATA IMACH( 8) /   35 /
!      DATA IMACH( 9) / O377777777777 /
!      DATA IMACH(10) /    2 /
!      DATA IMACH(11) /   27 /
!      DATA IMACH(12) / -127 /
!      DATA IMACH(13) /  127 /
!      DATA IMACH(14) /   63 /
!      DATA IMACH(15) / -127 /
!      DATA IMACH(16) /  127 /
!
!     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
!     THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86.
!
!      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(16) /  1024 /
!      DATA IMACH( 8) /  31 /
!      DATA IMACH( 9) / Z7FFFFFFF /
!      DATA IMACH(10) /  16 /
!      DATA IMACH(11) /   6 /
!      DATA IMACH(12) / -64 /
!      DATA IMACH(13) /  63 /
!      DATA IMACH(14) /  14 /
!      DATA IMACH(15) / -64 /
!      DATA IMACH(16) /  63 /
!
!     MACHINE CONSTANTS FOR THE INTERDATA 8/32
!     WITH THE UNIX SYSTEM FORTRAN 77 COMPILER.
!
!     FOR THE INTERDATA FORTRAN VII COMPILER REPLACE
!     THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S.
!
!      DATA IMACH( 1) /   5 /
!      DATA IMACH( 2) /   6 /
!      DATA IMACH( 3) /   6 /
!      DATA IMACH( 4) /   6 /
!      DATA IMACH( 5) /  32 /
!      DATA IMACH( 6) /   4 /
!      DATA IMACH( 7) /   2 /
!      DATA IMACH( 8) /  31 /
!      DATA IMACH( 9) / Z'7FFFFFFF' /
!      DATA IMACH(10) /  16 /
!      DATA IMACH(11) /   6 /
!      DATA IMACH(12) / -64 /
!      DATA IMACH(13) /  62 /
!      DATA IMACH(14) /  14 /
!      DATA IMACH(15) / -64 /
!      DATA IMACH(16) /  62 /
!
!     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
!
!      DATA IMACH( 1) /    5 /
!      DATA IMACH( 2) /    6 /
!      DATA IMACH( 3) /    7 /
!      DATA IMACH( 4) /    6 /
!      DATA IMACH( 5) /   36 /
!      DATA IMACH( 6) /    5 /
!      DATA IMACH( 7) /    2 /
!      DATA IMACH( 8) /   35 /
!      DATA IMACH( 9) / "377777777777 /
!      DATA IMACH(10) /    2 /
!      DATA IMACH(11) /   27 /
!      DATA IMACH(12) / -128 /
!      DATA IMACH(13) /  127 /
!      DATA IMACH(14) /   54 /
!      DATA IMACH(15) / -101 /
!      DATA IMACH(16) /  127 /
!
!     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
!
!      DATA IMACH( 1) /    5 /
!      DATA IMACH( 2) /    6 /
!      DATA IMACH( 3) /    7 /
!      DATA IMACH( 4) /    6 /
!      DATA IMACH( 5) /   36 /
!      DATA IMACH( 6) /    5 /
!      DATA IMACH( 7) /    2 /
!      DATA IMACH( 8) /   35 /
!      DATA IMACH( 9) / "377777777777 /
!      DATA IMACH(10) /    2 /
!      DATA IMACH(11) /   27 /
!      DATA IMACH(12) / -128 /
!      DATA IMACH(13) /  127 /
!      DATA IMACH(14) /   62 /
!      DATA IMACH(15) / -128 /
!      DATA IMACH(16) /  127 /
!
!     MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
!     32-BIT INTEGER ARITHMETIC.
!
!      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) / -127 /
!      DATA IMACH(13) /  127 /
!      DATA IMACH(14) /   56 /
!      DATA IMACH(15) / -127 /
!      DATA IMACH(16) /  127 /
!
!     MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
!     16-BIT INTEGER ARITHMETIC.
!
!      DATA IMACH( 1) /    5 /
!      DATA IMACH( 2) /    6 /
!      DATA IMACH( 3) /    7 /
!      DATA IMACH( 4) /    6 /
!      DATA IMACH( 5) /   16 /
!      DATA IMACH( 6) /    2 /
!      DATA IMACH( 7) /    2 /
!      DATA IMACH( 8) /   15 /
!      DATA IMACH( 9) / 32767 /
!      DATA IMACH(10) /    2 /
!      DATA IMACH(11) /   24 /
!      DATA IMACH(12) / -127 /
!      DATA IMACH(13) /  127 /
!      DATA IMACH(14) /   56 /
!      DATA IMACH(15) / -127 /
!      DATA IMACH(16) /  127 /
!
!     MACHINE CONSTANTS FOR THE PRIME 50 SERIES SYSTEMS
!     WTIH 32-BIT INTEGERS AND 64V MODE INSTRUCTIONS,
!     SUPPLIED BY IGOR BRAY.
!
!      DATA IMACH( 1) /            1 /
!      DATA IMACH( 2) /            1 /
!      DATA IMACH( 3) /            2 /
!      DATA IMACH( 4) /            1 /
!      DATA IMACH( 5) /           32 /
!      DATA IMACH( 6) /            4 /
!      DATA IMACH( 7) /            2 /
!      DATA IMACH( 8) /           31 /
!      DATA IMACH( 9) / :17777777777 /
!      DATA IMACH(10) /            2 /
!      DATA IMACH(11) /           23 /
!      DATA IMACH(12) /         -127 /
!      DATA IMACH(13) /         +127 /
!      DATA IMACH(14) /           47 /
!      DATA IMACH(15) /       -32895 /
!      DATA IMACH(16) /       +32637 /
!
!     MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000.
!
!      DATA IMACH( 1) /     0 /
!      DATA IMACH( 2) /     0 /
!      DATA IMACH( 3) /     7 /
!      DATA IMACH( 4) /     0 /
!      DATA IMACH( 5) /    32 /
!      DATA IMACH( 6) /     1 /
!      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 /
!
!     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
!
!     NOTE THAT THE PUNCH UNIT, I1MACH(3), HAS BEEN SET TO 7
!     WHICH IS APPROPRIATE FOR THE UNIVAC-FOR SYSTEM.
!     IF YOU HAVE THE UNIVAC-FTN SYSTEM, SET IT TO 1.
!
!      DATA IMACH( 1) /    5 /
!      DATA IMACH( 2) /    6 /
!      DATA IMACH( 3) /    7 /
!      DATA IMACH( 4) /    6 /
!      DATA IMACH( 5) /   36 /
!      DATA IMACH( 6) /    6 /
!      DATA IMACH( 7) /    2 /
!      DATA IMACH( 8) /   35 /
!      DATA IMACH( 9) / O377777777777 /
!      DATA IMACH(10) /    2 /
!      DATA IMACH(11) /   27 /
!      DATA IMACH(12) / -128 /
!      DATA IMACH(13) /  127 /
!      DATA IMACH(14) /   60 /
!      DATA IMACH(15) /-1024 /
!      DATA IMACH(16) / 1023 /
!
!     MACHINE CONSTANTS FOR VAX.
!
!      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) / -127 /
!      DATA IMACH(13) /  127 /
!      DATA IMACH(14) /   56 /
!      DATA IMACH(15) / -127 /
!      DATA IMACH(16) /  127 /
      IF (I .LT. 1  .OR.  I .GT. 16) GO TO 999
      I1MACH=IMACH(I)
!/6S
!/7S
      IF(I.EQ.6) I1MACH=1
!/
      RETURN
  999 WRITE(OUTPUT,1999) I
 1999 FORMAT(' I1MACH - I OUT OF BOUNDS',I10)
      STOP
      END

SHAR_EOF
fi # end of overwriting check
if test -f 'f90d1mach.f90'
then
	echo shar: will not over-write existing file "'f90d1mach.f90'"
else
cat << \SHAR_EOF > 'f90d1mach.f90'
! This is a portable version of D1MACH that should work for any compiler
! that complies with the Fortran 90 specifications for the intrinsic
! EPSILON.

      DOUBLE PRECISION FUNCTION D1MACH(I) 
      INTEGER I, I1MACH
!                                                                       
!  DOUBLE-PRECISION MACHINE CONSTANTS                                   
!                                                                       
!  D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.           
!                                                                       
!  D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.           
!                                                                       
!  D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.                 
!                                                                       
!  D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.                 
!                                                                       
!  D1MACH( 5) = LOG10(B)                                                
!                                                                       
      SELECT CASE (I) 
        CASE (1) 
           D1MACH = TINY(1D0) 
        CASE (2) 
           D1MACH = HUGE(1D0) 
        CASE (3) 
           D1MACH = EPSILON(1D0)/RADIX(1D0) 
        CASE (4) 
           D1MACH = EPSILON(1D0) 
        CASE (5) 
           D1MACH = RADIX(1D0) 
           D1MACH = LOG10(D1MACH) 
        CASE DEFAULT 
        WRITE(I1MACH(2),'(A,I10)') ' D1MACH - I OUT OF BOUNDS', I 
        STOP 
      END SELECT 
      END                                           

SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Src'
then
	mkdir 'Src'
fi
cd 'Src'
if test ! -d 'Sp'
then
	mkdir 'Sp'
fi
cd 'Sp'
if test -f 'intaritb.f90'
then
	echo shar: will not over-write existing file "'intaritb.f90'"
else
cat << \SHAR_EOF > 'intaritb.f90'
! This Fortran 90 module provides support for overloading operators
! using INTLIB.

MODULE INTERVAL_ARITHMETIC

   USE IVL_DEF
   IMPLICIT NONE

      TYPE(INTERVAL), PARAMETER:: REAL_LINE=INTERVAL(-HUGE(1D0),HUGE(1D0))
      DOUBLE PRECISION, PARAMETER :: ZERO = 0D0, ONE=1D0

 !Interface to the basic operations

      INTERFACE OPERATOR(+)
         MODULE PROCEDURE ADD_F90,&
                          REAL_PLUS_IVL_F90,&
                          IVL_PLUS_REAL_F90,&
                          IVL_PLUS_INTEGER_F90,&
                          INTEGER_PLUS_IVL_F90
      END INTERFACE

      INTERFACE OPERATOR(-)
         MODULE PROCEDURE SUB_F90,&
                          REAL_MINUS_IVL_F90,&
                          IVL_MINUS_REAL_F90,&
                          INTEGER_MINUS_IVL_F90,&
                          IVL_MINUS_INTEGER_F90
      END INTERFACE

      INTERFACE OPERATOR(-)
         MODULE PROCEDURE INEG_F90
      END INTERFACE

      INTERFACE OPERATOR(*)
         MODULE PROCEDURE MULT_F90,&
                          REAL_TIMES_IVL_F90,&
                          IVL_TIMES_REAL_F90,&
                          INTEGER_TIMES_IVL_F90,&
                          IVL_TIMES_INTEGER_F90
      END INTERFACE

      INTERFACE OPERATOR(/)
         MODULE PROCEDURE DIV_F90,&
                          REAL_OVER_IVL_F90,&
                          IVL_OVER_REAL_F90,&
                          INTEGER_OVER_IVL_F90,&
                          IVL_OVER_INTEGER_F90
      END INTERFACE

 !Interface to the elementary functions

      INTERFACE OPERATOR(**)
         MODULE PROCEDURE POWER_F90, IIPOWR_F90,&
                          REAL_TO_IVL_F90, IVL_TO_REAL_F90, &
                          IGR_TO_IVL_F90
      END INTERFACE

      INTERFACE ACOS
         MODULE PROCEDURE IACOS_F90
      END INTERFACE

      INTERFACE ACOT
         MODULE PROCEDURE IACOT_F90
      END INTERFACE

      INTERFACE ASIN
         MODULE PROCEDURE IASIN_F90
      END INTERFACE

      INTERFACE ATAN
         MODULE PROCEDURE IATAN_F90
      END INTERFACE

      INTERFACE COS
         MODULE PROCEDURE ICOS_F90
      END INTERFACE

      INTERFACE COT
         MODULE PROCEDURE ICOT_F90
      END INTERFACE

      INTERFACE EXP
         MODULE PROCEDURE IEXP_F90
      END INTERFACE

      INTERFACE LOG
         MODULE PROCEDURE ILOG_F90
      END INTERFACE

      INTERFACE SIN
         MODULE PROCEDURE ISIN_F90
      END INTERFACE

      INTERFACE SINH
         MODULE PROCEDURE ISINH_F90
      END INTERFACE

      INTERFACE SQRT
         MODULE PROCEDURE ISQRT_F90
      END INTERFACE

       INTERFACE TAN
         MODULE PROCEDURE ITAN_F90
      END INTERFACE

!Interface to utility operations

      INTERFACE OPERATOR(.IS.)
         MODULE PROCEDURE ICAP_F90
      END INTERFACE

      INTERFACE OPERATOR(.CH.)
         MODULE PROCEDURE IHULL_F90, IHULL_R_I, IHULL_I_R, IHULL_R_R, &
                          IHULL_N_N, IHULL_I_N, IHULL_N_I, &
                          IHULL_N_R, IHULL_R_N
      END INTERFACE

      INTERFACE OPERATOR(.SB.)
         MODULE PROCEDURE IILEI_F90
      END INTERFACE

      INTERFACE OPERATOR(.SP.)
         MODULE PROCEDURE IIGEI_F90
      END INTERFACE

      INTERFACE OPERATOR(.DJ.)
         MODULE PROCEDURE IDISJ_F90
      END INTERFACE

      INTERFACE OPERATOR(.IN.)
         MODULE PROCEDURE IRLEI_F90, IRILEI, IIILEI
      END INTERFACE

      INTERFACE OPERATOR(.LT.)
         MODULE PROCEDURE INTINTLT_F90, REALINTLT_F90, INTREALLT_F90,&
                                        IGRINTLT_F90,  INTIGRLT_F90
      END INTERFACE

      INTERFACE OPERATOR(.GT.)
         MODULE PROCEDURE INTINTGT_F90, REALINTGT_F90, INTREALGT_F90,&
                                        IGRINTGT_F90,  INTIGRGT_F90
      END INTERFACE

      INTERFACE OPERATOR(.LE.)
         MODULE PROCEDURE INTINTLE_F90, REALINTLE_F90, INTREALLE_F90,&
                                        IGRINTLE_F90,  INTIGRLE_F90
      END INTERFACE

      INTERFACE OPERATOR(.GE.)
         MODULE PROCEDURE INTINTGE_F90, REALINTGE_F90, INTREALGE_F90,&
                                        IGRINTGE_F90,  INTIGRGE_F90
      END INTERFACE

      INTERFACE OPERATOR(.NE.)
         MODULE PROCEDURE INTINTNE, REALINTNE, INTREALNE, IGRINTNE, &
                          INTIGRNE
      END INTERFACE
      
      INTERFACE OPERATOR(.EQ.)
         MODULE PROCEDURE INTINTEQ, REALINTEQ, INTREALEQ, IGRINTEQ, &
                          INTIGREQ
      END INTERFACE

      INTERFACE ABS
         MODULE PROCEDURE IVLABS_F90
      END INTERFACE

      INTERFACE IWID
         MODULE PROCEDURE IWID_F90
      END INTERFACE

      INTERFACE WID
         MODULE PROCEDURE IWID_F90
      END INTERFACE

      INTERFACE MAG
         MODULE PROCEDURE INTABS_F90
      END INTERFACE

      INTERFACE MAX
         MODULE PROCEDURE IVLMAX_F90, IVRMAX, RIVMAX, IVIMAX, IIVMAX
      END INTERFACE

      INTERFACE MIN
         MODULE PROCEDURE IVLMIN_F90, IVRMIN, RIVMIN, IVIMIN, IIVMIN
      END INTERFACE

      INTERFACE MIG
         MODULE PROCEDURE IMIG_F90
      END INTERFACE

      INTERFACE IMID
         MODULE PROCEDURE IMID_F90
      END INTERFACE

      INTERFACE MID
         MODULE PROCEDURE IMID_F90
      END INTERFACE

 ! Overloading assignment

      INTERFACE ASSIGNMENT (=)
         MODULE PROCEDURE INTEGER_TO_INTERVAL,&
                          DOUBLE_TO_INTERVAL
      END INTERFACE

 ! Type conversions

      INTERFACE IVL
         MODULE PROCEDURE IVL1_F90, IVL2_F90, IVL1I_F90, IVL2I_F90,&
                          IVL2DI_F90, IVL2ID_F90, IVL_IVL
      END INTERFACE

   ! Explicit conversion functions (mostly used internally for
   ! conversion to and from INTLIB argument lists).

   ! Double prec. array to interval     IDBLA
   ! Interval to double prec. array     DBLA
   ! Double to double array             DBLAD
   ! Integer to double array            DBLAN

 ! Additional functions, compatible with Fortran-SC

      INTERFACE SUP
         MODULE PROCEDURE SUP_F90
      END INTERFACE

      INTERFACE INF
         MODULE PROCEDURE INF_F90
      END INTERFACE

CONTAINS

 ! Fortran 90 version of the INTLIB routine RNDOUT, for efficiency
 ! (The elementary operations in INTLIB are redefined here, too,
 ! for efficiency.)
 
      SUBROUTINE RNDOUT_F90(X,RNDDWN,RNDUP) 
         IMPLICIT NONE
         TYPE(INTERVAL) :: X
         LOGICAL RNDDWN, RNDUP
         
         DOUBLE PRECISION MXULP, TTINY2, TOL0
         COMMON /MACH1/ MXULP, TTINY2, TOL0
 
         DOUBLE PRECISION TINY, TEST
         COMMON /MACH2/ TINY, TEST
         
         IF (RNDDWN) THEN
            IF (X%LOWER.GE.TEST) THEN
               X%LOWER = (1.D0 - MXULP ) * X%LOWER
            ELSE IF (X%LOWER.LE.-TEST) THEN
               X%LOWER = (1D0 + MXULP ) * X%LOWER
            ELSE IF (X%LOWER.LE.0.D0) THEN
               X%LOWER = -TEST
            ELSE
               X%LOWER = 0.D0
            END IF
         END IF

         IF (RNDUP) THEN
            IF (X%UPPER.GE.TEST) THEN
               X%UPPER = (1.D0 + MXULP )* X%UPPER
            ELSE IF (X%UPPER.LE.-TEST) THEN
               X%UPPER = (1.D0 - MXULP ) * X%UPPER
            ELSE IF(X%UPPER.GE.0D0) THEN
               X%UPPER = TEST
            ELSE
               X%UPPER = 0.D0
            ENDIF
        END IF
     END SUBROUTINE RNDOUT_F90
        
 ! Basic operation Fortran 77 calls

      FUNCTION ADD_F90(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL) :: ADD_F90
         TYPE(INTERVAL), INTENT(IN) :: A, B
         
         ADD_F90%LOWER = A%LOWER + B%LOWER
         ADD_F90%UPPER = A%UPPER + B%UPPER

         CALL RNDOUT_F90(ADD_F90, &
                         (A%LOWER.NE.0D0).AND.(B%LOWER.NE.0D0), &
                         (A%UPPER.NE.0D0).AND.(B%UPPER.NE.0D0) )
      END FUNCTION ADD_F90

      FUNCTION REAL_PLUS_IVL_F90(A, B)
         IMPLICIT NONE
         DOUBLE PRECISION, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         TYPE(INTERVAL) :: REAL_PLUS_IVL_F90
 
         REAL_PLUS_IVL_F90%LOWER = A + B%LOWER
         REAL_PLUS_IVL_F90%UPPER = A + B%UPPER

         CALL RNDOUT_F90(REAL_PLUS_IVL_F90, &
                         (A.NE.0D0).AND.(B%LOWER.NE.0D0), &
                         (A.NE.0D0).AND.(B%UPPER.NE.0D0) )
      END FUNCTION REAL_PLUS_IVL_F90

      FUNCTION IVL_PLUS_REAL_F90(A, B)
         IMPLICIT NONE
         DOUBLE PRECISION, INTENT(IN) :: B
         TYPE(INTERVAL), INTENT(IN) :: A
         TYPE(INTERVAL) :: IVL_PLUS_REAL_F90
         
         IVL_PLUS_REAL_F90%LOWER = B + A%LOWER
         IVL_PLUS_REAL_F90%UPPER = B + A%UPPER

         CALL RNDOUT_F90(IVL_PLUS_REAL_F90, &
                         (B.NE.0D0).AND.(A%LOWER.NE.0D0), &
                         (B.NE.0D0).AND.(A%UPPER.NE.0D0) )       
      END FUNCTION IVL_PLUS_REAL_F90

      FUNCTION INTEGER_PLUS_IVL_F90(A, B)
         IMPLICIT NONE
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         TYPE(INTERVAL) :: INTEGER_PLUS_IVL_F90
         DOUBLE PRECISION T
         
         T = DBLE(A)
         INTEGER_PLUS_IVL_F90%LOWER = T + B%LOWER
         INTEGER_PLUS_IVL_F90%UPPER = T + B%UPPER

         CALL RNDOUT_F90(INTEGER_PLUS_IVL_F90, &
                         (T.NE.0D0).AND.(B%LOWER.NE.0D0), &
                         (T.NE.0D0).AND.(B%UPPER.NE.0D0) )
      END FUNCTION INTEGER_PLUS_IVL_F90

      FUNCTION IVL_PLUS_INTEGER_F90(A, B)
         IMPLICIT NONE
         INTEGER, INTENT(IN) :: B
         TYPE(INTERVAL), INTENT(IN) :: A
         TYPE(INTERVAL) :: IVL_PLUS_INTEGER_F90
         DOUBLE PRECISION T

         T = DBLE(B)
         IVL_PLUS_INTEGER_F90%LOWER = T + A%LOWER
         IVL_PLUS_INTEGER_F90%UPPER = T + A%UPPER
          
         CALL RNDOUT_F90(IVL_PLUS_INTEGER_F90, &
                         (T.NE.0D0).AND.(A%LOWER.NE.0D0), &
                         (T.NE.0D0).AND.(A%UPPER.NE.0D0) )           
      END FUNCTION IVL_PLUS_INTEGER_F90

      FUNCTION SUB_F90(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL) :: SUB_F90
         TYPE(INTERVAL), INTENT(IN) :: A, B
         DOUBLE PRECISION TA1, TA2, TB1, TB2
         
         TA1 = A%LOWER;  TA2 = A%UPPER
         TB1 = B%LOWER;  TB2 = B%UPPER

         SUB_F90%LOWER = TA1 - TB2
         SUB_F90%UPPER = TA2 - TB1

         CALL RNDOUT_F90(SUB_F90, (TB2.NE.0D0), (TB1.NE.0D0) )
      END FUNCTION SUB_F90

      FUNCTION REAL_MINUS_IVL_F90(A,B)
         IMPLICIT NONE
         DOUBLE PRECISION, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         TYPE(INTERVAL):: REAL_MINUS_IVL_F90
         REAL_MINUS_IVL_F90 = SUB_F90(INTERVAL(A,A),B)
      END FUNCTION REAL_MINUS_IVL_F90

      FUNCTION IVL_MINUS_REAL_F90(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(IN) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         TYPE(INTERVAL):: IVL_MINUS_REAL_F90
         IVL_MINUS_REAL_F90 = SUB_F90(A,INTERVAL(B,B))
      END FUNCTION IVL_MINUS_REAL_F90

      FUNCTION INTEGER_MINUS_IVL_F90(A,B)
         IMPLICIT NONE
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         TYPE(INTERVAL):: INTEGER_MINUS_IVL_F90
         INTEGER_MINUS_IVL_F90 = SUB_F90(INTERVAL(A,A),B)
      END FUNCTION INTEGER_MINUS_IVL_F90

      FUNCTION IVL_MINUS_INTEGER_F90(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(IN) :: A
         INTEGER, INTENT(IN) :: B
         TYPE(INTERVAL):: IVL_MINUS_INTEGER_F90
         IVL_MINUS_INTEGER_F90 = SUB_F90(A,INTERVAL(B,B))
      END FUNCTION IVL_MINUS_INTEGER_F90

      FUNCTION INEG_F90(A)
         IMPLICIT NONE
         TYPE(INTERVAL) INEG_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         TYPE(INTERVAL) :: T

         T%LOWER = A%LOWER
         T%UPPER = A%UPPER
         INEG_F90%LOWER = -T%UPPER
         INEG_F90%UPPER = -T%LOWER
         CALL RNDOUT_F90 (INEG_F90, .TRUE.,.TRUE.)
      END FUNCTION INEG_F90

      FUNCTION MULT_F90(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL) MULT_F90
         TYPE(INTERVAL), INTENT(IN) :: A, B
         DOUBLE PRECISION TEMP, A1, B1
         TYPE(INTERVAL) :: AA, BB
        !  Pictures for cases.
         IF ((ZERO .LE. A%LOWER) .AND. (ZERO .LE. B%LOWER)) THEN
        !  Case 1   ---------------- 0 -----------------
        !        A:                    [==========]
        !        B:                         [===========]
            MULT_F90%LOWER = A%LOWER * B%LOWER
            MULT_F90%UPPER = A%UPPER * B%UPPER
         ELSE IF ((A%LOWER .LT. ZERO) .AND. (ZERO .LT. A%UPPER) &
                                .AND. (ZERO .LE. B%LOWER)) THEN
        !  Case 2   ---------------- 0 -----------------
        !        A:            [==================]
        !        B:                         [===========]
            MULT_F90%LOWER = A%LOWER * B%UPPER
            MULT_F90%UPPER = A%UPPER * B%UPPER
         ELSE IF ((A%UPPER .LE. ZERO) .AND. (ZERO .LE. B%LOWER)) THEN
        !  Case 3   ---------------- 0 -----------------
        !        A:    [==========]
        !        B:                         [===========]
            B1 = B%LOWER
            MULT_F90%LOWER = A%LOWER * B%UPPER
            MULT_F90%UPPER = A%UPPER * B1
         ELSE IF ((ZERO .LE. A%LOWER) .AND. (B%LOWER .LT. ZERO) &
                                .AND. (ZERO .LT. B%UPPER)) THEN
        !  Case 4   ---------------- 0 -----------------
        !        A:                     [==========]
        !        B:               [===========]
            MULT_F90%LOWER = A%UPPER * B%LOWER
            MULT_F90%UPPER = A%UPPER * B%UPPER
         ELSE IF ((A%UPPER .LE. ZERO) .AND. (B%LOWER .LT. ZERO) &
                                .AND. (ZERO .LT. B%UPPER)) THEN
        !  Case 5   ---------------- 0 -----------------
        !        A:    [==========]
        !        B           [===========]
            A1 = A%LOWER
            B1 = B%LOWER
            MULT_F90%LOWER = A%LOWER * B%UPPER
            MULT_F90%UPPER = A1 * B1
         ELSE IF ((ZERO .LE. A%LOWER) .AND. (B%UPPER .LE. ZERO)) THEN
        !  Case 6   ---------------- 0 -----------------
        !        A:                    [==========]
        !        B:  [===========]
            A1 = A%LOWER
            MULT_F90%LOWER = A%UPPER * B%LOWER
            MULT_F90%UPPER = A1 * B%UPPER
         ELSE IF ((A%LOWER .LT. ZERO) .AND. (ZERO .LT. A%UPPER) &
                                .AND. (B%UPPER .LE. ZERO)) THEN
        !  Case 7   ---------------- 0 -----------------
        !        A:              [==========]
        !        B:  [===========]
            A1 = A%LOWER
            B1 = B%LOWER
            MULT_F90%LOWER = A%UPPER * B%LOWER
            MULT_F90%UPPER = A1 * B1
         ELSE IF ((A%UPPER .LE. ZERO) .AND. (B%UPPER .LE. ZERO)) THEN
        !  Case 8   ---------------- 0 -----------------
        !        A:    [==========]
        !        B:   [===========]
            A1 = A%LOWER
            B1 = B%LOWER
            MULT_F90%LOWER = A%UPPER * B%UPPER
            MULT_F90%UPPER = A1 * B1
         ELSE
        !  Case 9   ---------------- 0 -----------------
        !        A:              [==========]
        !        B:              [===========]
        !  Must check two cases.
            AA%LOWER = A%LOWER
            AA%UPPER = A%UPPER
            BB%LOWER = B%LOWER
            BB%UPPER = B%UPPER
            MULT_F90%LOWER = AA%LOWER * BB%UPPER
            TEMP = AA%UPPER * BB%LOWER
            IF (TEMP .LT. MULT_F90%LOWER) THEN
               MULT_F90%LOWER = TEMP
            ELSE
            END IF
            MULT_F90%UPPER = AA%LOWER * BB%LOWER
            TEMP = AA%UPPER * BB%UPPER
            IF (TEMP .GT. MULT_F90%UPPER) THEN
              MULT_F90%UPPER = TEMP
            ELSE
            END IF
         END IF
         CALL RNDOUT_F90(MULT_F90,.TRUE.,.TRUE.)  
      END FUNCTION MULT_F90

      FUNCTION REAL_TIMES_IVL_F90(R, B)
         IMPLICIT NONE
         DOUBLE PRECISION, INTENT(IN) :: R
         TYPE(INTERVAL), INTENT(IN) :: B
         TYPE(INTERVAL) :: REAL_TIMES_IVL_F90
         DOUBLE PRECISION T1, T2
         LOGICAL RNDDWN, RNDUP
         IF ((R.EQ.0D0).OR.((B%LOWER.EQ.0D0).AND.& 
                            (B%UPPER.EQ.0D0))) THEN
            REAL_TIMES_IVL_F90%LOWER = 0D0
            REAL_TIMES_IVL_F90%UPPER = 0D0
            RETURN
         END IF
         T1 = B%LOWER
         T2 = B%UPPER
         RNDDWN = .TRUE.
         RNDUP = .TRUE.
         IF (T1.EQ.0D0) THEN
            IF (R.LT.0D0) THEN
               REAL_TIMES_IVL_F90%LOWER = R * T2
               REAL_TIMES_IVL_F90%UPPER = 0D0
               RNDUP = .FALSE.
            ELSE
               REAL_TIMES_IVL_F90%LOWER = 0D0
               REAL_TIMES_IVL_F90%UPPER = R * T2
               RNDDWN = .FALSE.
            END IF
         ELSE IF (T2.EQ.0D0) THEN
            IF (R.LT.0D0) THEN
               REAL_TIMES_IVL_F90%LOWER = 0D0
               REAL_TIMES_IVL_F90%UPPER = R * T1
               RNDDWN = .FALSE.
            ELSE
               REAL_TIMES_IVL_F90%LOWER = R * T1
               REAL_TIMES_IVL_F90%UPPER = 0D0
               RNDUP = .FALSE.
            END IF
         ELSE IF (R.GT.0D0) THEN
            REAL_TIMES_IVL_F90%LOWER = R * T1
            REAL_TIMES_IVL_F90%UPPER = R * T2
         ELSE
            REAL_TIMES_IVL_F90%LOWER = R * T2
            REAL_TIMES_IVL_F90%UPPER = R * T1
         END IF
         CALL RNDOUT_F90(REAL_TIMES_IVL_F90,RNDDWN,RNDUP)
      END FUNCTION REAL_TIMES_IVL_F90

      FUNCTION IVL_TIMES_REAL_F90(A, B)
         IMPLICIT NONE
         DOUBLE PRECISION, INTENT(IN) :: B
         TYPE(INTERVAL), INTENT(IN) :: A
         TYPE(INTERVAL) :: IVL_TIMES_REAL_F90
         IVL_TIMES_REAL_F90 = REAL_TIMES_IVL_F90(B,A)
      END FUNCTION IVL_TIMES_REAL_F90

      FUNCTION INTEGER_TIMES_IVL_F90(A, B)
         IMPLICIT NONE
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         TYPE(INTERVAL) :: INTEGER_TIMES_IVL_F90
         INTEGER_TIMES_IVL_F90 = REAL_TIMES_IVL_F90(DBLE(A),B)
      END FUNCTION INTEGER_TIMES_IVL_F90

      FUNCTION IVL_TIMES_INTEGER_F90(A, B)
         IMPLICIT NONE
         INTEGER, INTENT(IN) ::  B
         TYPE(INTERVAL), INTENT(IN) :: A
         TYPE(INTERVAL) :: IVL_TIMES_INTEGER_F90
         IVL_TIMES_INTEGER_F90 = REAL_TIMES_IVL_F90(DBLE(B),A)
      END FUNCTION IVL_TIMES_INTEGER_F90

      FUNCTION DIV_F90(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL) :: DIV_F90
         TYPE(INTERVAL), INTENT(IN) :: A, B
         DOUBLE PRECISION, DIMENSION(2) :: X
         INTERFACE
            SUBROUTINE ERRTST(X)
               DOUBLE PRECISION, DIMENSION(2) :: X
            END SUBROUTINE ERRTST
         END INTERFACE
         TYPE(INTERVAL) :: C

         INTEGER ISIG, IERR, IROUT, IPRTCL, ISEVER, IERPUN
         COMMON /ERFLGS/ ISIG, IERR, IROUT, IPRTCL, ISEVER, IERPUN
       
       ! Identifying code for this routine --
         IROUT = 3
         IERR = 0
         
       ! Do usual interval division if zero is not in the denominator
         IF (B%LOWER.GT.ZERO) THEN
            C%LOWER = ONE/B%UPPER
            C%UPPER = ONE/B%LOWER
            CALL RNDOUT_F90(C,.TRUE.,.TRUE.)
            DIV_F90 =  MULT_F90(A,C)
         ELSE IF (B%UPPER.LT.ZERO) THEN
            C%LOWER = ONE/B%UPPER
            C%UPPER = ONE/B%LOWER
            CALL RNDOUT_F90(C,.TRUE.,.TRUE.)
            DIV_F90 = MULT_F90(A,C)
         ELSE
            IERR = 6
            X = DBLAI(B)
            CALL ERRTST(X)
            DIV_F90 = REAL_LINE
         END IF
      END FUNCTION DIV_F90

      FUNCTION REAL_OVER_IVL_F90(A,B)
         IMPLICIT NONE
         DOUBLE PRECISION, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         TYPE(INTERVAL) :: REAL_OVER_IVL_F90
         REAL_OVER_IVL_F90 = DIV_F90(INTERVAL(A,A),B)
      END FUNCTION REAL_OVER_IVL_F90

      FUNCTION IVL_OVER_REAL_F90(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(IN) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         TYPE(INTERVAL) :: IVL_OVER_REAL_F90
         IVL_OVER_REAL_F90 = DIV_F90(A,INTERVAL(B,B))
      END FUNCTION IVL_OVER_REAL_F90

      FUNCTION INTEGER_OVER_IVL_F90(A,B)
         IMPLICIT NONE
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         TYPE(INTERVAL) :: INTEGER_OVER_IVL_F90
         INTEGER_OVER_IVL_F90 = DIV_F90(INTERVAL(A,A),B)
      END FUNCTION INTEGER_OVER_IVL_F90

      FUNCTION IVL_OVER_INTEGER_F90(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(IN) :: A
         INTEGER, INTENT(IN) :: B
         TYPE(INTERVAL):: IVL_OVER_INTEGER_F90
         IVL_OVER_INTEGER_F90 = DIV_F90(A,INTERVAL(B,B))
      END FUNCTION IVL_OVER_INTEGER_F90

 ! Elementary function Fortran 77 calls (interfaces to INTLIB)

      FUNCTION POWER_F90(A,N)
         IMPLICIT NONE
         TYPE(INTERVAL) POWER_F90
         TYPE(INTERVAL), INTENT(IN):: A
         INTEGER, INTENT(IN):: N
         INTERFACE
            SUBROUTINE POWER(A,N,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A
               INTEGER N
               DOUBLE PRECISION, DIMENSION(2) :: VALUE
            END SUBROUTINE POWER
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP
         CALL POWER (DBLAI(A),N,TMP)
         POWER_F90 = IDBLA(TMP)
      END FUNCTION POWER_F90

      FUNCTION IIPOWR_F90(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL) IIPOWR_F90
         TYPE(INTERVAL), INTENT(IN) :: A, B
         INTERFACE
            SUBROUTINE IIPOWR(A,B,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, B, VALUE
            END SUBROUTINE IIPOWR
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP
         CALL IIPOWR (DBLAI(A),DBLAI(B),TMP)
         IIPOWR_F90 = IDBLA(TMP)
      END FUNCTION IIPOWR_F90

      FUNCTION REAL_TO_IVL_F90(A,B)
         IMPLICIT NONE
         DOUBLE PRECISION, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         INTERFACE
            SUBROUTINE IIPOWR(A,B,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, B, VALUE
            END SUBROUTINE IIPOWR
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP
         TYPE(INTERVAL):: REAL_TO_IVL_F90
         CALL IIPOWR(DBLAD(A),DBLAI(B),TMP)
         REAL_TO_IVL_F90 = IDBLA(TMP)
      END FUNCTION REAL_TO_IVL_F90

      FUNCTION IVL_TO_REAL_F90(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(IN) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         TYPE(INTERVAL):: IVL_TO_REAL_F90
         INTERFACE
            SUBROUTINE IIPOWR(A,B,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, B, VALUE
            END SUBROUTINE IIPOWR
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP
         CALL IIPOWR(DBLAI(A),DBLAD(B),TMP)
         IVL_TO_REAL_F90 = IDBLA(TMP)
      END FUNCTION IVL_TO_REAL_F90

      FUNCTION IGR_TO_IVL_F90(A,B)
         IMPLICIT NONE
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         TYPE(INTERVAL):: IGR_TO_IVL_F90
         INTERFACE
            SUBROUTINE IIPOWR(A,B,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, B, VALUE
            END SUBROUTINE IIPOWR
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP
         CALL IIPOWR(DBLAN(A),DBLAI(B),TMP)
         IGR_TO_IVL_F90 = IDBLA(TMP)
      END FUNCTION IGR_TO_IVL_F90

      FUNCTION ICOS_F90(A)
         IMPLICIT NONE
         TYPE(INTERVAL) ICOS_F90
         TYPE(INTERVAL), INTENT(IN):: A
         INTERFACE
            SUBROUTINE ICOS(A,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, VALUE
            END SUBROUTINE ICOS
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP
         CALL ICOS (DBLAI(A),TMP)
         ICOS_F90 = IDBLA(TMP)
      END FUNCTION ICOS_F90

      FUNCTION IEXP_F90(A)
        IMPLICIT NONE
        TYPE(INTERVAL) IEXP_F90
         TYPE(INTERVAL), INTENT(IN):: A
         INTERFACE
            SUBROUTINE IEXP(A,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, VALUE
            END SUBROUTINE IEXP
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP
         CALL IEXP (DBLAI(A),TMP)
         IEXP_F90 = IDBLA(TMP)
      END FUNCTION IEXP_F90

      FUNCTION ILOG_F90(A)
         IMPLICIT NONE
         TYPE(INTERVAL) ILOG_F90
         TYPE(INTERVAL), INTENT(IN):: A
         INTERFACE
            SUBROUTINE ILOG(A,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, VALUE
            END SUBROUTINE ILOG
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP
         CALL ILOG (DBLAI(A),TMP)
         ILOG_F90 = IDBLA(TMP)
      END FUNCTION ILOG_F90

      FUNCTION ISIN_F90(A)
         IMPLICIT NONE
         TYPE(INTERVAL) ISIN_F90
         TYPE(INTERVAL), INTENT(IN):: A
         INTERFACE
            SUBROUTINE ISIN(A,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, VALUE
            END SUBROUTINE ISIN
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP
         CALL ISIN (DBLAI(A),TMP)
         ISIN_F90 = IDBLA(TMP)
      END FUNCTION ISIN_F90

      FUNCTION ITAN_F90(A)
         IMPLICIT NONE
         TYPE(INTERVAL) ITAN_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         INTERFACE
            SUBROUTINE ISIN(A,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, VALUE
            END SUBROUTINE ISIN
            SUBROUTINE ICOS(A,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, VALUE
            END SUBROUTINE ICOS
            SUBROUTINE IDIV(A,B,C)
               DOUBLE PRECISION, DIMENSION(2) :: A, B, C
            END SUBROUTINE IDIV
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP1, TMP2
         TMP2 = DBLAI(A)
         CALL ISIN (TMP2, TMP1)
         CALL ICOS (TMP2, TMP2)
         CALL IDIV (TMP1, TMP2, TMP1)
         ITAN_F90 = IDBLA(TMP1)
      END FUNCTION ITAN_F90

      FUNCTION ICOT_F90(A)
         IMPLICIT NONE
         TYPE(INTERVAL) ICOT_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         INTERFACE
            SUBROUTINE ISIN(A,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, VALUE
            END SUBROUTINE ISIN
            SUBROUTINE ICOS(A,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, VALUE
            END SUBROUTINE ICOS
            SUBROUTINE IDIV(A,B,C)
               DOUBLE PRECISION, DIMENSION(2) :: A, B, C
            END SUBROUTINE IDIV
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP1, TMP2
         TMP2 = DBLAI(A)
         CALL ISIN (TMP2, TMP1)
         CALL ICOS (TMP2, TMP2)
         CALL IDIV (TMP2, TMP1, TMP1)
         ICOT_F90 = IDBLA(TMP1)
      END FUNCTION ICOT_F90

      FUNCTION ISQRT_F90(A)
         IMPLICIT NONE
         TYPE(INTERVAL) ISQRT_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         INTERFACE
            SUBROUTINE ISQRT(A,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, VALUE
            END SUBROUTINE ISQRT
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP
         CALL ISQRT (DBLAI(A),TMP)
         ISQRT_F90 = IDBLA(TMP)
      END FUNCTION ISQRT_F90

      FUNCTION IATAN_F90(A)
         IMPLICIT NONE
         TYPE(INTERVAL) IATAN_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         INTERFACE
            SUBROUTINE IATAN(A,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, VALUE
            END SUBROUTINE IATAN
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP
         CALL IATAN (DBLAI(A),TMP)
         IATAN_F90 = IDBLA(TMP)
      END FUNCTION IATAN_F90

      FUNCTION IACOT_F90(A)
         IMPLICIT NONE
         TYPE(INTERVAL) IACOT_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         INTERFACE
            SUBROUTINE IACOT(A,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, VALUE
            END SUBROUTINE IACOT
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP
         CALL IACOT (DBLAI(A),TMP)
         IACOT_F90 = IDBLA(TMP)
      END FUNCTION IACOT_F90

      FUNCTION IASIN_F90(A)
         IMPLICIT NONE
         TYPE(INTERVAL) IASIN_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         INTERFACE
            SUBROUTINE IASIN(A,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, VALUE
            END SUBROUTINE IASIN
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP
         CALL IASIN (DBLAI(A),TMP)
         IASIN_F90 = IDBLA(TMP)
      END FUNCTION IASIN_F90

      FUNCTION IACOS_F90(A)
         IMPLICIT NONE
         TYPE(INTERVAL) IACOS_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         INTERFACE
            SUBROUTINE IACOS(A,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, VALUE
            END SUBROUTINE IACOS
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP
         CALL IACOS (DBLAI(A),TMP)
         IACOS_F90 = IDBLA(TMP)
      END FUNCTION IACOS_F90

      FUNCTION ISINH_F90(A)
         IMPLICIT NONE
         TYPE(INTERVAL) ISINH_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         INTERFACE
            SUBROUTINE ISINH(A,VALUE)
               DOUBLE PRECISION, DIMENSION(2) :: A, VALUE
            END SUBROUTINE ISINH
         END INTERFACE
         DOUBLE PRECISION, DIMENSION(2) :: TMP
         CALL ISINH (DBLAI(A),TMP)
         ISINH_F90 = IDBLA(TMP)
      END FUNCTION ISINH_F90

 ! Utility function Fortran 77 calls

      FUNCTION ICAP_F90(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL) ICAP_F90
         TYPE(INTERVAL), INTENT(IN) :: A, B
         TYPE(INTERVAL) :: T
         DOUBLE PRECISION, DIMENSION(2) :: X
         INTEGER ISIG, IERR, IROUT, IPRTCL, ISEVER, IERPUN
         COMMON /ERFLGS/ ISIG, IERR, IROUT, IPRTCL, ISEVER, IERPUN
         IROUT = 13
         IERR = 0
         T%LOWER = MAX (A%LOWER, B%LOWER)
         T%UPPER = MIN (A%UPPER, B%UPPER)
         ICAP_F90 = T
         IF (ICAP_F90%LOWER.GT.ICAP_F90%UPPER) THEN
            IERR=13
            X = DBLAI(ICAP_F90)
            CALL ERRTST(X)
            ICAP_F90 = IDBLA(X)
         END IF
      END FUNCTION ICAP_F90

      FUNCTION IHULL_F90(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL) IHULL_F90
         TYPE(INTERVAL), INTENT(IN) :: A, B
         TYPE(INTERVAL) :: T
         IF ( A%LOWER.GT.A%UPPER ) THEN
            IF (B%LOWER.GT.B%UPPER) THEN
               T=INTERVAL(MAX(A%UPPER,B%UPPER),MIN(A%LOWER,B%UPPER))
            ELSE
               T = B
            END IF
         ELSE IF ( B%LOWER.GT.B%UPPER ) THEN
            T = A
         ELSE
            T=INTERVAL(MIN(A%LOWER,B%LOWER),MAX(A%UPPER,B%UPPER))
         END IF
         IHULL_F90 = T
      END FUNCTION IHULL_F90

      FUNCTION IHULL_R_I(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL) IHULL_R_I
         DOUBLE PRECISION, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         IHULL_R_I = IHULL_F90(INTERVAL(A,A),B)
      END FUNCTION IHULL_R_I

      FUNCTION IHULL_I_R(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL) IHULL_I_R
         TYPE(INTERVAL), INTENT(IN) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         IHULL_I_R = IHULL_F90(A,INTERVAL(B,B))
      END FUNCTION IHULL_I_R

      FUNCTION IHULL_R_R(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL) IHULL_R_R
         DOUBLE PRECISION, INTENT(IN) :: A, B
         IHULL_R_R = IHULL_F90(INTERVAL(A,A),INTERVAL(B,B))
      END FUNCTION IHULL_R_R

      FUNCTION IHULL_N_I(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL) IHULL_N_I
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         IHULL_N_I = IHULL_F90(INTERVAL(A,A),B)
      END FUNCTION IHULL_N_I

      FUNCTION IHULL_I_N(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL) IHULL_I_N
         TYPE(INTERVAL), INTENT(IN) :: A
         INTEGER, INTENT(IN) :: B
         IHULL_I_N = IHULL_F90(A,INTERVAL(B,B))
      END FUNCTION IHULL_I_N

      FUNCTION IHULL_N_N(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL) IHULL_N_N
         INTEGER, INTENT(IN) :: A, B
         IHULL_N_N = IHULL_F90(INTERVAL(A,A),INTERVAL(B,B))
      END FUNCTION IHULL_N_N

      FUNCTION IHULL_N_R(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL) IHULL_N_R
         INTEGER, INTENT(IN) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         IHULL_N_R = IHULL_F90(INTERVAL(A,A),INTERVAL(B,B))
      END FUNCTION IHULL_N_R

      FUNCTION IHULL_R_N(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL) IHULL_R_N
         DOUBLE PRECISION, INTENT(IN) :: A
         INTEGER, INTENT(IN) :: B
         IHULL_R_N = IHULL_F90(INTERVAL(A,A),INTERVAL(B,B))
      END FUNCTION IHULL_R_N

      FUNCTION IILEI_F90(A,B) RESULT(L)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(IN):: A, B
         LOGICAL:: L
         L = (A%LOWER.GE.B%LOWER) .AND. (A%UPPER.LE.B%UPPER)
      END FUNCTION IILEI_F90

      FUNCTION IIGEI_F90(A,B)
         IMPLICIT NONE
         LOGICAL IIGEI_F90
         TYPE(INTERVAL), INTENT(IN) :: A, B
         IIGEI_F90 = IILEI_F90(B,A)
      END FUNCTION IIGEI_F90

      LOGICAL FUNCTION IRILEI(A,B)
         IMPLICIT NONE
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         IRILEI = IILEI_F90(INTERVAL(A,A),B)
      END FUNCTION IRILEI

      LOGICAL FUNCTION IIILEI(A,B)
         TYPE(INTERVAL), INTENT(IN) :: A, B
         IIILEI = A%LOWER.GT.B%LOWER .AND. A%UPPER.LT.B%UPPER
      END FUNCTION IIILEI

      FUNCTION IRLEI_F90(A,B) RESULT(L)
         IMPLICIT NONE
         DOUBLE PRECISION, INTENT(IN):: A
         TYPE(INTERVAL), INTENT(IN):: B
         LOGICAL:: L
         L = (A .GE. B%LOWER) .AND. (A .LE. B%UPPER)
      END FUNCTION IRLEI_F90

      FUNCTION IDISJ_F90(A,B) RESULT(L)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(IN):: A, B
         LOGICAL:: L
         L = (A%UPPER .LT. B%LOWER) .OR. (B%UPPER .LT. A%LOWER)
      END FUNCTION IDISJ_F90
      
 ! Interval-valued absolute value
 
      FUNCTION IVLABS_F90(A) RESULT(R)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(IN) :: A
         TYPE(INTERVAL) :: R
         TYPE(INTERVAL) :: TA, TMP
         TA = A
         TMP%LOWER = ABS(TA%LOWER)
         TMP%UPPER = ABS(TA%UPPER)
         IF (TMP%LOWER.LE.TMP%UPPER) THEN
            R = TMP
         ELSE
            R%LOWER = TMP%UPPER
            R%UPPER = TMP%LOWER
         END IF
         CALL RNDOUT_F90(R,.TRUE.,.TRUE.)
         IF ( ( TA%LOWER.LE.ZERO .AND. TA%UPPER.GE.ZERO ) .OR. &
                R%LOWER.LT. ZERO)  R%LOWER = ZERO
      END FUNCTION IVLABS_F90

 ! Not in INTLIB

      FUNCTION IVLMAX_F90(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(IN) :: A, B
         TYPE(INTERVAL) :: IVLMAX_F90
         IVLMAX_F90%LOWER = MAX(A%LOWER,B%LOWER)
         IVLMAX_F90%UPPER = MAX(A%UPPER,B%UPPER)
      END FUNCTION IVLMAX_F90

      FUNCTION IVRMAX(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(IN) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         TYPE(INTERVAL) :: IVRMAX
         IVRMAX%LOWER = MAX(A%LOWER,B)
         IVRMAX%UPPER = MAX(A%UPPER,B)
      END FUNCTION IVRMAX

      FUNCTION RIVMAX(A,B)
         IMPLICIT NONE
         DOUBLE PRECISION, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         TYPE(INTERVAL) :: RIVMAX
         RIVMAX%LOWER = MAX(B%LOWER,A)
         RIVMAX%UPPER = MAX(B%UPPER,A)
      END FUNCTION RIVMAX

      FUNCTION IVIMAX(A,B)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(IN) :: A
         INTEGER, INTENT(IN) :: B
         TYPE(INTERVAL) :: IVIMAX
         IVIMAX%LOWER = MAX(A%LOWER,DBLE(B))
         IVIMAX%UPPER = MAX(A%UPPER,DBLE(B))
      END FUNCTION IVIMAX

      FUNCTION IIVMAX(A,B)
         IMPLICIT NONE
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         TYPE(INTERVAL) :: IIVMAX
         IIVMAX%LOWER = MAX(B%LOWER,DBLE(A))
         IIVMAX%UPPER = MAX(B%UPPER,DBLE(A))
      END FUNCTION IIVMAX

  ! .LT.

      FUNCTION INTINTLT_F90(A,B)
         IMPLICIT NONE
         LOGICAL INTINTLT_F90
         TYPE(INTERVAL), INTENT(IN) :: A, B
         INTINTLT_F90 = A%UPPER.LT.B%LOWER
      END FUNCTION INTINTLT_F90

      FUNCTION REALINTLT_F90(A,B)
         IMPLICIT NONE
         LOGICAL REALINTLT_F90
         DOUBLE PRECISION, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         REALINTLT_F90 = A.LT.B%LOWER
      END FUNCTION REALINTLT_F90

      FUNCTION INTREALLT_F90(A,B)
         IMPLICIT NONE
         LOGICAL INTREALLT_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         INTREALLT_F90 = A%UPPER.LT.B
      END FUNCTION INTREALLT_F90

      FUNCTION IGRINTLT_F90(A,B)
         IMPLICIT NONE
         LOGICAL IGRINTLT_F90
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         IGRINTLT_F90 = A.LT.B%LOWER
      END FUNCTION IGRINTLT_F90

      FUNCTION INTIGRLT_F90(A,B)
         IMPLICIT NONE
         LOGICAL INTIGRLT_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         INTEGER, INTENT(IN) :: B
         INTIGRLT_F90 = A%UPPER.LT.B
      END FUNCTION INTIGRLT_F90

      FUNCTION IVLMIN_F90(A,B)
         TYPE(INTERVAL), INTENT(IN) :: A, B
         TYPE(INTERVAL) :: IVLMIN_F90
         IVLMIN_F90%LOWER = MIN(A%LOWER,B%LOWER)
         IVLMIN_F90%UPPER = MIN(A%UPPER,B%UPPER)
      END FUNCTION IVLMIN_F90

      FUNCTION IVRMIN(A,B)
         TYPE(INTERVAL), INTENT(IN) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         TYPE(INTERVAL) :: IVRMIN
         IVRMIN%LOWER = MIN(A%LOWER,B)
         IVRMIN%UPPER = MIN(A%UPPER,B)
      END FUNCTION IVRMIN

      FUNCTION RIVMIN(A,B)
         DOUBLE PRECISION, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         TYPE(INTERVAL) :: RIVMIN
         RIVMIN%LOWER = MIN(B%LOWER,A)
         RIVMIN%UPPER = MIN(B%UPPER,A)
      END FUNCTION RIVMIN

      FUNCTION IVIMIN(A,B)
         TYPE(INTERVAL), INTENT(IN) :: A
         INTEGER, INTENT(IN) :: B
         TYPE(INTERVAL) :: IVIMIN
         IVIMIN%LOWER = MIN(A%LOWER,DBLE(B))
         IVIMIN%UPPER = MIN(A%UPPER,DBLE(B))
      END FUNCTION IVIMIN

      FUNCTION IIVMIN(A,B)
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         TYPE(INTERVAL) :: IIVMIN
         IIVMIN%LOWER = MIN(B%LOWER,DBLE(A))
         IIVMIN%UPPER = MIN(B%UPPER,DBLE(A))
      END FUNCTION IIVMIN

  ! .GT.

      FUNCTION INTINTGT_F90(A,B)
         IMPLICIT NONE
         LOGICAL INTINTGT_F90
         TYPE(INTERVAL), INTENT(IN) :: A, B
         INTINTGT_F90 = A%LOWER.GT.B%UPPER
      END FUNCTION INTINTGT_F90

      FUNCTION REALINTGT_F90(A,B)
         IMPLICIT NONE
         LOGICAL REALINTGT_F90
         DOUBLE PRECISION, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         REALINTGT_F90 = A.GT.B%UPPER
      END FUNCTION REALINTGT_F90

      FUNCTION INTREALGT_F90(A,B)
         IMPLICIT NONE
         LOGICAL INTREALGT_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         INTREALGT_F90 = A%LOWER.GT.B
      END FUNCTION INTREALGT_F90

      FUNCTION IGRINTGT_F90(A,B)
         IMPLICIT NONE
         LOGICAL IGRINTGT_F90
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         IGRINTGT_F90 = A.GT.B%UPPER
      END FUNCTION IGRINTGT_F90

      FUNCTION INTIGRGT_F90(A,B)
         IMPLICIT NONE
         LOGICAL INTIGRGT_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         INTEGER, INTENT(IN) :: B
         INTIGRGT_F90 = A%LOWER.GT.B
      END FUNCTION INTIGRGT_F90

  ! .LE.

      FUNCTION INTINTLE_F90(A,B)
         IMPLICIT NONE
         LOGICAL INTINTLE_F90
         TYPE(INTERVAL), INTENT(IN) :: A, B
         INTINTLE_F90 = A%UPPER.LE.B%LOWER
      END FUNCTION INTINTLE_F90

      FUNCTION REALINTLE_F90(A,B)
         IMPLICIT NONE
         LOGICAL REALINTLE_F90
         DOUBLE PRECISION, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         REALINTLE_F90 = A.LE.B%LOWER
      END FUNCTION REALINTLE_F90

      FUNCTION INTREALLE_F90(A,B)
         IMPLICIT NONE
         LOGICAL INTREALLE_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         INTREALLE_F90 = A%UPPER.LE.B
      END FUNCTION INTREALLE_F90

      FUNCTION IGRINTLE_F90(A,B)
         IMPLICIT NONE
         LOGICAL IGRINTLE_F90
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         IGRINTLE_F90 = A.LE.B%LOWER
      END FUNCTION IGRINTLE_F90

      FUNCTION INTIGRLE_F90(A,B)
         IMPLICIT NONE
         LOGICAL INTIGRLE_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         INTEGER, INTENT(IN) :: B
         INTIGRLE_F90 = A%UPPER.LE.B
      END FUNCTION INTIGRLE_F90

  ! .GE.

      FUNCTION INTINTGE_F90(A,B)
         IMPLICIT NONE
         LOGICAL INTINTGE_F90
         TYPE(INTERVAL), INTENT(IN) :: A, B
         INTINTGE_F90 = A%LOWER.GE.B%UPPER
      END FUNCTION INTINTGE_F90

      FUNCTION REALINTGE_F90(A,B)
         IMPLICIT NONE
         LOGICAL REALINTGE_F90
         DOUBLE PRECISION, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         REALINTGE_F90 = A.GE.B%UPPER
      END FUNCTION REALINTGE_F90

      FUNCTION INTREALGE_F90(A,B)
         IMPLICIT NONE
         LOGICAL INTREALGE_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         INTREALGE_F90 = A%LOWER.GE.B
      END FUNCTION INTREALGE_F90

      FUNCTION IGRINTGE_F90(A,B)
         IMPLICIT NONE
         LOGICAL IGRINTGE_F90
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         IGRINTGE_F90 = A.GE.B%UPPER
      END FUNCTION IGRINTGE_F90

      FUNCTION INTIGRGE_F90(A,B)
         IMPLICIT NONE
         LOGICAL INTIGRGE_F90
         TYPE(INTERVAL), INTENT(IN) :: A
         INTEGER, INTENT(IN) :: B
         INTIGRGE_F90 = A%LOWER.GE.B
      END FUNCTION INTIGRGE_F90

  ! .NE.

      FUNCTION INTINTNE(A,B)
         IMPLICIT NONE
         LOGICAL INTINTNE
         TYPE(INTERVAL), INTENT(IN) :: A, B
         INTINTNE = (A%LOWER.NE.B%LOWER) .OR. (A%UPPER.NE.B%UPPER)
      END FUNCTION INTINTNE

      FUNCTION REALINTNE(A,B)
         IMPLICIT NONE
         LOGICAL REALINTNE
         DOUBLE PRECISION, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         REALINTNE = (A.NE.B%LOWER) .OR. (A.NE.B%UPPER)
      END FUNCTION REALINTNE

      FUNCTION INTREALNE(A,B)
         IMPLICIT NONE
         LOGICAL INTREALNE
         TYPE(INTERVAL), INTENT(IN) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         INTREALNE = (A%LOWER.NE.B) .OR. (A%UPPER.NE.B)
      END FUNCTION INTREALNE

      FUNCTION IGRINTNE(A,B)
         IMPLICIT NONE
         LOGICAL IGRINTNE
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         IGRINTNE = (DBLE(A).NE.B%LOWER) .OR. (DBLE(A).NE.B%UPPER)
      END FUNCTION IGRINTNE

      FUNCTION INTIGRNE(A,B)
         IMPLICIT NONE
         LOGICAL INTIGRNE
         TYPE(INTERVAL), INTENT(IN) :: A
         INTEGER, INTENT(IN) :: B
         INTIGRNE = (A%LOWER.NE.DBLE(B)) .OR. (A%UPPER.NE.DBLE(B))
      END FUNCTION INTIGRNE

  ! .EQ.

      FUNCTION INTINTEQ(A,B)
         IMPLICIT NONE
         LOGICAL INTINTEQ
         TYPE(INTERVAL), INTENT(IN) :: A, B
         INTINTEQ = (A%LOWER.EQ.B%LOWER) .AND. (A%UPPER.EQ.B%UPPER)
      END FUNCTION INTINTEQ

      FUNCTION REALINTEQ(A,B)
         IMPLICIT NONE
         LOGICAL REALINTEQ
         DOUBLE PRECISION, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         REALINTEQ = (A.EQ.B%LOWER) .AND. (A.EQ.B%UPPER)
      END FUNCTION REALINTEQ

      FUNCTION INTREALEQ(A,B)
         IMPLICIT NONE
         LOGICAL INTREALEQ
         TYPE(INTERVAL), INTENT(IN) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         INTREALEQ = (A%LOWER.EQ.B) .AND. (A%UPPER.EQ.B)
      END FUNCTION INTREALEQ

      FUNCTION IGRINTEQ(A,B)
         IMPLICIT NONE
         LOGICAL IGRINTEQ
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         IGRINTEQ = (DBLE(A).EQ.B%LOWER) .AND. (DBLE(A).EQ.B%UPPER)
      END FUNCTION IGRINTEQ

      FUNCTION INTIGREQ(A,B)
         IMPLICIT NONE
         LOGICAL INTIGREQ
         TYPE(INTERVAL), INTENT(IN) :: A
         INTEGER, INTENT(IN) :: B
         INTIGREQ = (A%LOWER.EQ.DBLE(B)) .AND. (A%UPPER.EQ.DBLE(B))
      END FUNCTION INTIGREQ

 ! Assignment of other data types to interval

      SUBROUTINE INTEGER_TO_INTERVAL (A,B)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(OUT) :: A
         INTEGER, INTENT(IN) :: B
         A = INTERVAL(B,B)
         IF (B.NE.0) CALL RNDOUT_F90(A,.TRUE.,.TRUE.)
       END SUBROUTINE INTEGER_TO_INTERVAL

      SUBROUTINE DOUBLE_TO_INTERVAL (A,B)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(OUT) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         A = INTERVAL(B,B)
         IF (B.NE.0) CALL RNDOUT_F90(A,.TRUE.,.TRUE.)
      END SUBROUTINE DOUBLE_TO_INTERVAL

   ! Internal data conversion routines for interfacing with INTLIB

      FUNCTION IDBLA(B) RESULT(A)
         IMPLICIT NONE
         TYPE(INTERVAL) :: A
         DOUBLE PRECISION, DIMENSION(2) :: B
         A%LOWER = B(1)
         A%UPPER = B(2)
      END FUNCTION IDBLA

      FUNCTION DBLAI(B) RESULT(A)
         IMPLICIT NONE
         DOUBLE PRECISION, DIMENSION(2) :: A
         TYPE(INTERVAL), INTENT(IN) :: B
         A(1) = B%LOWER
         A(2) = B%UPPER
      END FUNCTION DBLAI

      FUNCTION DBLAD(B) RESULT(A)
         IMPLICIT NONE
         DOUBLE PRECISION, DIMENSION(2) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         A(1) = B
         A(2) = B
      END FUNCTION DBLAD

      FUNCTION DBLAN(B) RESULT(A)
         IMPLICIT NONE
         DOUBLE PRECISION, DIMENSION(2) :: A
         INTEGER, INTENT(IN) :: B
         A(1) = DBLE(B)
         A(2) = DBLE(B)
      END FUNCTION DBLAN

      FUNCTION IVL1_F90(A) RESULT(R)
         IMPLICIT NONE
        DOUBLE PRECISION, INTENT(IN) :: A
         TYPE(INTERVAL) :: R
         R = INTERVAL(A,A)
         CALL RNDOUT_F90(R,.TRUE.,.TRUE.)
         END FUNCTION IVL1_F90

      FUNCTION IVL2_F90(A,B) RESULT(R)
         IMPLICIT NONE
         DOUBLE PRECISION, INTENT(IN) :: A, B
         TYPE(INTERVAL):: R
         R = INTERVAL(A,B)
         CALL RNDOUT_F90(R,.TRUE.,.TRUE.)
      END FUNCTION IVL2_F90

      FUNCTION IVL1I_F90(A) RESULT(R)
         IMPLICIT NONE
         INTEGER, INTENT(IN) :: A
         TYPE(INTERVAL) :: R
         R = INTERVAL(A,A)
         CALL RNDOUT_F90(R,.TRUE.,.TRUE.)
      END FUNCTION IVL1I_F90

      FUNCTION IVL2I_F90(A,B) RESULT(R)
         IMPLICIT NONE
         INTEGER, INTENT(IN) :: A, B
         TYPE(INTERVAL):: R
         R = INTERVAL(A,B)
         CALL RNDOUT_F90(R,.TRUE.,.TRUE.)
      END FUNCTION IVL2I_F90

      FUNCTION IVL2DI_F90(A,B) RESULT(R)
         IMPLICIT NONE
         DOUBLE PRECISION, INTENT(IN) :: A
         INTEGER, INTENT(IN) :: B
         TYPE(INTERVAL):: R
         R = INTERVAL(A,B)
         CALL RNDOUT_F90(R,.TRUE.,.TRUE.)
      END FUNCTION IVL2DI_F90

      FUNCTION IVL2ID_F90(A,B) RESULT(R)
         IMPLICIT NONE
         INTEGER, INTENT(IN) :: A
         DOUBLE PRECISION, INTENT(IN) :: B
         TYPE(INTERVAL):: R
         R = INTERVAL(A,B)
         CALL RNDOUT_F90(R,.TRUE.,.TRUE.)
      END FUNCTION IVL2ID_F90

      FUNCTION IVL_IVL(A) RESULT(R)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(IN) :: A
         TYPE(INTERVAL) :: R
         R = A
      END FUNCTION IVL_IVL

 ! mag, abs, max, iwid, mig, imid --

      FUNCTION INTABS_F90(A) RESULT(D)
         IMPLICIT NONE
         TYPE(INTERVAL) :: A
         DOUBLE PRECISION D
         
         TYPE(INTERVAL) :: T
         
         T%UPPER = ABS(A%LOWER)
         CALL RNDOUT_F90(T,.FALSE.,.TRUE.)
         T%LOWER=T%UPPER
         T%UPPER = ABS(A%UPPER)
         CALL RNDOUT_F90(T,.FALSE.,.TRUE.)
         D = MAX (T%LOWER, T%UPPER)
      END FUNCTION INTABS_F90

      FUNCTION IWID_F90(A) RESULT(B)
         IMPLICIT NONE
         TYPE(INTERVAL) :: A
         DOUBLE PRECISION B
         TYPE(INTERVAL) :: T
         B = A%UPPER - A%LOWER
         T = INTERVAL(B,B)
         CALL RNDOUT_F90(T,.FALSE.,.TRUE.)
         B = T%UPPER
      END FUNCTION IWID_F90

      FUNCTION IMIG_F90(A) RESULT(B)
         IMPLICIT NONE
         TYPE(INTERVAL) :: A
         DOUBLE PRECISION B
         TYPE(INTERVAL) :: T   
         IF ( ((A%LOWER.GT.0D0) .AND. (A%UPPER.GT.0D0)) &
           .OR.((A%LOWER.LT.0D0) .AND. (A%UPPER.LT.0D0)) ) THEN
            T%LOWER = ABS(A%UPPER)
            CALL RNDOUT_F90(T,.TRUE.,.FALSE.)
            T%UPPER=T%LOWER
            T%LOWER = ABS(A%LOWER)
            CALL RNDOUT_F90(T,.TRUE.,.FALSE.)
            B = MIN (T%LOWER, T%UPPER)
         ELSE
            B = 0D0
         END IF
      END FUNCTION IMIG_F90

      FUNCTION IMID_F90(B)
         IMPLICIT NONE
         TYPE(INTERVAL):: B
         DOUBLE PRECISION IMID_F90      
         IMID_F90 = B%LOWER + (B%UPPER - B%LOWER) / 2D0
      END FUNCTION IMID_F90

 ! Additional routines for compatibility with Fortran-SC syntax
 
      FUNCTION SUP_F90(X) RESULT(R)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(IN) :: X
         DOUBLE PRECISION R
         R = X%UPPER
      END FUNCTION SUP_F90

      FUNCTION INF_F90(X) RESULT(R)
         IMPLICIT NONE
         TYPE(INTERVAL), INTENT(IN) :: X
         DOUBLE PRECISION R
         R = X%LOWER
      END FUNCTION INF_F90

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

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

!*** utilfuns.f
      SUBROUTINE ICAP (A, B, RESULT)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!
!   Hong Jiang
!   Dept. of Math., Stat. and Computer Science
!   Marquette University
!   Milwaukee, WI 53233
!
!         and
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  DICAP computes the intersection of the two intervals A and B, and
!  places  the result in RESULT.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION A(2), B(2), RESULT(2)
!
!***********************************************************************
!
! Internal variable declarations --
!
!
!***********************************************************************
!
! Common block declarations --
!
!
      INTEGER ISIG, IERR, IROUT, IPRTCL, ISEVER, IERPUN
      COMMON /ERFLGS/ ISIG, IERR, IROUT, IPRTCL, ISEVER, IERPUN
!
! The above common block stores signalling information for error
! conditions.  In particular,
!
! ISIG   is set to 0 at the beginning of SIMINI, and is reset
!        to the severity level of the error (1, 2, or 3) when an error
!        condition occurs.  The user may reset the flag to zero after an
!        error condition, depending on the error, if ISEVER is set to
!        allow execution after an error of its severity.
!
! IERR   is the number of the error condition, if an error occurred in
!        the last routine in INTLIB with error checking.  If no errors
!        occurred in the last such routine called, then IERR is zero.
!        The specific error conditions associated with particular error
!        numbers is defined in the routine ERRTST.
!
! IROUT  is the code number of the package routine in which the error
!        occurred
!
! IPRTCL controls the level of error which is printed.  IPRCTL=0 prints
!        all levels, while IPRCTL=1 prints only errors of level 1 or
!        greater.  If IPRTCL=4, then no error information is printed.
!
! ISEVER gives the error level which will stop execution.  ISEVER=0
!        causes any error to stop execution, while ISEVER>=3 causes only
!        errors of severity 3 to stop execution.  Generally, severity 3
!        errors correspond to when it is impossible to assign a result
!        with any meaningful interpretation.
!
! IERPUN is the Fortran unit number to which errors should be printed.
!
!***********************************************************************
!
! Fortran-supplied functions and subroutines --
!
!   MIN, MAX
!
!***********************************************************************
!
! Package-supplied functions and subroutines --
!
!   ERRTST
!
!***********************************************************************
!
! Internal variable declarations --
!
      DOUBLE PRECISION T(2)
!
!***********************************************************************
!
! Beginning of executable statements --
!
!  Identifying code for this routine --
!
      IROUT = 13
      IERR = 0
!
      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
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      LOGICAL FUNCTION IDISJ (A, B)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!
!   Hong Jiang
!   Dept. of Math., Stat. and Computer Science
!   Marquette University
!   Milwaukee, WI 53233
!
!         and
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  This function returns .TRUE. if the intervals A and B are disjoint,
!  and .FALSE. otherwise.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION A(2), B(2)
!
!***********************************************************************
!
! Beginning of executable statements --
!
      IDISJ = (A(2) .LT. B(1)) .OR. (B(2) .LT. A(1))
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      SUBROUTINE IHULL (A, B, RESULT)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!
!   Hong Jiang
!   Dept. of Math., Stat. and Computer Science
!   Marquette University
!   Milwaukee, WI 53233
!
!         and
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  IHULL returns the convex-hull of the interval A and the interval B
!  in RESULT.  If one of the intervals is empty (lower bound is greater
!  than upper bound), then just the upper interval is returned.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION A(2), B(2), RESULT(2)
!
!***********************************************************************
!
! Internal variable declarations --
!
      DOUBLE PRECISION T(2)
!
!***********************************************************************
!
! Beginning of executable statements --
!
      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
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      LOGICAL FUNCTION IILEI (A, B)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!
!   Hong Jiang
!   Dept. of Math., Stat. and Computer Science
!   Marquette University
!   Milwaukee, WI 53233
!
!         and
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  IILEI sets its value to .TRUE. if interval A is in the closure of
!  interval B.  The value of IILEI is .FALSE. otherwise.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION A(2), B(2)
!
!***********************************************************************
!
! Beginning of executable statements --
!
      IILEI = (A(1) .GE. B(1)) .AND. (A(2) .LE. B(2))
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      LOGICAL FUNCTION IILTI (A, B)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!
!   Hong Jiang
!   Dept. of Math., Stat. and Computer Science
!   Marquette University
!   Milwaukee, WI 53233
!
!         and
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  IILE sets its value to .TRUE. if interval A is in the interior of
!  interval B.  The value of IILE is .FALSE. otherwise.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION A(2), B(2)
!
!***********************************************************************
!
! Beginning of executable statements --
!
      IILTI = (A(1) .GT. B(1)) .AND. (A(2) .LT. B(2))
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      DOUBLE PRECISION FUNCTION IINF (A)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!
!   Hong Jiang
!   Dept. of Math., Stat. and Computer Science
!   Marquette University
!   Milwaukee, WI 53233
!
!         and
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  IINF returns the lower endpoint of the interval A.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION A(2)
!
!***********************************************************************
!
! Beginning of executable statements --
!
      IINF = A(1)
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      DOUBLE PRECISION FUNCTION IMID (X)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!
!   Hong Jiang
!   Dept. of Math., Stat. and Computer Science
!   Marquette University
!   Milwaukee, WI 53233
!
!         and
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  IMID returns a floating point approximation to the midpoint
!  of the interval A, using available floating point arithmetic.  The
!  value returned by this routine can be considered to DEFINE the
!  midpoint.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION X(2)
!
!***********************************************************************
!
! Common block declarations --
!
      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)
!
      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
!
      INTEGER ITINY2, JTINY2
!
      COMMON /IMATH/   ITINY2, JTINY2
!
!
!  The above common blocks hold mathematical constants which are used in
!  the elementary function routines.
!
!  Variable descriptions
!
!  MAXX        a double precision representation of the largest
!              representable integer
!
!  MXLGM1      an approximation to the logarithm of .25 times the
!              largest representable machine number, minus 1. This
!              should be a rigorous lower bound on the logarithm of the
!              largest representable machine number.  Its default
!              computation in SIMINI is to use the Fortran LOG function.
!              This should be changed if the Fortran LOG function is not
!              sufficiently accurate.
!
!
!  A           an interval enclosure for 1/pi
!
!  PI          an interval enclosure for pi
!
!  PI2         an interval enclosure for pi/2
!
!  PI3         an interval enclosure for pi/3
!
!  PI4         an interval enclosure for pi/4
!
!  PI6         an interval enclosure for pi/6
!
!  PI8         an interval enclosure for pi/8
!
!  E           an interval enclosure for E
!
!  E14         an interval enclosure for e^{1/4}
!
!  ESXTNT      an interval enclosure for e^{1/16}
!
!  SQT10       an interval enclosure for SQRT(10)
!
!  TINY2       the maximum of the smallest machine number and the
!              reciprocal of the largest machine number.  This quantity
!              is checked to avoid overflow in certain places.
!
!  CBTEP       an approximation to the cube root of six times the
!              largest distance between numbers.
!
!  SQT3        an interval enclosure for the square root of 3.
!
!  ODSQT3      the reciprocal of the square root of 3 times
!              1+ 100*MXULP, used in argument reduction in the
!              arctangent routine.
!
!  ITINY2      the logarithm base 10 of TINY2 truncated to an integer.
!
!  JTINY2      sixteen times the logarithm base e of TINY2 truncated to
!              an integer.
!
!  See the statements in the routine SIMINI for the definitions of
!  the other constants.
!
!***********************************************************************
!
! Beginning of executable statements --
!
      IMID = X(1) + (X(2) - X(1)) / TWO(1)
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      DOUBLE PRECISION FUNCTION IMIG (A)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  IMIG returns the mignitude of the interval A. Since ABS is not
!  assumed to give an exact result, the result is rounded down.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION A(2)
!
!***********************************************************************
!
! Internal variable declarations --
!
      DOUBLE PRECISION TMP(2)
!
!***********************************************************************
!
! Fortran-supplied functions and subroutines --
!
!  ABS, MIN
!
!***********************************************************************
!
! Package-supplied functions and subroutines --
!
!  RNDOUT
!
!***********************************************************************
!
! Beginning of executable statements --
!
      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
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      SUBROUTINE INEG (A, RESULT)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!
!   Hong Jiang
!   Dept. of Math., Stat. and Computer Science
!   Marquette University
!   Milwaukee, WI 53233
!
!         and
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  INEG performs interval unary negation.  The result is rounded out in
!  case the negatives of the endpoints are not representable.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION A(2), RESULT(2)
!
!***********************************************************************
!
! Package-supplied functions and subroutines --
!
!  RNDOUT
!
!***********************************************************************
!
! Internal variable declarations --
!
      DOUBLE PRECISION T(2)
!
!***********************************************************************
!
! Beginning of executable statements --
!
      T(1) = A(1)
      T(2) = A(2)
      RESULT(1) = -T(2)
      RESULT(2) = -T(1)
      CALL RNDOUT (RESULT, .TRUE.,.TRUE.)
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      DOUBLE PRECISION FUNCTION INTABS (A)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!
!   Hong Jiang
!   Dept. of Math., Stat. and Computer Science
!   Marquette University
!   Milwaukee, WI 53233
!
!         and
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  INTABS returns the maximum absolute value of the interval A.
!  Since ABS is not assumed to give an exact result, the result is
!  rounded up.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION A(2)
!
!***********************************************************************
!
! Internal variable declarations --
!
      DOUBLE PRECISION TMP(2)
!
!***********************************************************************
!
! Fortran-supplied functions and subroutines --
!
!  ABS, MAX
!
!***********************************************************************
!
! Package-supplied functions and subroutines --
!
!  RNDOUT
!
!***********************************************************************
!
! Beginning of executable statements --
!
      TMP(2) = ABS(A(1))
      CALL RNDOUT(TMP,.FALSE.,.TRUE.)
      TMP(1)=TMP(2)
      TMP(2) = ABS(A(2))
      CALL RNDOUT(TMP,.FALSE.,.TRUE.)
!
      INTABS = MAX (TMP(1), TMP(2))
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      LOGICAL FUNCTION IRLEI (A, B)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!
!   Hong Jiang
!   Dept. of Math., Stat. and Computer Science
!   Marquette University
!   Milwaukee, WI 53233
!
!         and
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  IRLEI sets its value to .TRUE. if double precision value A is in the
!  closure of the interval B. The value of IRLEI is set to .FALSE.
!  otherwise.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION A, B(2)
!
!***********************************************************************
!
! Beginning of executable statements --
!
      IRLEI = (A .GE. B(1)) .AND. (A .LE. B(2))
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      LOGICAL FUNCTION IRLTI (A, B)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!
!   Hong Jiang
!   Dept. of Math., Stat. and Computer Science
!   Marquette University
!   Milwaukee, WI 53233
!
!         and
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  IRLEI sets its value to .TRUE. if double precision value A is in the
!  interior of the interval B. The value of IRLEI is set to .FALSE.
!  otherwise.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION A, B(2)
!
!***********************************************************************
!
! Beginning of executable statements --
!
      IRLTI = (A .GT. B(1)) .AND. (A .LT. B(2))
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      DOUBLE PRECISION FUNCTION ISUP (A)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!
!   Hong Jiang
!   Dept. of Math., Stat. and Computer Science
!   Marquette University
!   Milwaukee, WI 53233
!
!         and
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  IINF returns the upper endpoint of the interval A.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION A(2)
!
!***********************************************************************
!
! Beginning of executable statements --
!
      ISUP = A(2)
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      SUBROUTINE IVL1 (R, RESULT)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!
!   Hong Jiang
!   Dept. of Math., Stat. and Computer Science
!   Marquette University
!   Milwaukee, WI 53233
!
!         and
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  IVL1 constructs an interval RESULT from the double precision variable
!  R.  This is done by using simulated directed roundings with RNDOUT.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION R, RESULT(2)
!
!***********************************************************************
!
! Package-supplied functions and subroutines --
!
!  RNDOUT
!
!***********************************************************************
!
! Internal variable declarations --
!
      DOUBLE PRECISION T
!
!***********************************************************************
!
! Beginning of executable statements --
!
      T = R
      RESULT(1) = T
      RESULT(2) = T
      IF (T.NE.0D0) CALL RNDOUT (RESULT, .TRUE.,.TRUE.)
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      SUBROUTINE IVL2 (R1, R2, RESULT)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!
!   Hong Jiang
!   Dept. of Math., Stat. and Computer Science
!   Marquette University
!   Milwaukee, WI 53233
!
!         and
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  IVL2 constructs an interval RESULT, roughly equal to [R1,R2], from
!  the double precision variables R1 and R2.  This is done by using
!  simulated directed roundings with RNDOUT.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION R1, R2, RESULT(2)
!
!***********************************************************************
!
! Package-supplied functions and subroutines --
!
!  RNDOUT
!
!***********************************************************************
!
! Internal variable declarations --
!
      DOUBLE PRECISION T(2)
      LOGICAL RNDUP, RNDDWN
!
!***********************************************************************
!
! Common block declarations --
!
!
      INTEGER ISIG, IERR, IROUT, IPRTCL, ISEVER, IERPUN
      COMMON /ERFLGS/ ISIG, IERR, IROUT, IPRTCL, ISEVER, IERPUN
!
! The above common block stores signalling information for error
! conditions.  In particular,
!
! ISIG   is set to 0 at the beginning of SIMINI, and is reset
!        to the severity level of the error (1, 2, or 3) when an error
!        condition occurs.  The user may reset the flag to zero after an
!        error condition, depending on the error, if ISEVER is set to
!        allow execution after an error of its severity.
!
! IERR   is the number of the error condition, if an error occurred in
!        the last routine in INTLIB with error checking.  If no errors
!        occurred in the last such routine called, then IERR is zero.
!        The specific error conditions associated with particular error
!        numbers is defined in the routine ERRTST.
!
! IROUT  is the code number of the package routine in which the error
!        occurred
!
! IPRTCL controls the level of error which is printed.  IPRCTL=0 prints
!        all levels, while IPRCTL=1 prints only errors of level 1 or
!        greater.  If IPRTCL=4, then no error information is printed.
!
! ISEVER gives the error level which will stop execution.  ISEVER=0
!        causes any error to stop execution, while ISEVER>=3 causes only
!        errors of severity 3 to stop execution.  Generally, severity 3
!        errors correspond to when it is impossible to assign a result
!        with any meaningful interpretation.
!
! IERPUN is the Fortran unit number to which errors should be printed.
!
!***********************************************************************
!
! Beginning of executable statements --
!
!  Identifying code for this routine --
!
      IROUT = 15
      IERR = 0
!
      T(1) = R1
      T(2) = R2
      RESULT(1) = T(1)
      RESULT(2) = T(2)
!
      IF (RESULT(2).LT.RESULT(1)) THEN
         IERR = 1
         CALL ERRTST(RESULT)
      END IF
!
      RNDDWN = T(1).NE.0D0
      RNDUP = T(2).NE.0D0
      CALL RNDOUT (RESULT, RNDDWN, RNDUP)
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      SUBROUTINE IVLABS (A, RESULT)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  IVLABS returns a rigorous bound on the range of the absolute value
!  of the interval A.  The intrinsic function ABS is assumed to be
!  accurate to the same accuracy as the elementary operations.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION A(2), RESULT(2)
!
!***********************************************************************
!
! Internal variable declarations --
!
      DOUBLE PRECISION TMP(2), TA(2)
!
!***********************************************************************
!
! Fortran-supplied functions and subroutines --
!
!  ABS
!
!***********************************************************************
!
! Package-supplied functions and subroutines --
!
!  RNDOUT
!
!***********************************************************************
!
! Beginning of executable statements --
!
      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)        &
     &   RESULT(1) = 0D0

      RETURN
      END
!***********************************************************************
!***********************************************************************
      SUBROUTINE IVLI (A, RESULT)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   January 8, 1993
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  IVLI places the contents of interval A in RESULT.
!
!***********************************************************************
!
! Argument declarations --
!
      DOUBLE PRECISION A(2), RESULT(2)
!
!***********************************************************************
!
! Beginning of executable statements --
!
      RESULT(1) = A(1)
      RESULT(2) = A(2)
!
      RETURN
      END
!***********************************************************************
!***********************************************************************
      DOUBLE PRECISION FUNCTION IWID (A)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
! Written by:
!
!   R. Baker Kearfott
!
!   Department of Mathematics
!   U.S.L. Box 4-1010
!   Lafayette, LA 70504
!
!   April 11, 1992
!
! Part of the interval elementary function library
! (Utility function)
!
!***********************************************************************
!
! Function --
!
!  IWID returns the width of the interval A, rounded up.
!
!***********************************************