C     ALGORITHM 602, COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 3,
C     SEP., 1983, P. 355-357.
      SUBROUTINE HURRY(A, L, M, SUM, F, N, ERR, SIGMAS, DA)             HUR   10
C
C***********************************************************************
C*                                                                     *
C* COMPUTES THE ACCELERATED SUM OF A SERIES OR LIMIT OF A SEQUENCE.    *
C*                                                                     *
C* ARGUMENTS:                                                          *
C*    A      = ARRAY OF ELEMENTS OF SERIES OR SEQUENCE                 *
C*    L      = LEAST NUMBER OF TERMS TO BE USED                        *
C*    M      = MOST NUMBER OF TERMS TO BE USED                         *
C*    SUM    = .TRUE. FOR A SUM, .FALSE. FOR A SEQUENCE                *
C*    F      = FINAL VALUE RETURNED                                    *
C*    N      = NUMBER OF TERMS GIVING 'BEST' RESULT                    *
C*    ERR    = ESTIMATED UNCERTAINTY OF RESULT                         *
C*    SIGMAS = .TRUE. IF UNCERTAINTIES PROVIDED IN DA                  *
C*    DA     = ARRAY OF UNCERTAINTIES OF ELEMENTS                      *
C*                                                                     *
C* INPUTS FROM DRIVER PROGRAM: A, L, M, SUM, SIGMAS, DA                *
C*                                                                     *
C* OUTPUTS TO DRIVER PROGRAM: F, N, ERR                                *
C*                                                                     *
C***********************************************************************
C
C
      LOGICAL SUM, SIGMAS, BETTER, BEGUN, BEFORE, LARGE
      INTEGER I, J, L, M, N
      DOUBLE PRECISION A(M), DA(M), RESULT(50), TRUNC(50), DRDA(50),
     * NOISE(50), W1(50), W2(50), WW1(50,50), WW2(50,50), F, ERR, S,
     * HUGE, SMALL, FACTOR, TEST
C
      DATA HUGE, SMALL, FACTOR /1.0D75,0.01D0,3.0D0/
C
C          HUGE IS USED TO START THE TEST FOR LEAST TRUNCATION
C          ERROR SO FAR.  ANY ALLOWABLE NUMBER LARGER THAN ANY
C          CONCEIVABLE TRUNCATION ERROR WILL DO.  SEE SECTION 2
C          OF ACCOMPANYING ARTICLE FOR EXPLANATION OF SMALL AND
C          FACTOR.
C
      DO 100 I=L,M
C        (COMPUTE RESULT OF ACCELERATING I TERMS)
        IF (I.EQ.L) CALL WHIZ(A, I, W1, W2, SUM, RESULT(I), S, SIGMAS,
     *   DRDA, WW1, WW2, 50)
        IF (I.NE.L) CALL WHIZ1(A, I, W1, W2, SUM, RESULT(I), S, SIGMAS,
     *   DRDA, WW1, WW2, 50)
        NOISE(I) = 0.0D0
C
        IF (.NOT.SIGMAS) GO TO 20
        DO 10 J=1,I
C
C                        ***  WARNING  ***
C              THE NEXT EXECUTABLE STATEMENT CAN AND DOES CAUSE
C              UNDERFLOWS -- THE APPROPRIATE FIX IS MACHINE
C              DEPENDENT.  IF THE MACHINE SETS UNDERFLOWED
C              QUANTITY TO ZERO, NO HARM RESULTS.
C
          NOISE(I) = NOISE(I) + (DA(J)*DRDA(J))**2
   10   CONTINUE
        IF (NOISE(I).GT.0.0) NOISE(I) = DSQRT(NOISE(I))
C
   20   CONTINUE
C              (CHECK TRUNCATION ERROR AND CONVERGENCE)
        IF (I.LE.L) GO TO 30
        TRUNC(I) = DABS(RESULT(I)-RESULT(I-1))
        BETTER = (TRUNC(I).LT.TRUNC(I-1)) .OR. (TRUNC(I).LT.SMALL*
     *   DABS(RESULT(I)))
        BEGUN = BEGUN .OR. (BETTER .AND. BEFORE)
        GO TO 40
   30   CONTINUE
        TRUNC(I) = 0.0D0
        BETTER = .FALSE.
        BEGUN = .FALSE.
        TEST = HUGE
   40   CONTINUE
        BEFORE = BETTER
C
        IF (BEGUN) GO TO 50
        N = I
        GO TO 90
   50   CONTINUE
C              (TEST NUMBER OF TERMS GIVING BEST RESULT SO FAR)
        IF (TRUNC(I).GE.TEST) GO TO 60
        N = I
        TEST = TRUNC(I)
        GO TO 70
   60   CONTINUE
        IF (N.NE.I-1) TEST = (TEST+TRUNC(I))/2.0D0
   70   CONTINUE
        IF (I.EQ.N) GO TO 80
C                 (IS NOISE FROM FROM TERMS LARGE YET?)
        LARGE = TEST.LT.FACTOR*NOISE(I)
        IF (LARGE) GO TO 110
   80   CONTINUE
   90   CONTINUE
  100 CONTINUE
C
  110 CONTINUE
      F = RESULT(N)
      ERR = DMAX1(TRUNC(N),NOISE(N))
C
      RETURN
      END
      SUBROUTINE WHIZ(A, N, QNUM, QDEN, SUM, VALUE, S, DERIVS, DVDA,    WHI   10
     * DQNUM, DQDEN, M)
C
C***********************************************************************
C*                                                                     *
C* THE U ALGORITHM FOR ACCELERATING A SERIES OR A LIMIT OF A SEQUENCE. *
C*                                                                     *
C* ARGUMENTS:                                                          *
C*    A      = ARRAY OF ELEMENTS OF SERIES OR SEQUENCE                 *
C*    N      = NUMBER OF ELEMENTS IN A                                 *
C*    QNUM   = BACKWARD DIAGONAL OF NUMERATOR ARRAY, AT LEAST N LONG   *
C*    QDEN   = BACKWARD DIAGONAL OF DENOMINATOR ARRAY, AT LEAST N LONG *
C*    SUM    = .TRUE. FOR A SUM, .FALSE. FOR A SEQUENCE                *
C*    VALUE  = ACCELERATED VALUE OF A SUM OR LIMIT                     *
C*    S      = PARTIAL SUM OF SERIES                                   *
C*    DERIVS = .TRUE. IF DERIVATIVES ARE TO BE CALCULATED              *
C*    DVDA   = ARRAY OF CALCULATED DERIVATIVES, D(VALUE)/DA            *
C*    DQNUM  = WORKING STORAGE ARRAY, N*N LONG                         *
C*             (USED ONLY IF DERIVS = .TRUE.)                          *
C*    DQDEN  = WORKING STORAGE ARRAY, N*N LONG                         *
C*             (USED ONLY IF DERIVS = .TRUE.)                          *
C*    M      = FIRST DIMENSION OF DQNUM AND DQDEN ARRAYS               *
C*             (DQNUM AND DQDEN ARE DIMENSIONED (M,N), AND M MUST BE   *
C*              AT LEAST AS LARGE AS THE LARGEST N TO BE USED)         *
C*                                                                     *
C* INPUTS FROM DRIVER, PASSED BY HURRY:  A, N, SUM, DERIVS             *
C*                                                                     *
C* INPUT FROM HURRY:  M                                                *
C*                                                                     *
C* OUTPUTS PASSED TO WHIZ1 BY HURRY:  QNUM, QDEN, VALUE, S, DVDA,      *
C*                                    DQNUM, DQDEN                     *
C*                                                                     *
C***********************************************************************
C
      LOGICAL SUM, DERIVS
      INTEGER N, M, NEXT, I, L
      DOUBLE PRECISION A(N), QNUM(N), QDEN(N), VALUE, S, DVDA(N),
     * DQNUM(M,N), DQDEN(M,N), TERM, FNEXT, FL, RATIO, FJ, FACTOR, C,
     * DTERM, DS
C
      S = 0.0D0
C
      DO 130 NEXT=1,N
C           (GET NEXT DIAGONAL)
        IF (SUM) GO TO 10
        TERM = A(NEXT) - S
        S = A(NEXT)
        GO TO 20
   10   CONTINUE
        TERM = A(NEXT)
        S = A(NEXT) + S
   20   CONTINUE
        L = NEXT - 1
        FNEXT = FLOAT(NEXT)
        QDEN(NEXT) = 1.0D0/(TERM*FNEXT**2)
        QNUM(NEXT) = S*QDEN(NEXT)
        IF (.NOT.DERIVS) GO TO 80
        DO 70 I=1,NEXT
          IF (I.NE.NEXT) GO TO 30
          DTERM = 1.0D0
          DS = 1.0D0
          GO TO 60
   30     CONTINUE
          IF (SUM) GO TO 40
          IF (I.EQ.L) DTERM = -1.0D0
          IF (I.NE.L) DTERM = 0.0D0
          DS = 0.0D0
          GO TO 50
   40     CONTINUE
          DTERM = 0.0D0
          DS = 1.0D0
   50     CONTINUE
   60     CONTINUE
          DQDEN(I,NEXT) = -QDEN(I)*DTERM/TERM
          DQNUM(I,NEXT) = DQDEN(I,NEXT)*S + QDEN(NEXT)*DS
   70   CONTINUE
   80   CONTINUE
        IF (NEXT.LE.1) GO TO 120
        FACTOR = 1.0D0
        FL = FLOAT(L)
        RATIO = FL/FNEXT
        LPLUS1 = L + 1
        DO 110 K=1,L
          J = LPLUS1 - K
          FJ = FLOAT(J)
          C = FACTOR*FJ/FNEXT
          FACTOR = FACTOR*RATIO
          QDEN(J) = QDEN(J+1) - C*QDEN(J)
          QNUM(J) = QNUM(J+1) - C*QNUM(J)
          IF (.NOT.DERIVS) GO TO 100
          DO 90 I=1,L
            DQDEN(I,J) = DQDEN(I,J+1) - C*DQDEN(I,J)
            DQNUM(I,J) = DQNUM(I,J+1) - C*DQNUM(I,J)
   90     CONTINUE
          DQDEN(NEXT,J) = DQDEN(NEXT,J+1)
          DQNUM(NEXT,J) = DQNUM(NEXT,J+1)
  100     CONTINUE
  110   CONTINUE
  120   CONTINUE
  130 CONTINUE
C
      VALUE = QNUM(1)/QDEN(1)
      IF (.NOT.DERIVS) GO TO 150
      DO 140 I=1,N
        DVDA(I) = (DQNUM(I,1)-VALUE*DQDEN(I,1))/QDEN(1)
  140 CONTINUE
  150 CONTINUE
C
      RETURN
      END
      SUBROUTINE WHIZ1(A, N, QNUM, QDEN, SUM, VALUE, S, DERIVS, DVDA,   WHI   10
     * DQNUM, DQDEN, M)
C
C***********************************************************************
C*                                                                     *
C* THE U ALGORITHM FOR ACCELERATING A SERIES OR A LIMIT OF A SEQUENCE. *
C* THIS SUBROUTINE IS USED TO GET THE N-TERMS RESULT FROM THE          *
C* RESULT OF N-1 TERMS.  WHIZ1 DIFFERS FROM WHIZ IN THAT:              *
C*  (1) THE ALGORITHM IS RUN FOR NEXT=N RATHER THAN FOR NEXT=1 TO N    *
C*  (2) S IS NOT ZEROED AT THE START OF THE SUBROUTINE                 *
C*                                                                     *
C* ARGUMENTS:                                                          *
C*    A      = ARRAY OF ELEMENTS OF SERIES OR SEQUENCE                 *
C*    N      = NUMBER OF ELEMENTS IN A                                 *
C*    QNUM   = BACKWARD DIAGONAL OF NUMERATOR ARRAY, AT LEAST N LONG   *
C*    QDEN   = BACKWARD DIAGONAL OF DENOMINATOR ARRAY, AT LEAST N LONG *
C*    SUM    = .TRUE. FOR A SUM, .FALSE. FOR A LIMIT                   *
C*    VALUE  = ACCELERATED VALUE OF A SUM OR LIMIT                     *
C*    S      = SIMPLE SUM OF SERIES                                    *
C*    DERIVS = .TRUE. IF DERIVATIVES ARE TO BE CALCULATED              *
C*    DVDA   = ARRAY OF CALCULATED DERIVATIVES, D(VALUE)/DA            *
C*    DQNUM  = WORKING STORAGE ARRAY, N*N LONG                         *
C*             (USED ONLY IF DERIVS = .TRUE.)                          *
C*    DQDEN  = WORKING STORAGE ARRAY, N*N LONG                         *
C*             (USED ONLY IF DERIVS = .TRUE.)                          *
C*    M      = FIRST DIMENSION OF DQNUM AND DQDEN ARRAYS               *
C*             (DQNUM AND DQDEN ARE DIMENSIONED (M,N), AND M MUST BE   *
C*              AT LEAST AS LARGE AS THE LARGEST N TO BE USED)         *
C*                                                                     *
C* INPUTS FROM DRIVER, PASSED BY HURRY:  A, N, SUM, DERIVS             *
C*                                                                     *
C* INPUT FROM HURRY:  M                                                *
C*                                                                     *
C* INPUTS FROM WHIZ, PASSED BY HURRY:  QNUM, QDEN, VALUE, S, DVDA,     *
C*                                     DQNUM, DQDEN                    *
C*                                                                     *
C* OUTPUTS TO HURRY:  VALUE, DVDA                                      *
C*                                                                     *
C***********************************************************************
C
      LOGICAL SUM, DERIVS
      INTEGER N, M, NEXT, I, L
      DOUBLE PRECISION A(N), QNUM(N), QDEN(N), VALUE, S, DVDA(N),
     * DQNUM(M,N), DQDEN(M,N), TERM, FNEXT, FL, RATIO, FJ, FACTOR, C,
     * DTERM, DS
C
      NEXT = N
      IF (SUM) GO TO 10
      TERM = A(NEXT) - S
      S = A(NEXT)
      GO TO 20
   10 CONTINUE
      TERM = A(NEXT)
      S = A(NEXT) + S
   20 CONTINUE
      L = NEXT - 1
      FNEXT = FLOAT(NEXT)
      QDEN(NEXT) = 1.0/(TERM*FNEXT**2)
      QNUM(NEXT) = S*QDEN(NEXT)
      IF (.NOT.DERIVS) GO TO 80
      DO 70 I=1,NEXT
        IF (I.NE.NEXT) GO TO 30
        DTERM = 1.0D0
        DS = 1.0D0
        GO TO 60
   30   CONTINUE
        IF (SUM) GO TO 40
        IF (I.EQ.L) DTERM = -1.0D0
        IF (I.NE.L) DTERM = 0.0D0
        DS = 0.0D0
        GO TO 50
   40   CONTINUE
        DTERM = 0.0D0
        DS = 1.0D0
   50   CONTINUE
   60   CONTINUE
        DQDEN(I,NEXT) = -QDEN(I)*DTERM/TERM
        DQNUM(I,NEXT) = DQDEN(I,NEXT)*S + QDEN(NEXT)*DS
   70 CONTINUE
   80 CONTINUE
      IF (NEXT.LE.1) GO TO 120
      FACTOR = 1.0D0
      FL = FLOAT(L)
      RATIO = FL/FNEXT
      LPLUS1 = L + 1
      DO 110 K=1,L
        J = LPLUS1 - K
        FJ = FLOAT(J)
        C = FACTOR*FJ/FNEXT
        FACTOR = FACTOR*RATIO
        QDEN(J) = QDEN(J+1) - C*QDEN(J)
        QNUM(J) = QNUM(J+1) - C*QNUM(J)
        IF (.NOT.DERIVS) GO TO 100
        DO 90 I=1,L
          DQDEN(I,J) = DQDEN(I,J+1) - C*DQDEN(I,J)
          DQNUM(I,J) = DQNUM(I,J+1) - C*DQNUM(I,J)
   90   CONTINUE
        DQDEN(NEXT,J) = DQDEN(NEXT,J+1)
        DQNUM(NEXT,J) = DQNUM(NEXT,J+1)
  100   CONTINUE
  110 CONTINUE
  120 CONTINUE
C
      VALUE = QNUM(1)/QDEN(1)
      IF (.NOT.DERIVS) GO TO 140
      DO 130 I=1,N
        DVDA(I) = (DQNUM(I,1)-VALUE*DQDEN(I,1))/QDEN(1)
  130 CONTINUE
  140 CONTINUE
C
      RETURN
      END
C                   XACCEL--DRIVER PROGRAM FOR HURRY                    MAN   10
C                                                                       MAN   20
C***********************************************************************MAN   30
C*                                                                     *MAN   40
C* DEMONSTRATES THE USE OF SUBROUTINE HURRY FOR ACCELERATING SUMS AND  *MAN   50
C* LIMITS OF SEQUENCES.                                                *MAN   60
C*                                                                     *MAN   70
C* INPUT:                                                              *MAN   80
C* THIS 'DRIVER' PROGRAM ALLOWS THE USER TO TEST THE PERFORMANCE OF    *MAN   90
C* HURRY FOR ONE OR MORE BUILT-IN SEQUENCES.  FOR EACH TEST DESIRED,   *MAN  100
C* THE USER PROVIDES A DATA CARD IN THE FOLLOWING FORMAT:              *MAN  110
C*                                                                     *MAN  120
C*    ITEM             COLUMNS        FORMAT         RANGE             *MAN  130
C*                                                                     *MAN  140
C*    EXAMPLE CODE      1 - 10          I10          1 - 10            *MAN  150
C*    MINIMUM NUMBER   11 - 20          I10          2 - 50            *MAN  160
C*       OF TERMS                                                      *MAN  170
C*    MAXIMUM NUMBER   21 - 30          I10          2 - 50            *MAN  180
C*       OF TERMS                                                      *MAN  190
C*    X                31 - 40          G10.0          -               *MAN  200
C*    NOISE            41 - 50          G10.0          -               *MAN  210
C*                                                                     *MAN  220
C* ALL ENTRIES MUST BE RIGHT-JUSTIFIED, OF COURSE.  THE 'EXAMPLE       *MAN  230
C* CODE' IS GIVEN IN THE LIST OF SEQUENCES BELOW.  THE MINIMUM AND     *MAN  240
C* MAXIMUM NUMBER OF TERMS TO BE USED IN HURRY MUST LIE IN THE RANGE   *MAN  250
C* 2 - 50 (INCLUSIVE) AND THE MAXIMUM MUST BE GREATER THAN OR EQUAL    *MAN  260
C* TO THE MINIMUM.  X IS THE POINT AT WHICH THE SERIES OR SEQUENCE IS  *MAN  270
C* TO BE EVALUATED AND SHOULD, THEREFORE, LIE IN ITS DOMAIN.           *MAN  280
C* (EXAMPLES 1, 2, 5, AND 9 MAKE USE OF THE X INPUT; THE OTHERS DO     *MAN  290
C* NOT.)  DIGITAL 'NOISE' IS EXPLAINED IN THE ACCOMPANYING ARTICLE.    *MAN  300
C*                                                                     *MAN  310
C* TEST SEQUENCES AND SERIES PROVIDED:                                 *MAN  320
C*                                                                     *MAN  330
C*    EXAMPLE CODE    SEQUENCE/SERIES                                  *MAN  340
C*                                                                     *MAN  350
C*          1         MACLAURIN SERIES FOR EXP(X)                      *MAN  360
C*          2         CONTINUED FRACTION SEQUENCE FOR EXP(X)           *MAN  370
C*          3         ALTERNATING, LINEARLY CONVERGENT SERIES          *MAN  380
C*          4         ZETA(2) LOGARITHMIC CONVERGENCE                  *MAN  390
C*          5         SERIES FOR -LOG(1-X) (CONVERGENT -1 <= X < 1)    *MAN  400
C*          6         ALTERNATING ASYMPTOTIC SERIES 1                  *MAN  410
C*          7         ALTERNATING ASYMPTOTIC SERIES 2                  *MAN  420
C*          8         ALTERNATING ASYMPTOTIC SERIES 3                  *MAN  430
C*          9         MONOTONE ASYMPTOTIC SERIES, USE X=25, 50, OR 100 *MAN  440
C*         10         SEQUENCE FOR EULER'S CONSTANT                    *MAN  450
C*                                                                     *MAN  460
C* FOR MORE INFORMATION ON THESE SEQUENCES AND SERIES, SEE THE BODY    *MAN  470
C* OF THE PROGRAM.                                                     *MAN  480
C*                                                                     *MAN  490
C* FOR MORE INFORMATION ON THE PROGRAM, CONTACT:                       *MAN  500
C*                             DAVID A SMITH                           *MAN  510
C*                             MATHEMATICS DEPARTMENT                  *MAN  520
C*                             DUKE UNIVERSITY                         *MAN  530
C*                             DURHAM, NC  27706                       *MAN  540
C*                                                                     *MAN  550
C***********************************************************************MAN  560
C* THIS PROGRAM CONTAINS SUBROUTINES FROM THE IMSL LIBRARY, A PROPRIE- *MAN  570
C* TARY PACKAGE FROM INTERNATIONAL MATHEMATICAL & STATISTICAL LIBRAR-  *MAN  580
C* IES, INC., HOUSTON, TEXAS.  THESE ROUTINES MAY NOT BE REDISTRIBUTED *MAN  590
C* OR REMOVED FROM THIS SOFTWARE FOR USE IN OTHER SOFTWARE DEVELOPMENT.*MAN  600
C* IMSL ROUTINES INCLUDED ARE GGNQF, MDNRIS, MERFI, UERTST, UGETIO.    *MAN  610
C*                                                                     *MAN  620
C***********************************************************************MAN  630
C                                                                       MAN  640
C                                                                       MAN  650
      INTEGER EXAMPL, MIN, MAX, I                                       MAN  660
      DOUBLE PRECISION NOISE, X, TERM(50), DA(50), ANSWER, AI, BI, AYE, MAN  670
     * TOP(3), BOT(3), FACTOR, PTSUM, TINY, ENORM, PI, ARG, RESULT,     MAN  680
     * ESTERR, ACTERR, DIGITS, ANS, FUZZ, DSEED                         MAN  690
      LOGICAL OK, SIGMAS, SUM                                           MAN  700
C                                                                       MAN  710
      DIGITS(ARG,ENORM) = -DLOG10(DMAX1(DABS(ARG/ENORM),TINY))          MAN  720
C                                                                       MAN  730
      TINY = 1.0D-16                                                    MAN  740
      DSEED = 9973.D0                                                   MAN  750
      PI = 3.141592653589793D0                                          MAN  760
      FUZZ = 1.0D-5                                                     MAN  770
C                                                                       MAN  780
C                                                                       MAN  790
C        (THE ENTIRE PROGRAM IS REPEATED FOR EACH INPUT CARD)           MAN  800
C                                                                       MAN  810
   10 READ(5,99999,END=270) EXAMPL, MIN, MAX, X, NOISE                  MAN  820
      WRITE (6,99998) EXAMPL, X, NOISE, MIN, MAX                        MAN  830
C                                                                       MAN  840
      SIGMAS = .FALSE.                                                  MAN  850
C                                                                       MAN  860
C                                                                       MAN  870
C        CHECK INPUT VALUES                                             MAN  880
C                                                                       MAN  890
      OK = (MIN.GE.2) .AND. (MIN.LE.MAX) .AND. (MAX.LE.50) .AND.        MAN  900
     * (EXAMPL.GE.1) .AND. (EXAMPL.LE.10)                               MAN  910
      IF (OK) GO TO (20, 40, 80, 100, 120, 140, 160, 180, 200, 220),    MAN  920
     * EXAMPL                                                           MAN  930
      WRITE (6,99997)                                                   MAN  940
      GO TO 10                                                          MAN  950
C                                                                       MAN  960
C                                                                       MAN  970
C**********EXAMPLE 1:  MACLAURIN SERIES FOR EXP(X)                      MAN  980
C                                                                       MAN  990
   20 TERM(1) = 1.0D0                                                   MAN 1000
      DO 30 I=2,MAX                                                     MAN 1010
        AYE = FLOAT(I)                                                  MAN 1020
        TERM(I) = TERM(I-1)*X/(AYE-1.0D0)                               MAN 1030
   30 CONTINUE                                                          MAN 1040
      SUM = .TRUE.                                                      MAN 1050
      ANSWER = DEXP(X)                                                  MAN 1060
      GO TO 240                                                         MAN 1070
C                                                                       MAN 1080
C                                                                       MAN 1090
C**********EXAMPLE 2:  CONTINUED-FRACTION SEQUENCE FOR EXP(X)           MAN 1100
C                                                                       MAN 1110
   40 TERM(1) = 1.0D0                                                   MAN 1120
      TOP(1) = 1.0D0                                                    MAN 1130
      TOP(2) = 0.0D0                                                    MAN 1140
      TOP(3) = 1.0D0                                                    MAN 1150
      BOT(1) = 0.0D0                                                    MAN 1160
      BOT(2) = 1.0D0                                                    MAN 1170
      BOT(3) = 1.0D0                                                    MAN 1180
      DO 70 I=2,MAX                                                     MAN 1190
        IF (MOD(I,2).EQ.0) GO TO 50                                     MAN 1200
        AI = X                                                          MAN 1210
        BI = 2.0D0                                                      MAN 1220
        GO TO 60                                                        MAN 1230
   50   CONTINUE                                                        MAN 1240
        AI = -X                                                         MAN 1250
        BI = FLOAT(I-1)                                                 MAN 1260
   60   CONTINUE                                                        MAN 1270
        TOP(1) = TOP(2)                                                 MAN 1280
        TOP(2) = TOP(3)                                                 MAN 1290
        BOT(1) = BOT(2)                                                 MAN 1300
        BOT(2) = BOT(3)                                                 MAN 1310
        TOP(3) = BI*TOP(2) + AI*TOP(1)                                  MAN 1320
        BOT(3) = BI*BOT(2) + AI*BOT(1)                                  MAN 1330
        TERM(I) = TOP(3)/BOT(3)                                         MAN 1340
   70 CONTINUE                                                          MAN 1350
      SUM = .FALSE.                                                     MAN 1360
      ANSWER = DEXP(X)                                                  MAN 1370
      GO TO 240                                                         MAN 1380
C                                                                       MAN 1390
C                                                                       MAN 1400
C**********EXAMPLE 3:  ALTERNATING LINEARLY CONVERGENT SERIES           MAN 1410
C                                                                       MAN 1420
   80 FACTOR = -1.0D0                                                   MAN 1430
      DO 90 I=1,MAX                                                     MAN 1440
        AYE = FLOAT(I)                                                  MAN 1450
        FACTOR = -FACTOR                                                MAN 1460
        TERM(I) = FACTOR/DSQRT(AYE)                                     MAN 1470
   90 CONTINUE                                                          MAN 1480
      SUM = .TRUE.                                                      MAN 1490
      ANSWER = 0.60489864342162D0                                       MAN 1500
      GO TO 240                                                         MAN 1510
C                                                                       MAN 1520
C                                                                       MAN 1530
C**********EXAMPLE 4:  ZETA(2) LOGARITHMIC CONVERGENCE                  MAN 1540
C                                                                       MAN 1550
  100 DO 110 I=1,MAX                                                    MAN 1560
        AYE = FLOAT(I)                                                  MAN 1570
        TERM(I) = 1.0D0/(AYE*AYE)                                       MAN 1580
  110 CONTINUE                                                          MAN 1590
      SUM = .TRUE.                                                      MAN 1600
      ANSWER = 1.644934066848226D0                                      MAN 1610
      GO TO 240                                                         MAN 1620
C                                                                       MAN 1630
C                                                                       MAN 1640
C**********EXAMPLE 5:  SERIES FOR -LOG(1-X)                             MAN 1650
C                      (CONVERGES ON -1<=X<1)                           MAN 1660
  120 FACTOR = 1.D0                                                     MAN 1670
      DO 130 I=1,MAX                                                    MAN 1680
        AYE = FLOAT(I)                                                  MAN 1690
        FACTOR = FACTOR*X                                               MAN 1700
        TERM(I) = FACTOR/AYE                                            MAN 1710
  130 CONTINUE                                                          MAN 1720
      SUM = .TRUE.                                                      MAN 1730
      ANSWER = -DLOG(1.0D0-X)                                           MAN 1740
      GO TO 240                                                         MAN 1750
C                                                                       MAN 1760
C                                                                       MAN 1770
C**********EXAMPLE 6:  ALTERNATING ASYMPTOTIC SERIES 1                  MAN 1780
  140 TERM(1) = 3.0D0/(PI*PI)                                           MAN 1790
      DO 150 I=2,MAX                                                    MAN 1800
        AYE = FLOAT(I)                                                  MAN 1810
        TERM(I) = -TERM(I-1)*(4.0D0*AYE-1.0D0)/(PI*PI)                  MAN 1820
  150 CONTINUE                                                          MAN 1830
      SUM = .TRUE.                                                      MAN 1840
      ANSWER = 0.19259404877D0                                          MAN 1850
      GO TO 240                                                         MAN 1860
C                                                                       MAN 1870
C                                                                       MAN 1880
C**********EXAMPLE 7:  ALTERNATING ASYMPTOTIC SERIES 2                  MAN 1890
C                                                                       MAN 1900
  160 TERM(1) = 1.0D0                                                   MAN 1910
      DO 170 I=2,MAX                                                    MAN 1920
        AYE = FLOAT(I)                                                  MAN 1930
        TERM(I) = -TERM(I-1)*(2.0D0*AYE-3.0D0)/2.0D0                    MAN 1940
  170 CONTINUE                                                          MAN 1950
      SUM = .TRUE.                                                      MAN 1960
      ANSWER = 0.7578721564D0                                           MAN 1970
      GO TO 240                                                         MAN 1980
C                                                                       MAN 1990
C                                                                       MAN 2000
C**********EXAMPLE 8:  ALTERNATING ASYMPTOTIC SERIES 3                  MAN 2010
C                                                                       MAN 2020
  180 FACTOR = -1.0D0                                                   MAN 2030
      DO 190 I=1,MAX                                                    MAN 2040
        AYE = FLOAT(I)                                                  MAN 2050
        FACTOR = -FACTOR*(AYE+1.0D0)                                    MAN 2060
        TERM(I) = FACTOR*AYE                                            MAN 2070
  190 CONTINUE                                                          MAN 2080
      SUM = .TRUE.                                                      MAN 2090
      ANSWER = 0.210957917D0                                            MAN 2100
      GO TO 240                                                         MAN 2110
C                                                                       MAN 2120
C                                                                       MAN 2130
C**********EXAMPLE 9:  MONOTONE AYMPTOTIC SERIES                        MAN 2140
C                      (X=25, 50, OR 100)                               MAN 2150
C                                                                       MAN 2160
  200 FACTOR = 1.D0                                                     MAN 2170
      DO 210 I=1,MAX                                                    MAN 2180
        AYE = FLOAT(I)                                                  MAN 2190
        FACTOR = FACTOR*(AYE+1.0D0)/X                                   MAN 2200
        TERM(I) = FACTOR*AYE*AYE                                        MAN 2210
  210 CONTINUE                                                          MAN 2220
      SUM = .TRUE.                                                      MAN 2230
      IF (DABS(X-25.0D0).LT.FUZZ) ANSWER = 0.140396959326971D0          MAN 2240
      IF (DABS(X-50.0D0).LT.FUZZ) ANSWER = 0.051707744368624644D0       MAN 2250
      IF (DABS(X-100.0D0).LT.FUZZ) ANSWER = 0.226372038599530D-1        MAN 2260
      GO TO 240                                                         MAN 2270
C                                                                       MAN 2280
C                                                                       MAN 2290
C**********EXAMPLE 10:  SEQUENCE FOR EULER'S CONTANT                    MAN 2300
C                                                                       MAN 2310
  220 PTSUM = 0.0D0                                                     MAN 2320
      DO 230 I=1,MAX                                                    MAN 2330
        AYE = FLOAT(I)                                                  MAN 2340
        PTSUM = PTSUM + 1.0D0/AYE                                       MAN 2350
        TERM(I) = PTSUM - DLOG(AYE)                                     MAN 2360
  230 CONTINUE                                                          MAN 2370
      SUM = .FALSE.                                                     MAN 2380
      ANSWER = 0.5772156649015329D0                                     MAN 2390
C                                                                       MAN 2400
C                                                                       MAN 2410
C        SEQUENCES CREATED...NOW SEE IF NOISE NEEDED:                   MAN 2420
C                                                                       MAN 2430
  240 CONTINUE                                                          MAN 2440
      IF (NOISE.EQ.0.0D0) GO TO 260                                     MAN 2450
      SIGMAS = .TRUE.                                                   MAN 2460
      DO 250 I=1,MAX                                                    MAN 2470
        DA(I) = NOISE*TERM(I)                                           MAN 2480
        IF (NOISE.GT.0.0) TERM(I) = TERM(I) + DA(I)*GGNQF(DSEED)        MAN 2490
C                                                                       MAN 2500
C        THE PRECEDING STATEMENT USES THE IMSL ROUTINE GGNQF,           MAN 2510
C        WHICH CALLS IMSL ROUTINES MDNRIS, MERFI, UERTST, UGETIO.       MAN 2520
C        SEE IMSL NOTICE IN PROGRAM HEADING.                            MAN 2530
C                                                                       MAN 2540
  250 CONTINUE                                                          MAN 2550
C                                                                       MAN 2560
C                                                                       MAN 2570
  260 CONTINUE                                                          MAN 2580
C                                                                       MAN 2590
C                                                                       MAN 2600
      CALL HURRY(TERM, MIN, MAX, SUM, RESULT, NUSED, ESTERR, SIGMAS, DA)MAN 2610
C                                                                       MAN 2620
C                                                                       MAN 2630
      ANS = ANSWER                                                      MAN 2640
      IF (ANS.EQ.0.0D0) ANS = 1.0D0                                     MAN 2650
      ACTERR = RESULT - ANSWER                                          MAN 2660
      DEST = DIGITS(ESTERR,ANS)                                         MAN 2670
      DACT = DIGITS(ACTERR,ANS)                                         MAN 2680
C                                                                       MAN 2690
      WRITE (6,99996) NUSED, DEST, RESULT, DACT, ANSWER                 MAN 2700
C                                                                       MAN 2710
      GO TO 10                                                          MAN 2720
C                                                                       MAN 2730
      WRITE (6,99995)                                                   MAN 2740
C                                                                       MAN 2750
      STOP                                                              MAN 2760
99999 FORMAT (3I10, 2G10.0)                                             MAN 2770
99998 FORMAT (////9H EXAMPLE:, I3, 5H; X =, F10.5, 9H, NOISE =,         MAN 2780
     * G10.3//4X, 23H NUMBER OF TERMS: MIN =, I3/21X, 6H MAX =, I3)     MAN 2790
99997 FORMAT (44H ERROR IN INPUT VALUE OF MIN, MAX, OR EXAMPL)          MAN 2800
99996 FORMAT (18X, 9H ACTUAL =, I3//4X, 27H DIGITS ESTIMATED ACCURACY:, MAN 2810
     * F7.2, 5X, 17H COMPUTED RESULT:, G25.16/4X, 18H DIGITS ACTUAL ACC,MAN 2820
     * 6HURACY:, 3X, F7.2, 5X, 16H CORRECT ANSWER:, G26.16)             MAN 2830
99995 FORMAT (////23H ALL EXAMPLES PROCESSED)                           MAN 2840
      END                                                               MAN 2850
C   IMSL ROUTINE NAME   - GGNQF                                         GGN   10
C                                                                       GGN   20
C-----------------------------------------------------------------------GGN   30
C                                                                       GGN   40
C   COMPUTER            - IBM/SINGLE                                    GGN   50
C                                                                       GGN   60
C   LATEST REVISION     - JUNE 1, 1980                                  GGN   70
C                                                                       GGN   80
C   PURPOSE             - NORMAL RANDOM DEVIATE GENERATOR - FUNCTION    GGN   90
C                           FORM OF GGNML                               GGN  100
C                                                                       GGN  110
C   USAGE               - FUNCTION GGNQF (DSEED)                        GGN  120
C                                                                       GGN  130
C   ARGUMENTS    GGNQF  - RESULTANT NORMAL (0,1) DEVIATE.               GGN  140
C                DSEED  - INPUT/OUTPUT DOUBLE PRECISION VARIABLE        GGN  150
C                           ASSIGNED AN INTEGER VALUE IN THE            GGN  160
C                           EXCLUSIVE RANGE (1.D0, 2147483647.D0).      GGN  170
C                           DSEED IS REPLACED BY A NEW VALUE TO BE      GGN  180
C                           USED IN A SUBSEQUENT CALL.                  GGN  190
C                                                                       GGN  200
C   PRECISION/HARDWARE  - SINGLE/ALL                                    GGN  210
C                                                                       GGN  220
C   REQD. IMSL ROUTINES - MDNRIS,MERFI,UERTST,UGETIO                    GGN  230
C                                                                       GGN  240
C   NOTATION            - INFORMATION ON SPECIAL NOTATION AND           GGN  250
C                           CONVENTIONS IS AVAILABLE IN THE MANUAL      GGN  260
C                           INTRODUCTION OR THROUGH IMSL ROUTINE UHELP  GGN  270
C                                                                       GGN  280
C   COPYRIGHT           - 1980 BY IMSL, INC. ALL RIGHTS RESERVED.       GGN  290
C                                                                       GGN  300
C   WARRANTY            - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN GGN  310
C                           APPLIED TO THIS CODE. NO OTHER WARRANTY,    GGN  320
C                           EXPRESSED OR IMPLIED, IS APPLICABLE.        GGN  330
C                                                                       GGN  340
C-----------------------------------------------------------------------GGN  350
C                                                                       GGN  360
      REAL FUNCTION GGNQF(DSEED)                                        GGN  370
C                                  SPECIFICATIONS FOR ARGUMENTS
      DOUBLE PRECISION DSEED
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER IER
      DOUBLE PRECISION D2P31M, D2PN31
C                                  D2P31M = (2**31) - 1
C                                  D2PN31 = (2**31)
      DATA D2P31M /2147483647.D0/, D2PN31 /2147483648.D0/
C                                  GENERATE A RANDOM (0,1) DEVIATE.
C                                  16807 = (7**5)
C                                  FIRST EXECUTABLE STATEMENT
      DSEED = DMOD(16807.D0*DSEED,D2P31M)
      GGNQF = DSEED/D2PN31
C                                  TRANSFORM TO NORMAL DEVIATE
      CALL MDNRIS(GGNQF, GGNQF, IER)
      RETURN
      END
C   IMSL ROUTINE NAME   - MDNRIS                                        MDN   10
C                                                                       MDN   20
C-----------------------------------------------------------------------MDN   30
C                                                                       MDN   40
C   COMPUTER            - IBM/SINGLE                                    MDN   50
C                                                                       MDN   60
C   LATEST REVISION     - JUNE 1, 1981                                  MDN   70
C                                                                       MDN   80
C   PURPOSE             - INVERSE STANDARD NORMAL (GAUSSIAN)            MDN   90
C                           PROBABILITY DISTRIBUTION FUNCTION           MDN  100
C                                                                       MDN  110
C   USAGE               - CALL MDNRIS (P,Y,IER)                         MDN  120
C                                                                       MDN  130
C   ARGUMENTS    P      - INPUT VALUE IN THE EXCLUSIVE RANGE (0.0,1.0)  MDN  140
C                Y      - OUTPUT VALUE OF THE INVERSE NORMAL (0,1)      MDN  150
C                           PROBABILITY DISTRIBUTION FUNCTION           MDN  160
C                IER    - ERROR PARAMETER (OUTPUT)                      MDN  170
C                         TERMINAL ERROR                                MDN  180
C                           IER = 129 INDICATES P LIES OUTSIDE THE LEGALMDN  190
C                             RANGE. PLUS OR MINUS MACHINE INFINITY IS  MDN  200
C                             GIVEN AS THE RESULT (SIGN IS THE SIGN OF  MDN  210
C                             THE FUNCTION VALUE OF THE NEAREST LEGAL   MDN  220
C                             ARGUMENT).                                MDN  230
C                                                                       MDN  240
C   PRECISION/HARDWARE  - SINGLE/ALL                                    MDN  250
C                                                                       MDN  260
C   REQD. IMSL ROUTINES - MERFI,UERTST,UGETIO                           MDN  270
C                                                                       MDN  280
C   NOTATION            - INFORMATION ON SPECIAL NOTATION AND           MDN  290
C                           CONVENTIONS IS AVAILABLE IN THE MANUAL      MDN  300
C                           INTRODUCTION OR THROUGH IMSL ROUTINE UHELP  MDN  310
C                                                                       MDN  320
C   COPYRIGHT           - 1978 BY IMSL, INC. ALL RIGHTS RESERVED.       MDN  330
C                                                                       MDN  340
C   WARRANTY            - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN MDN  350
C                           APPLIED TO THIS CODE. NO OTHER WARRANTY,    MDN  360
C                           EXPRESSED OR IMPLIED, IS APPLICABLE.        MDN  370
C                                                                       MDN  380
C-----------------------------------------------------------------------MDN  390
C                                                                       MDN  400
      SUBROUTINE MDNRIS(P, Y, IER)                                      MDN  410
C                                  SPECIFICATIONS FOR ARGUMENTS
      REAL P, Y
      INTEGER IER
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      REAL EPS, G0, G1, G2, G3, H0, H1, H2, A, W, WI, SN, SD
      REAL SIGMA, SQRT2, X, XINF
      DATA XINF /Z7FFFFFFF/
      DATA SQRT2 /1.414214/
      DATA EPS /Z3C100000/
      DATA G0 /.1851159E-3/, G1 /-.2028152E-2/
      DATA G2 /-.1498384/, G3 /.1078639E-1/
      DATA H0 /.9952975E-1/, H1 /.5211733/
      DATA H2 /-.6888301E-1/
C                                  FIRST EXECUTABLE STATEMENT
      IER = 0
      IF (P.GT.0.0 .AND. P.LT.1.0) GO TO 10
      IER = 129
      SIGMA = SIGN(1.0,P)
      Y = SIGMA*XINF
      GO TO 30
   10 IF (P.LE.EPS) GO TO 20
      X = 1.0 - (P+P)
      CALL MERFI(X, Y, IER)
      Y = -SQRT2*Y
      GO TO 40
C                                  P TOO SMALL, COMPUTE Y DIRECTLY
   20 A = P + P
      W = SQRT(-ALOG(A+(A-A*A)))
C                                  USE A RATIONAL FUNCTION IN 1./W
      WI = 1./W
      SN = ((G3*WI+G2)*WI+G1)*WI
      SD = ((WI+H2)*WI+H1)*WI + H0
      Y = W + W*(G0+SN/SD)
      Y = -Y*SQRT2
      GO TO 40
   30 CONTINUE
      CALL UERTST(IER, 6HMDNRIS)
   40 RETURN
      END
C   IMSL ROUTINE NAME   - MERFI                                         MER   10
C                                                                       MER   20
C-----------------------------------------------------------------------MER   30
C                                                                       MER   40
C   COMPUTER            - IBM/SINGLE                                    MER   50
C                                                                       MER   60
C   LATEST REVISION     - JANUARY 1, 1978                               MER   70
C                                                                       MER   80
C   PURPOSE             - INVERSE ERROR FUNCTION                        MER   90
C                                                                       MER  100
C   USAGE               - CALL MERFI (P,Y,IER)                          MER  110
C                                                                       MER  120
C   ARGUMENTS    P      - INPUT VALUE IN THE EXCLUSIVE RANGE (-1.0,1.0) MER  130
C                Y      - OUTPUT VALUE OF THE INVERSE ERROR FUNCTION    MER  140
C                IER    - ERROR PARAMETER (OUTPUT)                      MER  150
C                         TERMINAL ERROR                                MER  160
C                           IER = 129 INDICATES P LIES OUTSIDE THE LEGALMER  170
C                             RANGE. PLUS OR MINUS MACHINE INFINITY IS  MER  180
C                             GIVEN AS THE RESULT (SIGN IS THE SIGN OF  MER  190
C                             THE FUNCTION VALUE OF THE NEAREST LEGAL   MER  200
C                             ARGUMENT).                                MER  210
C                                                                       MER  220
C   PRECISION/HARDWARE  - SINGLE/ALL                                    MER  230
C                                                                       MER  240
C   REQD. IMSL ROUTINES - UERTST,UGETIO                                 MER  250
C                                                                       MER  260
C   NOTATION            - INFORMATION ON SPECIAL NOTATION AND           MER  270
C                           CONVENTIONS IS AVAILABLE IN THE MANUAL      MER  280
C                           INTRODUCTION OR THROUGH IMSL ROUTINE UHELP  MER  290
C                                                                       MER  300
C   COPYRIGHT           - 1978 BY IMSL, INC. ALL RIGHTS RESERVED.       MER  310
C                                                                       MER  320
C   WARRANTY            - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN MER  330
C                           APPLIED TO THIS CODE. NO OTHER WARRANTY,    MER  340
C                           EXPRESSED OR IMPLIED, IS APPLICABLE.        MER  350
C                                                                       MER  360
C-----------------------------------------------------------------------MER  370
C                                                                       MER  380
      SUBROUTINE MERFI(P, Y, IER)                                       MER  390
C                                  SPECIFICATIONS FOR ARGUMENTS
      REAL P, Y
      INTEGER IER
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      REAL A, B, X, Z, W, WI, SN, SD, F, Z2, RINFM, A1, A2, A3, B0, B1,
     * B2, B3, C0, C1, C2, C3, D0, D1, D2, E0, E1, E2, E3, F0, F1, F2,
     * G0, G1, G2, G3, H0, H1, H2, SIGMA
      DATA A1 /-.5751703/, A2 /-1.896513/, A3 /-.5496261E-1/
      DATA B0 /-.1137730/, B1 /-3.293474/, B2 /-2.374996/
      DATA B3 /-1.187515/
      DATA C0 /-.1146666/, C1 /-.1314774/, C2 /-.2368201/
      DATA C3 /.5073975E-1/
      DATA D0 /-44.27977/, D1 /21.98546/, D2 /-7.586103/
      DATA E0 /-.5668422E-1/, E1 /.3937021/, E2 /-.3166501/
      DATA E3 /.6208963E-1/
      DATA F0 /-6.266786/, F1 /4.666263/, F2 /-2.962883/
      DATA G0 /.1851159E-3/, G1 /-.2028152E-2/
      DATA G2 /-.1498384/, G3 /.1078639E-1/
      DATA H0 /.9952975E-1/, H1 /.5211733/
      DATA H2 /-.6888301E-1/
      DATA RINFM /Z7FFFFFFF/
C                                  FIRST EXECUTABLE STATEMENT
      IER = 0
      X = P
      SIGMA = SIGN(1.0,X)
C                                  TEST FOR INVALID ARGUMENT
      IF (.NOT.(X.GT.-1. .AND. X.LT.1.)) GO TO 50
      Z = ABS(X)
      IF (Z.LE..85) GO TO 30
      A = 1. - Z
      B = Z
C                                  REDUCED ARGUMENT IS IN (.85,1.),
C                                     OBTAIN THE TRANSFORMED VARIABLE
      W = SQRT(-ALOG(A+A*B))
      IF (W.LT.2.5) GO TO 20
      IF (W.LT.4.) GO TO 10
C                                  W GREATER THAN 4., APPROX. F BY A
C                                     RATIONAL FUNCTION IN 1./W
      WI = 1./W
      SN = ((G3*WI+G2)*WI+G1)*WI
      SD = ((WI+H2)*WI+H1)*WI + H0
      F = W + W*(G0+SN/SD)
      GO TO 40
C                                  W BETWEEN 2.5 AND 4., APPROX. F
C                                     BY A RATIONAL FUNCTION IN W
   10 SN = ((E3*W+E2)*W+E1)*W
      SD = ((W+F2)*W+F1)*W + F0
      F = W + W*(E0+SN/SD)
      GO TO 40
C                                  W BETWEEN 1.13222 AND 2.5, APPROX.
C                                     F BY A RATIONAL FUNCTION IN W
   20 SN = ((C3*W+C2)*W+C1)*W
      SD = ((W+D2)*W+D1)*W + D0
      F = W + W*(C0+SN/SD)
      GO TO 40
C                                  Z BETWEEN 0. AND .85, APPROX. F
C                                     BY A RATIONAL FUNCTION IN Z
   30 Z2 = Z*Z
      F = Z + Z*(B0+A1*Z2/(B1+Z2+A2/(B2+Z2+A3/(B3+Z2))))
C                                  FORM THE SOLUTION BY MULT. F BY
C                                     THE PROPER SIGN
   40 Y = SIGMA*F
      IER = 0
      GO TO 60
C                                  ERROR EXIT. SET SOLUTION TO PLUS
C                                     (OR MINUS) INFINITY
   50 IER = 129
      Y = SIGMA*RINFM
      CONTINUE
      CALL UERTST(IER, 6HMERFI )
   60 RETURN
      END
C   IMSL ROUTINE NAME   - UERTST                                        UER   10
C                                                                       UER   20
C-----------------------------------------------------------------------UER   30
C                                                                       UER   40
C   COMPUTER            - IBM/SINGLE                                    UER   50
C                                                                       UER   60
C   LATEST REVISION     - JANUARY 1, 1978                               UER   70
C                                                                       UER   80
C   PURPOSE             - PRINT A MESSAGE REFLECTING AN ERROR CONDITION UER   90
C                                                                       UER  100
C   USAGE               - CALL UERTST (IER,NAME)                        UER  110
C                                                                       UER  120
C   ARGUMENTS    IER    - ERROR PARAMETER. (INPUT)                      UER  130
C                           IER = I+J WHERE                             UER  140
C                             I = 128 IMPLIES TERMINAL ERROR,           UER  150
C                             I =  64 IMPLIES WARNING WITH FIX, AND     UER  160
C                             I =  32 IMPLIES WARNING.                  UER  170
C                             J = ERROR CODE RELEVANT TO CALLING        UER  180
C                                 ROUTINE.                              UER  190
C                NAME   - A SIX CHARACTER LITERAL STRING GIVING THE     UER  200
C                           NAME OF THE CALLING ROUTINE. (INPUT)        UER  210
C                                                                       UER  220
C   PRECISION/HARDWARE  - SINGLE/ALL                                    UER  230
C                                                                       UER  240
C   REQD. IMSL ROUTINES - UGETIO                                        UER  250
C                                                                       UER  260
C   NOTATION            - INFORMATION ON SPECIAL NOTATION AND           UER  270
C                           CONVENTIONS IS AVAILABLE IN THE MANUAL      UER  280
C                           INTRODUCTION OR THROUGH IMSL ROUTINE UHELP  UER  290
C                                                                       UER  300
C   REMARKS      THE ERROR MESSAGE PRODUCED BY UERTST IS WRITTEN        UER  310
C                ONTO THE STANDARD OUTPUT UNIT. THE OUTPUT UNIT         UER  320
C                NUMBER CAN BE DETERMINED BY CALLING UGETIO AS          UER  330
C                FOLLOWS..   CALL UGETIO(1,NIN,NOUT).                   UER  340
C                THE OUTPUT UNIT NUMBER CAN BE CHANGED BY CALLING       UER  350
C                UGETIO AS FOLLOWS..                                    UER  360
C                                NIN = 0                                UER  370
C                                NOUT = NEW OUTPUT UNIT NUMBER          UER  380
C                                CALL UGETIO(3,NIN,NOUT)                UER  390
C                SEE THE UGETIO DOCUMENT FOR MORE DETAILS.              UER  400
C                                                                       UER  410
C   COPYRIGHT           - 1978 BY IMSL, INC. ALL RIGHTS RESERVED.       UER  420
C                                                                       UER  430
C   WARRANTY            - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN UER  440
C                           APPLIED TO THIS CODE. NO OTHER WARRANTY,    UER  450
C                           EXPRESSED OR IMPLIED, IS APPLICABLE.        UER  460
C                                                                       UER  470
C-----------------------------------------------------------------------UER  480
C                                                                       UER  490
      SUBROUTINE UERTST(IER, NAME)                                      UER  500
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER IER
      INTEGER*2          NAME(3)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER*2          NAMSET(3),NAMEQ(3)
      DATA NAMSET /2HUE,2HRS,2HET/
      DATA NAMEQ /2H  ,2H  ,2H  /
C                                  FIRST EXECUTABLE STATEMENT
      DATA LEVEL /4/, IEQDF /0/, IEQ /1H=/
      IF (IER.GT.999) GO TO 50
      IF (IER.LT.-32) GO TO 70
      IF (IER.LE.128) GO TO 10
      IF (LEVEL.LT.1) GO TO 60
C                                  PRINT TERMINAL MESSAGE
      CALL UGETIO(1, NIN, IOUNIT)
      IF (IEQDF.EQ.1) WRITE (IOUNIT,99999) IER, NAMEQ, IEQ, NAME
      IF (IEQDF.EQ.0) WRITE (IOUNIT,99999) IER, NAME
      GO TO 60
   10 IF (IER.LE.64) GO TO 20
      IF (LEVEL.LT.2) GO TO 60
C                                  PRINT WARNING WITH FIX MESSAGE
      CALL UGETIO(1, NIN, IOUNIT)
      IF (IEQDF.EQ.1) WRITE (IOUNIT,99998) IER, NAMEQ, IEQ, NAME
      IF (IEQDF.EQ.0) WRITE (IOUNIT,99998) IER, NAME
      GO TO 60
   20 IF (IER.LE.32) GO TO 30
C                                  PRINT WARNING MESSAGE
      IF (LEVEL.LT.3) GO TO 60
      CALL UGETIO(1, NIN, IOUNIT)
      IF (IEQDF.EQ.1) WRITE (IOUNIT,99997) IER, NAMEQ, IEQ, NAME
      IF (IEQDF.EQ.0) WRITE (IOUNIT,99997) IER, NAME
      GO TO 60
   30 CONTINUE
C                                  CHECK FOR UERSET CALL
      DO 40 I=1,3
        IF (NAME(I).NE.NAMSET(I)) GO TO 50
   40 CONTINUE
      LEVOLD = LEVEL
      LEVEL = IER
      IER = LEVOLD
      IF (LEVEL.LT.0) LEVEL = 4
      IF (LEVEL.GT.4) LEVEL = 4
      GO TO 60
   50 CONTINUE
      IF (LEVEL.LT.4) GO TO 60
C                                  PRINT NON-DEFINED MESSAGE
      CALL UGETIO(1, NIN, IOUNIT)
      IF (IEQDF.EQ.1) WRITE (IOUNIT,99996) IER, NAMEQ, IEQ, NAME
      IF (IEQDF.EQ.0) WRITE (IOUNIT,99996) IER, NAME
   60 IEQDF = 0
      RETURN
99999 FORMAT (19H *** TERMINAL ERROR, 10X, 7H(IER = , I3, 10H) FROM IMS,
     * 10HL ROUTINE , 3A2, A1, 3A2)
99998 FORMAT (36H *** WARNING WITH FIX ERROR  (IER = , I3, 9H) FROM IM,
     * 11HSL ROUTINE , 3A2, A1, 3A2)
99997 FORMAT (18H *** WARNING ERROR, 11X, 7H(IER = , I3, 11H) FROM IMSL,
     * 9H ROUTINE , 3A2, A1, 3A2)
99996 FORMAT (20H *** UNDEFINED ERROR, 9X, 7H(IER = , I5, 10H) FROM IMS,
     * 10HL ROUTINE , 3A2, A1, 3A2)
C                                  SAVE P FOR P = R CASE
C                                    P IS THE PAGE NAME
C                                    R IS THE ROUTINE NAME
   70 IEQDF = 1
      DO 80 I=1,3
        NAMEQ(I) = NAME(I)
   80 CONTINUE
      RETURN
      END
C   IMSL ROUTINE NAME   - UGETIO                                        UGE   10
C                                                                       UGE   20
C-----------------------------------------------------------------------UGE   30
C                                                                       UGE   40
C   COMPUTER            - IBM/SINGLE                                    UGE   50
C                                                                       UGE   60
C   LATEST REVISION     - JUNE 1, 1981                                  UGE   70
C                                                                       UGE   80
C   PURPOSE             - TO RETRIEVE CURRENT VALUES AND TO SET NEW     UGE   90
C                           VALUES FOR INPUT AND OUTPUT UNIT            UGE  100
C                           IDENTIFIERS.                                UGE  110
C                                                                       UGE  120
C   USAGE               - CALL UGETIO(IOPT,NIN,NOUT)                    UGE  130
C                                                                       UGE  140
C   ARGUMENTS    IOPT   - OPTION PARAMETER. (INPUT)                     UGE  150
C                           IF IOPT=1, THE CURRENT INPUT AND OUTPUT     UGE  160
C                           UNIT IDENTIFIER VALUES ARE RETURNED IN NIN  UGE  170
C                           AND NOUT, RESPECTIVELY.                     UGE  180
C                           IF IOPT=2, THE INTERNAL VALUE OF NIN IS     UGE  190
C                           RESET FOR SUBSEQUENT USE.                   UGE  200
C                           IF IOPT=3, THE INTERNAL VALUE OF NOUT IS    UGE  210
C                           RESET FOR SUBSEQUENT USE.                   UGE  220
C                NIN    - INPUT UNIT IDENTIFIER.                        UGE  230
C                           OUTPUT IF IOPT=1, INPUT IF IOPT=2.          UGE  240
C                NOUT   - OUTPUT UNIT IDENTIFIER.                       UGE  250
C                           OUTPUT IF IOPT=1, INPUT IF IOPT=3.          UGE  260
C                                                                       UGE  270
C   PRECISION/HARDWARE  - SINGLE/ALL                                    UGE  280
C                                                                       UGE  290
C   REQD. IMSL ROUTINES - NONE REQUIRED                                 UGE  300
C                                                                       UGE  310
C   NOTATION            - INFORMATION ON SPECIAL NOTATION AND           UGE  320
C                           CONVENTIONS IS AVAILABLE IN THE MANUAL      UGE  330
C                           INTRODUCTION OR THROUGH IMSL ROUTINE UHELP  UGE  340
C                                                                       UGE  350
C   REMARKS      EACH IMSL ROUTINE THAT PERFORMS INPUT AND/OR OUTPUT    UGE  360
C                OPERATIONS CALLS UGETIO TO OBTAIN THE CURRENT UNIT     UGE  370
C                IDENTIFIER VALUES. IF UGETIO IS CALLED WITH IOPT=2 OR  UGE  380
C                IOPT=3, NEW UNIT IDENTIFIER VALUES ARE ESTABLISHED.    UGE  390
C                SUBSEQUENT INPUT/OUTPUT IS PERFORMED ON THE NEW UNITS. UGE  400
C                                                                       UGE  410
C   COPYRIGHT           - 1978 BY IMSL, INC. ALL RIGHTS RESERVED.       UGE  420
C                                                                       UGE  430
C   WARRANTY            - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN UGE  440
C                           APPLIED TO THIS CODE. NO OTHER WARRANTY,    UGE  450
C                           EXPRESSED OR IMPLIED, IS APPLICABLE.        UGE  460
C                                                                       UGE  470
C-----------------------------------------------------------------------UGE  480
C                                                                       UGE  490
      SUBROUTINE UGETIO(IOPT, NIN, NOUT)                                UGE  500
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER IOPT, NIN, NOUT
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER NIND, NOUTD
      DATA NIND /5/, NOUTD /6/
C                                  FIRST EXECUTABLE STATEMENT
      IF (IOPT.EQ.3) GO TO 20
      IF (IOPT.EQ.2) GO TO 10
      IF (IOPT.NE.1) GO TO 30
      NIN = NIND
      NOUT = NOUTD
      GO TO 30
   10 NIND = NIN
      GO TO 30
   20 NOUTD = NOUT
   30 RETURN
      END
         1         2        35     -10.0   -1.D-15                      DAT   10
         3         2        20             -1.D-15                      DAT   20
         5         2        20      -1.0   -1.D-15                      DAT   30
         5         2        20      -2.0   -1.D-15                      DAT   40
         6         2        25             -1.D-15                      DAT   50
         7         2        25             -1.D-15                      DAT   60
         8         2        25             -1.D-15                      DAT   70
         1         2        20       1.0   -1.D-15                      DAT   80
         1         2        35      10.0   -1.D-15                      DAT   90
         1         2        35      25.0   -1.D-15                      DAT  100
         4         2        25             -1.D-15                      DAT  110
         5         2        30       0.5   -1.D-15                      DAT  120
         9         2        25     100.0   -1.D-15                      DAT  130
         9         2        25      50.0   -1.D-15                      DAT  140
         9         2        25      25.0   -1.D-15                      DAT  150
         2         2        35      10.0   -1.D-15                      DAT  160
        10         2        25             -1.D-15                      DAT  170
         1         8        35     -10.0    1.D-12                      DAT  180
         3         8        20              1.D-12                      DAT  190
         5         8        20      -1.0    1.D-12                      DAT  200
         5         8        20      -2.0    1.D-12                      DAT  210
         6         8        25              1.D-12                      DAT  220
         7         8        25              1.D-12                      DAT  230
         8         8        25              1.D-12                      DAT  240
         1         8        20       1.0    1.D-12                      DAT  250
         1         8        35      10.0    1.D-12                      DAT  260
         1         8        35      25.0    1.D-12                      DAT  270
         4         8        25              1.D-12                      DAT  280
         5         8        30       0.5    1.D-12                      DAT  290
         9         8        25     100.0    1.D-12                      DAT  300
         9         8        25      50.0    1.D-12                      DAT  310
         9         8        25      25.0    1.D-12                      DAT  320
         2         8        35      10.0    1.D-12                      DAT  330
        10         8        25              1.D-12                      DAT  340
         1         6        35     -10.0    1.D-08                      DAT  350
         3         6        20              1.D-08                      DAT  360
         5         6        20      -1.0    1.D-08                      DAT  370
         5         6        20      -2.0    1.D-08                      DAT  380
         6         6        25              1.D-08                      DAT  390
         7         6        25              1.D-08                      DAT  400
         8         6        25              1.D-08                      DAT  410
         1         6        20       1.0    1.D-08                      DAT  420
         1         6        35      10.0    1.D-08                      DAT  430
         4         6        25              1.D-08                      DAT  440
         5         6        30       0.5    1.D-08                      DAT  450
         9         6        25     100.0    1.D-08                      DAT  460
         9         6        25      50.0    1.D-08                      DAT  470
         2         6        35      10.0    1.D-08                      DAT  480
        10         6        25              1.D-08                      DAT  490
         3         2        20              1.D-04                      DAT  500
         5         2        20      -1.0    1.D-04                      DAT  510
         5         2        20      -2.0    1.D-04                      DAT  520
         6         2        25              1.D-04                      DAT  530
         7         2        25              1.D-04                      DAT  540
         1         2        20       1.0    1.D-04                      DAT  550
         4         2        25              1.D-04                      DAT  560
         5         2        30       0.5    1.D-04                      DAT  570
         9         2        25     100.0    1.D-04                      DAT  580
         9         2        25      50.0    1.D-04                      DAT  590
         9         2        25      25.0    1.D-04                      DAT  600
        10         2        25              1.D-04                      DAT  610
         9         2        25     100.0   -1.D-12                      DAT  620
         9         2        25      50.0   -1.D-12                      DAT  630
         9         2        25      25.0   -1.D-12                      DAT  640
         9         2        25      50.0   -1.D-08                      DAT  650
         9         2        25      25.0   -1.D-08                      DAT  660
         9         2        25      25.0   -1.D-04                      DAT  670