C     ALGORITHM 619, COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 3,
C     SEP., 1984, P. 348-353.
C   DRIVER PROGRAM FOR DLAINV
C   THE INVERSE LAPLACE TRANSFORM OF 1/(P**2+1) IS COMPUTED
C   FOR T=0.1,1,2,3,4,5,10,20,30,40,50,60,70,80,90 AND 100
C   (THESE VALUES ARE STORED IN THE ARRAY TVAL).
C   THE REQUESTED TOLERANCES ARE EPSAB=EPSRE=1.0D-4, 1.0D-8
C   AND 1.0D-12 (THESE VALUES ARE STORED IN THE ARRAY E)
C   THE EXACT INVERSE LAPLACE TRANSFORM IS SIN(T).
C   ALSO THE EXACT ERROR IS COMPUTED
C
      DOUBLE PRECISION C,DIF,E,ESTERR,EXA,RESULT,T,TVAL
      DOUBLE PRECISION DABS,DSIN
      DIMENSION E(3),TVAL(16)
      EXTERNAL FUN
C
C   THE ARRAY E CONTAINS THE REQUESTED TOLERANCES
C
      DATA E(1),E(2),E(3)/1.0D-4,1.0D-8,1.0D-12/
C
C   THE ARRAY TVAL CONTAINS THE VALUES OF THE INDEPENDENT
C   VARIABLE OF THE INVERSE LAPLACE TRANSFORM
C
      DATA TVAL(1),TVAL(2),TVAL(3),TVAL(4),TVAL(5),TVAL(6),
     *  TVAL(7),TVAL(8),TVAL(9),TVAL(10),TVAL(11),TVAL(12),
     *  TVAL(13),TVAL(14),TVAL(15),TVAL(16) /0.1D0,1.0D0,2.0D0,
     *  3.0D0,4.0D0,5.0D0,1.0D1,2.0D1,3.0D1,4.0D1,5.0D1,6.0D1,
     *  7.0D1,8.0D1,9.0D1,1.0D2/
C
C   C IS THE ABSCISSA OF CONVERGENCE
C
      C = 0.0D0
      WRITE(6,99997)
      DO 20 I=1,3
        WRITE(6,99998) E(I)
        DO 10 K=1,16
          T = TVAL(K)
C
C   COMPUTATION OF THE APPROXIMATE INVERSE LAPLACE TRANSFORM
C
          CALL DLAINV(FUN,T,C,E(I),E(I),RESULT,ESTERR,NUM,IER)
C
C   COMPUTATION OF THE EXACT INVERSE LAPLACE TRANSFORM
C
          EXA = DSIN(T)
C
C   COMPUTATION OF THE EXACT ERROR
C
          DIF = DABS(EXA-RESULT)
          WRITE(6,99999) T,EXA,RESULT,DIF,ESTERR,NUM,IER
   10   CONTINUE
   20 CONTINUE
99997 FORMAT(1H0,4X,1HT,5X,1H*,8X,11HEXACT VALUE,13X,8HCOMPUTED,
     *  6H VALUE,9X,5HERROR,7X,6HESTERR,4X,3HNUM,4X,3HIER//1H ,
     *  98(1H*)/)
99998 FORMAT(1H0,16HEPSAB = EPSRE = ,D9.2/)
99999 FORMAT(5X,F5.1,1X,1H*,1X,2(D23.16,3X),2(D9.2,3X),I3,I6)
      STOP
      END
       SUBROUTINE DLAINV(FUN,T,C,EPSRE,EPSAB,RESULT,ESTERR,NUM,
     *   IER)
C
C ......................................................................
C
C   1. DLAINV
C        INVERSION OF LAPLACE TRANSFORM USING THE DURBIN FORMULA
C        IN COMBINATION WITH THE EPSILON ALGORITHM
C
C   2. PURPOSE
C        THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO THE
C        INVERSE LAPLACE TRANSFORM F(T) OF FUN, FOR THE VALUE
C        OF THE INDEPENDENT VARIABLE EQUAL TO T, HOPEFULLY
C        SATISFYING FOLLOWING CLAIM FOR ACCURACY :
C        ABS(F(T)-RESULT).LE.MAX(EPSAB,EPSR*ABS(F(T)))
C
C   3. CALLING SEQUENCE
C        CALL DLAINV(FUN,T,C,EPSRE,EPSAB,RESULT,ESTERR,NUM,IER)
C
C      INPUT PARAMETERS
C        FUN    - DOUBLE PRECISION
C                 SUBROUTINE DEFINING THE LAPLACE TRANSFORM AS
C                 A COMPLEX FUNCTION.  THE CALLING SEQUENCE OF
C                 FUN IS : CALL FUN(A,B,C,D) WHERE
C                 A : DOUBLE PRECISION
C                     REAL PART OF THE INDEPENDENT VARIABLE
C                     OF THE LAPLACE TRANSFORM (INPUT)
C                 B : DOUBLE PRECISION
C                     IMAGINARY PART OF THE INDEPENDENT
C                     VARIABLE OF THE LAPLACE TRANSFORM (INPUT)
C                 C : DOUBLE PRECISION
C                     REAL PART OF THE VALUE OF THE LAPLACE
C                     TRANSFORM (OUTPUT)
C                 D : DOUBLE PRECISION
C                     IMAGINARY PART OF THE VALUE OF THE
C                     LAPLACE TRANSFORM (OUTPUT)
C                 THE ACTUAL NAME FOR FUN NEEDS TO BE DECLARED
C                 EXTERNAL IN THE DRIVER PROGRAM.
C
C        T      - DOUBLE PRECISION
C                 VALUE OF THE INDEPENDENT VARIABLE FOR WHICH THE
C                 INVERSE LAPLACE TRANSFORM HAS TO BE COMPUTED.
C                 T SHOULD BE GREATER THAN ZERO.
C
C        C      - DOUBLE PRECISION
C                 ABSCISSA OF CONVERGENCE OF THE LAPLACE TRANSFORM
C
C        EPSRE  - DOUBLE PRECISION
C                 RELATIVE ACCURACY REQUESTED
C
C        EPSAB  - DOUBLE PRECISION
C                 ABSOLUTE ACCURACY REQUESTED.
C                 THE ROUTINE TRIES TO SATISFY THE LEAST STRINGENT
C                 OF BOTH ACCURACY REQUIREMENT.
C
C      OUTPUT PARAMETERS
C        RESULT - DOUBLE PRECISION
C                 INVERSE LAPLACE TRANSFORM
C
C        ESTERR - DOUBLE PRECISION
C                 ESTIMATE OF THE ABSOLUTE ERROR ABS(F(T)-RESULT)
C
C        NUM    - INTEGER
C                 NUMBER OF EVALUATIONS OF FUN
C
C        IER    - INTEGER
C                 PARAMETER GIVING INFORMATION ON THE TERMINATION
C                 OF THE ALGORITHM
C                 IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
C                         ROUTINE
C                 IER = 1 THE COMPUTATIONS ARE TERMINATED BECAUSE
C                         THE BOUND ON THE NUMBER OF EVALUATIONS
C                         OF FUN HAS BEEN ACHIEVED.  THIS BOUND
C                         IS EQUAL TO 8*MAX+5 WHERE  MAX  IS A
C                         NUMBER INITIALIZED IN A DATA
C                         STATEMENT.  ONE CAN ALLOW MORE FUNCTION
C                         EVALUATIONS BY INCREASING THE VALUE OF
C                         MAX IN THE DATA-STATEMENT.
C                 IER = 2 THE VALUE OF T IS LESS THAN OR EQUAL
C                         TO ZERO.
C
C   4. SUBROUTINES OR FUNCTIONS NEEDED
C               - FUN : USER PROVIDED SUBROUTINE
C               - DQEXT : EPSILON ALGORITHM
C               - D1MACH : THIS SUBPROGRAM IS CALLED BY
C                          DQEXT, AND PROVIDES MACHINE CONSTANTS
C               - DABS, DATAN, DEXP, DMAX1, DSIN : FORTRAN FUNCTIONS
C
C ......................................................................
C
      DOUBLE PRECISION AIM,AK,ARE,ARG,BB,C,DABS,DATAN,DEXP,DMAX1,
     *  DSIN,EPSAB,EPSRE,ESTERR,FIM,FRE,PID16,R,RESULT,RES3LA,REX,
     *  SI,T
      INTEGER I,IER,K,KC,KK,KS,M,NEX,NRES,NUM
      DIMENSION SI(32),RES3LA(3),REX(52)
C
C   THE ARRAY SI CONTAINS VALUES OF THE SINE AND COSINE FUNCTIONS
C   REQUIRED IN THE DURBIN FORMULA. SI(8) AND SI(16) ARE GIVEN IN
C   THE FOLLOWING DATA STATEMENT. THE OTHER VALUES ARE COMPUTED.
C
      DATA SI(8),SI(16)/ 1.0D+00,0.0D+00/
C
C   MAX IS A BOUND ON THE NUMBER OF TERMS USED IN THE DURBIN
C   FORMULA.
C
      DATA MAX/500/
C
C   TEST ON VALIDITY OF THE INPUT PARAMETER T
C
      IER = 2
      RESULT = 0.0D+00
      ESTERR = 1.0D+00
      NUM = 0
      IF (T.LE.0.0D+00) GO TO 999
C
C   PID16 IS EQUAL TO PI/16
C
      PID16 = DATAN(1.0D+00)/4.0D+00
C
C   COMPUTATION OF THE ELEMENTS OF THE ARRAY SI
C
      AK = 1.0D+00
      DO 10 K=1,7
        SI(K) = DSIN(AK*PID16)
        AK = AK+1.0D+00
        KK = 16-K
        SI(KK) = SI(K)
10    CONTINUE
      IER = 0
      NRES = 0
      DO 20 K=17,32
        SI(K) = -SI(K-16)
20    CONTINUE
C
C   INITIALIZATION OF THE SUMMATION OF THE DURBIN FORMULA.
C
      ARG = PID16/T
      ARE = C+2.0D+00/T
      AIM = 0.0D+00
      BB = DEXP(ARE*T)/(1.6D+01*T)
      CALL FUN (ARE,AIM,FRE,FIM)
      NUM = 5
      R = 5.0D-01*FRE
      NEX = 0
      KC = 8
      KS = 0
C
C   MAIN LOOP FOR THE SUMMATION
C
      DO 40 I=1,MAX
        M = 8
        IF (I.EQ.1) M = 12
        DO 30 K=1,M
          AIM = AIM+ARG
          KC = KC+1
          KS = KS+1
          IF (KC.GT.32) KC = 1
          IF (KS.GT.32) KS = 1
          CALL FUN(ARE,AIM,FRE,FIM)
          R = R+FRE*SI(KC)-FIM*SI(KS)
30      CONTINUE
        NUM = NUM+8
        NEX = NEX+1
        REX(NEX) = R
C
C   EXTRAPOLATION USING THE EPSILON ALGORITHM
C
        IF(NEX.GE.3) CALL DQEXT(NEX,REX,RESULT,ESTERR,RES3LA,NRES)
        IF(NRES.LT.4) GO TO 40
C
C   COMPUTATION OF INTERMEDIATE RESULT AND ESTIMATE OF THE
C   ABSOLUTE ERROR
C
        RESULT = RESULT * BB
        ESTERR = ESTERR * BB
        IF (ESTERR.LT.DMAX1(EPSAB,EPSRE*DABS(RESULT)).AND.DABS(R*BB-
     *   RESULT).LT.5.0D-01*DABS(RESULT)) GO TO 999
40    CONTINUE
C
C   SET ERROR FLAG IN THE CASE THAT THE NUMBER OF TERMS IN THE
C   SUMMATION IS EQUAL TO MAX
C
      IER = 1
999   RETURN
      END
      SUBROUTINE DQEXT(N,EPSTAB,RESULT,ABSERR,RES3LA,NRES)
C
C 1.        DQEXT
C           EPSILON ALGORITHM
C              STANDARD FORTRAN SUBROUTINE
C              DOUBLE PRECISION VERSION
C
C 2.        PURPOSE
C              THE ROUTINE DETERMINES THE LIMIT OF A GIVEN SEQUENCE OF
C              APPROXIMATIONS, BY MEANS OF THE EPSILON ALGORITHM
C              OF P. WYNN.
C              AN ESTIMATE OF THE ABSOLUTE ERROR IS ALSO GIVEN.
C              THE CONDENSED EPSILON TABLE IS COMPUTED. ONLY THOSE
C              ELEMENTS NEEDED FOR THE COMPUTATION OF THE NEXT DIAGONAL
C              ARE PRESERVED.
C
C 3.        CALLING SEQUENCE
C              CALL DQEXT(N,EPSTAB,RESULT,ABSERR,RES3LA,NRES)
C
C           PARAMETERS
C              N      - INTEGER
C                       EPSTAB(N) CONTAINS THE NEW ELEMENT IN THE
C                       FIRST COLUMN OF THE EPSILON TABLE.
C
C              EPSTAB - DOUBLE PRECISION
C                       VECTOR OF DIMENSION 52 CONTAINING THE ELEMENTS
C                       OF THE TWO LOWER DIAGONALS OF THE
C                       TRIANGULAR EPSILON TABLE
C                       THE ELEMENTS ARE NUMBERED STARTING AT THE
C                       RIGHT-HAND CORNER OF THE TRIANGLE.
C
C              RESULT - DOUBLE PRECISION
C                       RESULTING APPROXIMATION TO THE INTEGRAL
C
C              ABSERR - DOUBLE PRECISION
C                       ESTIMATE OF THE ABSOLUTE ERROR COMPUTED FROM
C                       RESULT AND THE 3 PREVIOUS RESULTS
C
C              RES3LA - DOUBLE PRECISION
C                       VECTOR OF DIMENSION 3 CONTAINING THE LAST 3
C                       RESULTS
C
C              NRES   - INTEGER
C                       NUMBER OF CALLS TO THE ROUTINE
C                       (SHOULD BE ZERO AT FIRST CALL)
C
C 4.        SUBROUTINES OR FUNCTIONS NEEDED
C                     - D1MACH
C                     - FORTRAN DABS, DMAX1
C
C     ..................................................................
C
      DOUBLE PRECISION ABSERR,DABS,DELTA1,DELTA2,DELTA3,DMAX1,EPMACH,
     *  EPSINF,EPSTAB,ERROR,ERR1,ERR2,ERR3,E0,E1,E1ABS,E2,E3,OFLOW,
     *  RES,RESULT,RES3LA,SS,TOL1,TOL2,TOL3
      INTEGER I,IB,IB2,IE,INDX,K1,K2,K3,LIMEXP,N,NEWELM,NRES,NUM
      DIMENSION EPSTAB(52),RES3LA(3)
C
C           LIST OF MAJOR VARIABLES
C           -----------------------
C
C           E0     - THE 4 ELEMENTS ON WHICH THE
C           E1       COMPUTATION OF A NEW ELEMENT IN
C           E2       THE EPSILON TABLE IS BASED
C           E3                 E0
C                        E3    E1    NEW
C                              E2
C           NEWELM - NUMBER OF ELEMENTS TO BE COMPUTED IN THE NEW
C                    DIAGONAL
C           ERROR  - ERROR = ABS(E1-E0)+ABS(E2-E1)+ABS(NEW-E2)
C           RESULT - THE ELEMENT IN THE NEW DIAGONAL WITH LEAST VALUE
C                    OF ERROR
C
C           MACHINE DEPENDENT CONSTANTS
C           ---------------------------
C
C           EPMACH IS THE LARGEST RELATIVE SPACING.
C           OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
C           LIMEXP IS THE MAXIMUM NUMBER OF ELEMENTS THE EPSILON
C           TABLE CAN CONTAIN. IF THIS NUMBER IS REACHED, THE UPPER
C           DIAGONAL OF THE EPSILON TABLE IS DELETED.
C
C***FIRST EXECUTABLE STATEMENT
      OFLOW = D1MACH(2)
      EPMACH = D1MACH(4)
      NRES = NRES+1
      ABSERR = OFLOW
      RESULT = EPSTAB(N)
      IF(N.LT.3) GO TO 100
      LIMEXP = 50
      EPSTAB(N+2) = EPSTAB(N)
      NEWELM = (N-1)/2
      EPSTAB(N) = OFLOW
      NUM = N
      K1 = N
      DO 40 I = 1,NEWELM
        K2 = K1-1
        K3 = K1-2
        RES = EPSTAB(K1+2)
        E0 = EPSTAB(K3)
        E1 = EPSTAB(K2)
        E2 = RES
        E1ABS = DABS(E1)
        DELTA2 = E2-E1
        ERR2 = DABS(DELTA2)
        TOL2 = DMAX1(DABS(E2),E1ABS)*EPMACH
        DELTA3 = E1-E0
        ERR3 = DABS(DELTA3)
        TOL3 = DMAX1(E1ABS,DABS(E0))*EPMACH
        IF(ERR2.GT.TOL2.OR.ERR3.GT.TOL3) GO TO 10
C
C           IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE
C           ACCURACY, CONVERGENCE IS ASSUMED.
C           RESULT = E2
C           ABSERR = ABS(E1-E0)+ABS(E2-E1)
C
        RESULT = RES
        ABSERR = ERR2+ERR3
C***JUMP OUT OF DO-LOOP
        GO TO 100
   10   E3 = EPSTAB(K1)
        EPSTAB(K1) = E1
        DELTA1 = E1-E3
        ERR1 = DABS(DELTA1)
        TOL1 = DMAX1(E1ABS,DABS(E3))*EPMACH
C
C           IF TWO ELEMENTS ARE VERY CLOSE TO EACH OTHER, OMIT
C           A PART OF THE TABLE BY ADJUSTING THE VALUE OF N
C
        IF(ERR1.LE.TOL1.OR.ERR2.LE.TOL2.OR.ERR3.LE.TOL3) GO TO 20
        SS = 1.0D+00/DELTA1+1.0D+00/DELTA2-1.0D+00/DELTA3
        EPSINF = DABS(SS*E1)
C
C           TEST TO DETECT IRREGULAR BEHAVIOUR IN THE TABLE, AND
C           EVENTUALLY OMIT A PART OF THE TABLE ADJUSTING THE VALUE
C           OF N.
C
        IF(EPSINF.GT.1.0D-04) GO TO 30
   20   N = I+I-1
C***JUMP OUT OF DO-LOOP
        GO TO 50
C
C           COMPUTE A NEW ELEMENT AND EVENTUALLY ADJUST
C           THE VALUE OF RESULT.
C
   30   RES = E1+1.0D+00/SS
        EPSTAB(K1) = RES
        K1 = K1-2
        ERROR = ERR2+DABS(RES-E2)+ERR3
        IF(ERROR.GT.ABSERR) GO TO 40
        ABSERR = ERROR
        RESULT = RES
   40 CONTINUE
C
C           SHIFT THE TABLE.
C
   50 IF(N.EQ.LIMEXP) N = 2*(LIMEXP/2)-1
      IB = 1
      IF((NUM/2)*2.EQ.NUM) IB = 2
      IE = NEWELM+1
      DO 60 I=1,IE
        IB2 = IB+2
        EPSTAB(IB) = EPSTAB(IB2)
        IB = IB2
   60 CONTINUE
      IF(NUM.EQ.N) GO TO 80
      INDX = NUM-N+1
      DO 70 I = 1,N
        EPSTAB(I)= EPSTAB(INDX)
        INDX = INDX+1
   70 CONTINUE
   80 IF(NRES.GE.4) GO TO 90
      RES3LA(NRES) = RESULT
      ABSERR = OFLOW
      GO TO 100
C
C           COMPUTE ERROR ESTIMATE
C
   90 ABSERR = DABS(RESULT-RES3LA(3))+DABS(RESULT-RES3LA(2))
     *  +DABS(RESULT-RES3LA(1))
      RES3LA(1) = RES3LA(2)
      RES3LA(2) = RES3LA(3)
      RES3LA(3) = RESULT
  100 ABSERR = DMAX1(ABSERR,5.0D+00*EPMACH*DABS(RESULT))
      RETURN
      END
      SUBROUTINE FUN(X,Y,A,B)
C
C   IN THIS SUBROUTINE THE FUNCTION 1/(P**2+1) IS EVALUATED
C   WHERE P=X+IY. THIS FUNCTION IS THE LAPLACE TRANSFORM OF
C   SIN(T)
C
      DOUBLE PRECISION X,Y,A,B,C,D
C  A+IB = 1/((X+IY)**2+1)
      C = X*X-Y*Y+1.0D0
      D = C*C + 4.0D0*X*X*Y*Y
      A = C/D
      B = -2.0D0*X*Y/D
      RETURN
      END
      DOUBLE PRECISION FUNCTION D1MACH(I)
C***BEGIN PROLOGUE  D1MACH
C***DATE WRITTEN   750101   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  R1
C***KEYWORDS  MACHINE CONSTANTS
C***AUTHOR  FOX, P. A., (BELL LABS)
C           HALL, A. D., (BELL LABS)
C           SCHRYER, N. L., (BELL LABS)
C***PURPOSE  RETURNS DOUBLE PRECISION MACHINE DEPENDENT CONSTANTS
C***DESCRIPTION
C     D1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS
C     FOR THE LOCAL MACHINE ENVIRONMENT.  IT IS A FUNCTION
C     SUBPROGRAM WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED
C     AS FOLLOWS, FOR EXAMPLE
C
C          D = D1MACH(I)
C
C     WHERE I=1,...,5.  THE (OUTPUT) VALUE OF D ABOVE IS
C     DETERMINED BY THE (INPUT) VALUE OF I.  THE RESULTS FOR
C     VARIOUS VALUES OF I ARE DISCUSSED BELOW.
C
C  DOUBLE-PRECISION MACHINE CONSTANTS
C  D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
C  D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
C  D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
C  D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
C  D1MACH( 5) = LOG10(B)
C***REFERENCES  FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A
C                 PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL
C                 SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188.
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  D1MACH
C
      INTEGER SMALL(4)
      INTEGER LARGE(4)
      INTEGER RIGHT(4)
      INTEGER DIVER(4)
      INTEGER LOG10(4)
C
      DOUBLE PRECISION DMACH(5)
C
      EQUIVALENCE (DMACH(1),SMALL(1))
      EQUIVALENCE (DMACH(2),LARGE(1))
      EQUIVALENCE (DMACH(3),RIGHT(1))
      EQUIVALENCE (DMACH(4),DIVER(1))
      EQUIVALENCE (DMACH(5),LOG10(1))
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C     DATA SMALL(1) / ZC00800000 /
C     DATA SMALL(2) / Z000000000 /
C
C     DATA LARGE(1) / ZDFFFFFFFF /
C     DATA LARGE(2) / ZFFFFFFFFF /
C
C     DATA RIGHT(1) / ZCC5800000 /
C     DATA RIGHT(2) / Z000000000 /
C
C     DATA DIVER(1) / ZCC6800000 /
C     DATA DIVER(2) / Z000000000 /
C
C     DATA LOG10(1) / ZD00E730E7 /
C     DATA LOG10(2) / ZC77800DC0 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
C
C     DATA SMALL(1) / O1771000000000000 /
C     DATA SMALL(2) / O0000000000000000 /
C
C     DATA LARGE(1) / O0777777777777777 /
C     DATA LARGE(2) / O0007777777777777 /
C
C     DATA RIGHT(1) / O1461000000000000 /
C     DATA RIGHT(2) / O0000000000000000 /
C
C     DATA DIVER(1) / O1451000000000000 /
C     DATA DIVER(2) / O0000000000000000 /
C
C     DATA LOG10(1) / O1157163034761674 /
C     DATA LOG10(2) / O0006677466732724 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
C
C     DATA SMALL(1) / O1771000000000000 /
C     DATA SMALL(2) / O7770000000000000 /
C
C     DATA LARGE(1) / O0777777777777777 /
C     DATA LARGE(2) / O7777777777777777 /
C
C     DATA RIGHT(1) / O1461000000000000 /
C     DATA RIGHT(2) / O0000000000000000 /
C
C     DATA DIVER(1) / O1451000000000000 /
C     DATA DIVER(2) / O0000000000000000 /
C
C     DATA LOG10(1) / O1157163034761674 /
C     DATA LOG10(2) / O0006677466732724 /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C
C     DATA SMALL(1) / 00604000000000000000B /
C     DATA SMALL(2) / 00000000000000000000B /
C
C     DATA LARGE(1) / 37767777777777777777B /
C     DATA LARGE(2) / 37167777777777777777B /
C
C     DATA RIGHT(1) / 15604000000000000000B /
C     DATA RIGHT(2) / 15000000000000000000B /
C
C     DATA DIVER(1) / 15614000000000000000B /
C     DATA DIVER(2) / 15010000000000000000B /
C
C     DATA LOG10(1) / 17164642023241175717B /
C     DATA LOG10(2) / 16367571421742254654B /
C
C     MACHINE CONSTANTS FOR THE CRAY 1
C
C     DATA SMALL(1) / 200004000000000000000B /
C     DATA SMALL(2) / 00000000000000000000B /
C
C     DATA LARGE(1) / 577777777777777777777B /
C     DATA LARGE(2) / 000007777777777777777B /
C
C     DATA RIGHT(1) / 377214000000000000000B /
C     DATA RIGHT(2) / 000000000000000000000B /
C
C     DATA DIVER(1) / 377224000000000000000B /
C     DATA DIVER(2) / 000000000000000000000B /
C
C     DATA LOG10(1) / 377774642023241175717B /
C     DATA LOG10(2) / 000007571421742254654B /
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD -
C     STATIC DMACH(5)
C
C     DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/
C     DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/
C     DATA LOG10/40423K,42023K,50237K,74776K/
C
C     MACHINE CONSTANTS FOR THE HARRIS 220
C
C     DATA SMALL(1),SMALL(2) / [20000000, [00000201 /
C     DATA LARGE(1),LARGE(2) / [37777777, [37777577 /
C     DATA RIGHT(1),RIGHT(2) / [20000000, [00000333 /
C     DATA DIVER(1),DIVER(2) / [20000000, [00000334 /
C     DATA LOG10(1),LOG10(2) / [23210115, [10237777 /
C
C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.
C
C     DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 /
C     DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 /
C     DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 /
C     DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 /
C     DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /
C
C      MACHINE CONSTANTS FOR THE HP 2100
C      THREE WORD DOUBLE PRECISION OPTION WITH FTN4
C
C      DATA SMALL(1), SMALL(2), SMALL(3) / 40000B,       0,       1 /
C      DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B /
C      DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B,       0,    265B /
C      DATA DIVER(1), DIVER(2), DIVER(3) / 40000B,       0,    276B /
C      DATA LOG10(1), LOG10(2), LOG10(3) / 46420B,  46502B,  77777B /
C
C
C      MACHINE CONSTANTS FOR THE HP 2100
C      FOUR WORD DOUBLE PRECISION OPTION WITH FTN4
C
C      DATA SMALL(1), SMALL(2) /  40000B,       0 /
C      DATA SMALL(3), SMALL(4) /       0,       1 /
C      DATA LARGE(1), LARGE(2) /  77777B, 177777B /
C      DATA LARGE(3), LARGE(4) / 177777B, 177776B /
C      DATA RIGHT(1), RIGHT(2) /  40000B,       0 /
C      DATA RIGHT(3), RIGHT(4) /       0,    225B /
C      DATA DIVER(1), DIVER(2) /  40000B,       0 /
C      DATA DIVER(3), DIVER(4) /       0,    227B /
C      DATA LOG10(1), LOG10(2) /  46420B,  46502B /
C      DATA LOG10(3), LOG10(4) /  76747B, 176377B /
C
C
C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND
C     THE PERKIN ELMER (INTERDATA) 7/32.
C
C     DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 /
C     DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF /
C     DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 /
C     DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 /
C     DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
C
C     DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 /
C     DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 /
C     DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 /
C     DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 /
C     DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
C
C     DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 /
C     DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 /
C     DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 /
C     DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 /
C     DATA LOG10(1),LOG10(2) / "177464202324, "476747767461 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1),SMALL(2) /    8388608,           0 /
C     DATA LARGE(1),LARGE(2) / 2147483647,          -1 /
C     DATA RIGHT(1),RIGHT(2) /  612368384,           0 /
C     DATA DIVER(1),DIVER(2) /  620756992,           0 /
C     DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /
C
C     DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 /
C     DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 /
C     DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 /
C     DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 /
C     DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1),SMALL(2) /    128,      0 /
C     DATA SMALL(3),SMALL(4) /      0,      0 /
C
C     DATA LARGE(1),LARGE(2) /  32767,     -1 /
C     DATA LARGE(3),LARGE(4) /     -1,     -1 /
C
C     DATA RIGHT(1),RIGHT(2) /   9344,      0 /
C     DATA RIGHT(3),RIGHT(4) /      0,      0 /
C
C     DATA DIVER(1),DIVER(2) /   9472,      0 /
C     DATA DIVER(3),DIVER(4) /      0,      0 /
C
C     DATA LOG10(1),LOG10(2) /  16282,   8348 /
C     DATA LOG10(3),LOG10(4) / -31493, -12296 /
C
C     DATA SMALL(1),SMALL(2) / O000200, O000000 /
C     DATA SMALL(3),SMALL(4) / O000000, O000000 /
C
C     DATA LARGE(1),LARGE(2) / O077777, O177777 /
C     DATA LARGE(3),LARGE(4) / O177777, O177777 /
C
C     DATA RIGHT(1),RIGHT(2) / O022200, O000000 /
C     DATA RIGHT(3),RIGHT(4) / O000000, O000000 /
C
C     DATA DIVER(1),DIVER(2) / O022400, O000000 /
C     DATA DIVER(3),DIVER(4) / O000000, O000000 /
C
C     DATA LOG10(1),LOG10(2) / O037632, O020232 /
C     DATA LOG10(3),LOG10(4) / O102373, O147770 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER
C
C     DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 /
C     DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 /
C     DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 /
C     DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 /
C     DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FOR COMPILER
C
C     DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 /
C     DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 /
C     DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 /
C     DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 /
C     DATA LOG10(1), LOG10(2) / O177746420232, O411757177572/
C
C
C     MACHINE CONSTANTS FOR VAX 11/780
C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
C    ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***
C    *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C     DATA SMALL(1), SMALL(2) /        128,           0 /
C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
C     DATA RIGHT(1), RIGHT(2) /       9344,           0 /
C     DATA DIVER(1), DIVER(2) /       9472,           0 /
C     DATA LOG10(1), LOG10(2) /  546979738,  -805665541 /
C
C     DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 /
C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 /
C     DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 /
C     DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFFA84FB /
C
C***FIRST EXECUTABLE STATEMENT  D1MACH
C
      D1MACH = DMACH(I)
      RETURN
C
      END