[CONTACT]

[ABOUT]

[POLICY]

[ADVERTISE]

C ALGORITHM COLLECTED ALGORITHMS FROM

Found at: ftp.icm.edu.pl:70/packages/netlib/toms-2014-06-10/708

C      ALGORITHM 708, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 18, NO. 3, SEPTEMBER, 1992, PP. 360-373z.
C     PROGRAM BTST (OUTPUT, TAPE6=OUTPUT)
C-----------------------------------------------------------------------
C
C     SAMPLE PROGRAM USING BRATIO. GIVEN THE NONNEGATIVE VALUES
C     A, B, X, Y WHERE A AND B ARE NOT BOTH 0 AND X + Y = 1. THEN
C
C              CALL BRATIO (A, B, X, Y, W, W1, IERR)
C
C     COMPUTES THE VALUES
C
C                W = I (A,B)  AND W1 = 1 - I (A,B).
C                     X                     X
C
C     IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
C     IF NO INPUT ERRORS ARE DETECTED THEN IERR IS SET TO 0 AND
C     W AND W1 ARE COMPUTED. FOR MORE DETAILS SEE THE IN-LINE
C     DOCUMENTATION OF BRATIO.
C
C     THE LAST FUNCTION IN THIS PACKAGE, IPMPAR, MUST BE DEFINED
C     FOR THE PARTICULAR COMPUTER BEING USED. FOR DETAILS SEE THE
C     IN-LINE DOCUMENTATION OF IPMPAR.
C
C     NO DATA IS READ. THE OUTPUT FOR THE PROGRAM IS WRITTEN ON
C     UNIT 6. THE FIRST STATMENT OF THIS TEXT MAY BE USED TO
C     BEGIN THE PROGRAM FOR THE CDC 6000-7000 SERIES COMPUTERS.
C-----------------------------------------------------------------------
      WRITE (6,1)
    1 FORMAT(11H1   X     Y,11X,1HW,14X,2HW1/)
    2 FORMAT(2F6.2,2E16.6)
C
      A = 5.3
      B = 10.1
      X = 1.E-2
      DO 10 L = 1,50
         Y = 0.5 + (0.5 - X)
         CALL BRATIO (A, B, X, Y, W, W1, IERR)
         IF (IERR .NE. 0) STOP
         WRITE (6,2) X, Y, W, W1
         X = X + 1.E-2
   10 CONTINUE
      STOP
      END 
      SUBROUTINE BRATIO (A, B, X, Y, W, W1, IERR) 
C-----------------------------------------------------------------------
C
C            EVALUATION OF THE INCOMPLETE BETA FUNCTION IX(A,B)
C
C                     --------------------
C
C     IT IS ASSUMED THAT A AND B ARE NONNEGATIVE, AND THAT X .LE. 1
C     AND Y = 1 - X.  BRATIO ASSIGNS W AND W1 THE VALUES
C
C                      W  = IX(A,B)
C                      W1 = 1 - IX(A,B) 
C
C     IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
C     IF NO INPUT ERRORS ARE DETECTED THEN IERR IS SET TO 0 AND
C     W AND W1 ARE COMPUTED. OTHERWISE, IF AN ERROR IS DETECTED,
C     THEN W AND W1 ARE ASSIGNED THE VALUE 0 AND IERR IS SET TO
C     ONE OF THE FOLLOWING VALUES ...
C
C        IERR = 1  IF A OR B IS NEGATIVE
C        IERR = 2  IF A = B = 0
C        IERR = 3  IF X .LT. 0 OR X .GT. 1
C        IERR = 4  IF Y .LT. 0 OR Y .GT. 1
C        IERR = 5  IF X + Y .NE. 1
C        IERR = 6  IF X = A = 0
C        IERR = 7  IF Y = B = 0
C
C--------------------
C     WRITTEN BY ALFRED H. MORRIS, JR.
C        NAVAL SURFACE WARFARE CENTER
C        DAHLGREN, VIRGINIA
C     REVISED ... NOV 1991
C-----------------------------------------------------------------------
      REAL LAMBDA
C-----------------------------------------------------------------------
C
C     ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST 
C            FLOATING POINT NUMBER FOR WHICH 1.0 + EPS .GT. 1.0
C
                       EPS = SPMPAR(1)
C
C-----------------------------------------------------------------------
      W = 0.0
      W1 = 0.0
      IF (A .LT. 0.0 .OR. B .LT. 0.0) GO TO 300
      IF (A .EQ. 0.0 .AND. B .EQ. 0.0) GO TO 310
      IF (X .LT. 0.0 .OR. X .GT. 1.0) GO TO 320
      IF (Y .LT. 0.0 .OR. Y .GT. 1.0) GO TO 330
      Z = ((X + Y) - 0.5) - 0.5
      IF (ABS(Z) .GT. 3.0*EPS) GO TO 340
C
      IERR = 0
      IF (X .EQ. 0.0) GO TO 200
      IF (Y .EQ. 0.0) GO TO 210
      IF (A .EQ. 0.0) GO TO 211
      IF (B .EQ. 0.0) GO TO 201
C
      EPS = AMAX1(EPS, 1.E-15)
      IF (AMAX1(A,B) .LT. 1.E-3*EPS) GO TO 230
C
      IND = 0
      A0 = A
      B0 = B
      X0 = X
      Y0 = Y
      IF (AMIN1(A0, B0) .GT. 1.0) GO TO 30
C
C             PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1
C
      IF (X .LE. 0.5) GO TO 10
      IND = 1
      A0 = B
      B0 = A
      X0 = Y
      Y0 = X
C
   10 IF (B0 .LT. AMIN1(EPS,EPS*A0)) GO TO 80
      IF (A0 .LT. AMIN1(EPS,EPS*B0) .AND. B0*X0 .LE. 1.0) GO TO 90
      IF (AMAX1(A0, B0) .GT. 1.0) GO TO 20
      IF (A0 .GE. AMIN1(0.2, B0)) GO TO 100
      IF (X0**A0 .LE. 0.9) GO TO 100
      IF (X0 .GE. 0.3) GO TO 110
      N = 20
      GO TO 130
C
   20 IF (B0 .LE. 1.0) GO TO 100
      IF (X0 .GE. 0.3) GO TO 110
      IF (X0 .GE. 0.1) GO TO 21
      IF ((X0*B0)**A0 .LE. 0.7) GO TO 100
   21 IF (B0 .GT. 15.0) GO TO 131
      N = 20
      GO TO 130
C
C             PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1
C
   30 IF (A .GT. B) GO TO 31
         LAMBDA = A - (A + B)*X
         GO TO 32
   31 LAMBDA = (A + B)*Y - B
   32 IF (LAMBDA .GE. 0.0) GO TO 40
      IND = 1
      A0 = B
      B0 = A
      X0 = Y
      Y0 = X
      LAMBDA = ABS(LAMBDA)
C
   40 IF (B0 .LT. 40.0 .AND. B0*X0 .LE. 0.7) GO TO 100
      IF (B0 .LT. 40.0) GO TO 140
      IF (A0 .GT. B0) GO TO 50
         IF (A0 .LE. 100.0) GO TO 120
         IF (LAMBDA .GT. 0.03*A0) GO TO 120
         GO TO 180
   50 IF (B0 .LE. 100.0) GO TO 120
      IF (LAMBDA .GT. 0.03*B0) GO TO 120
      GO TO 180
C
C            EVALUATION OF THE APPROPRIATE ALGORITHM
C
   80 W = FPSER(A0, B0, X0, EPS)
      W1 = 0.5 + (0.5 - W)
      GO TO 220
C
   90 W1 = APSER(A0, B0, X0, EPS)
      W = 0.5 + (0.5 - W1)
      GO TO 220
C
  100 W = BPSER(A0, B0, X0, EPS)
      W1 = 0.5 + (0.5 - W)
      GO TO 220
C
  110 W1 = BPSER(B0, A0, Y0, EPS)
      W = 0.5 + (0.5 - W1)
      GO TO 220
C
  120 W = BFRAC(A0, B0, X0, Y0, LAMBDA, 15.0*EPS) 
      W1 = 0.5 + (0.5 - W)
      GO TO 220
C
  130 W1 = BUP(B0, A0, Y0, X0, N, EPS)
      B0 = B0 + N
  131 CALL BGRAT(B0, A0, Y0, X0, W1, 15.0*EPS, IERR1)
      W = 0.5 + (0.5 - W1)
      GO TO 220
C
  140 N = B0
      B0 = B0 - N
      IF (B0 .NE. 0.0) GO TO 141
         N = N - 1
         B0 = 1.0
  141 W = BUP(B0, A0, Y0, X0, N, EPS)
      IF (X0 .GT. 0.7) GO TO 150
      W = W + BPSER(A0, B0, X0, EPS)
      W1 = 0.5 + (0.5 - W)
      GO TO 220
C
  150 IF (A0 .GT. 15.0) GO TO 151
         N = 20
         W = W + BUP(A0, B0, X0, Y0, N, EPS)
         A0 = A0 + N
  151 CALL BGRAT(A0, B0, X0, Y0, W, 15.0*EPS, IERR1)
      W1 = 0.5 + (0.5 - W)
      GO TO 220
C
  180 W = BASYM(A0, B0, LAMBDA, 100.0*EPS)
      W1 = 0.5 + (0.5 - W)
      GO TO 220
C
C               TERMINATION OF THE PROCEDURE
C
  200 IF (A .EQ. 0.0) GO TO 350
  201 W = 0.0
      W1 = 1.0
      RETURN
C
  210 IF (B .EQ. 0.0) GO TO 360
  211 W = 1.0
      W1 = 0.0
      RETURN
C
  220 IF (IND .EQ. 0) RETURN
      T = W
      W = W1
      W1 = T
      RETURN
C
C           PROCEDURE FOR A AND B .LT. 1.E-3*EPS
C
  230 W = B/(A + B) 
      W1 = A/(A + B)
      RETURN
C
C                       ERROR RETURN
C
  300 IERR = 1
      RETURN
  310 IERR = 2
      RETURN
  320 IERR = 3
      RETURN
  330 IERR = 4
      RETURN
  340 IERR = 5
      RETURN
  350 IERR = 6
      RETURN
  360 IERR = 7
      RETURN
      END 
      REAL FUNCTION FPSER (A, B, X, EPS)
C-----------------------------------------------------------------------
C
C                 EVALUATION OF I (A,B) 
C                                X
C
C          FOR B .LT. MIN(EPS,EPS*A) AND X .LE. 0.5.
C
C-----------------------------------------------------------------------
C
C                  SET  FPSER = X**A
C
      FPSER = 1.0
      IF (A .LE. 1.E-3*EPS) GO TO 10
      FPSER = 0.0
      T = A*ALOG(X) 
      IF (T .LT. EXPARG(1)) RETURN
      FPSER = EXP(T)
C
C                NOTE THAT 1/B(A,B) = B 
C
   10 FPSER = (B/A)*FPSER
      TOL = EPS/A
      AN = A + 1.0
      T = X
      S = T/AN
   20    AN = AN + 1.0
         T = X*T
         C = T/AN
         S = S + C
         IF (ABS(C) .GT. TOL) GO TO 20
C
      FPSER = FPSER*(1.0 + A*S)
      RETURN
      END 
      REAL FUNCTION APSER (A, B, X, EPS)
C-----------------------------------------------------------------------
C     APSER YIELDS THE INCOMPLETE BETA RATIO I(SUB(1-X))(B,A) FOR
C     A .LE. MIN(EPS,EPS*B), B*X .LE. 1, AND X .LE. 0.5. USED WHEN
C     A IS VERY SMALL. USE ONLY IF ABOVE INEQUALITIES ARE SATISFIED.
C-----------------------------------------------------------------------
      REAL J
C--------------------
      DATA G/.577215664901533/
C--------------------
      BX = B*X
      T = X - BX
      IF (B*EPS .GT. 2.E-2) GO TO 10
         C = ALOG(X) + PSI(B) + G + T
         GO TO 20
   10 C = ALOG(BX) + G + T
C
   20 TOL = 5.0*EPS*ABS(C)
      J = 1.0
      S = 0.0
   30    J = J + 1.0
         T = T*(X - BX/J)
         AJ = T/J
         S = S + AJ 
         IF (ABS(AJ) .GT. TOL) GO TO 30 
C
      APSER = -A*(C + S)
      RETURN
      END 
      REAL FUNCTION BPSER(A, B, X, EPS) 
C-----------------------------------------------------------------------
C     POWER SERIES EXPANSION FOR EVALUATING IX(A,B) WHEN B .LE. 1
C     OR B*X .LE. 0.7.  EPS IS THE TOLERANCE USED.
C-----------------------------------------------------------------------
      REAL N
C
      BPSER = 0.0
      IF (X .EQ. 0.0) RETURN
C-----------------------------------------------------------------------
C            COMPUTE THE FACTOR X**A/(A*BETA(A,B))
C-----------------------------------------------------------------------
      A0 = AMIN1(A,B)
      IF (A0 .LT. 1.0) GO TO 10
         Z = A*ALOG(X) - BETALN(A,B)
         BPSER = EXP(Z)/A
         GO TO 70
   10 B0 = AMAX1(A,B)
      IF (B0 .GE. 8.0) GO TO 60
      IF (B0 .GT. 1.0) GO TO 40
C
C            PROCEDURE FOR A0 .LT. 1 AND B0 .LE. 1
C
      BPSER = X**A
      IF (BPSER .EQ. 0.0) RETURN
C
      APB = A + B
      IF (APB .GT. 1.0) GO TO 20
         Z = 1.0 + GAM1(APB)
         GO TO 30
   20 U = DBLE(A) + DBLE(B) - 1.D0
      Z = (1.0 + GAM1(U))/APB 
C
   30 C = (1.0 + GAM1(A))*(1.0 + GAM1(B))/Z
      BPSER = BPSER*C*(B/APB) 
      GO TO 70
C
C         PROCEDURE FOR A0 .LT. 1 AND 1 .LT. B0 .LT. 8
C
   40 U = GAMLN1(A0)
      M = B0 - 1.0
      IF (M .LT. 1) GO TO 50
      C = 1.0
      DO 41 I = 1,M 
         B0 = B0 - 1.0
   41    C = C*(B0/(A0 + B0)) 
      U = ALOG(C) + U
C
   50 Z = A*ALOG(X) - U
      B0 = B0 - 1.0 
      APB = A0 + B0 
      IF (APB .GT. 1.0) GO TO 51
         T = 1.0 + GAM1(APB)
         GO TO 52
   51 U = DBLE(A0) + DBLE(B0) - 1.D0
      T = (1.0 + GAM1(U))/APB 
   52 BPSER = EXP(Z)*(A0/A)*(1.0 + GAM1(B0))/T
      GO TO 70
C
C            PROCEDURE FOR A0 .LT. 1 AND B0 .GE. 8
C
   60 U = GAMLN1(A0) + ALGDIV(A0,B0)
      Z = A*ALOG(X) - U
      BPSER = (A0/A)*EXP(Z)
   70 IF (BPSER .EQ. 0.0 .OR. A .LE. 0.1*EPS) RETURN
C-----------------------------------------------------------------------
C                     COMPUTE THE SERIES
C-----------------------------------------------------------------------
      SUM = 0.0
      N = 0.0
      C = 1.0
      TOL = EPS/A
  100    N = N + 1.0
         C = C*(0.5 + (0.5 - B/N))*X
         W = C/(A + N)
         SUM = SUM + W
         IF (ABS(W) .GT. TOL) GO TO 100 
      BPSER = BPSER*(1.0 + A*SUM)
      RETURN
      END 
      REAL FUNCTION BUP(A, B, X, Y, N, EPS)
C-----------------------------------------------------------------------
C     EVALUATION OF IX(A,B) - IX(A+N,B) WHERE N IS A POSITIVE INTEGER.
C     EPS IS THE TOLERANCE USED.
C-----------------------------------------------------------------------
      REAL L
C
C          OBTAIN THE SCALING FACTOR EXP(-MU) AND 
C             EXP(MU)*(X**A*Y**B/BETA(A,B))/A
C
      APB = A + B
      AP1 = A + 1.0 
      MU = 0
      D = 1.0
      IF (N .EQ. 1 .OR. A .LT. 1.0) GO TO 10
      IF (APB .LT. 1.1*AP1) GO TO 10
         MU = ABS(EXPARG(1))
         K = EXPARG(0)
         IF (K .LT. MU) MU = K
         T = MU
         D = EXP(-T)
C
   10 BUP = BRCMP1(MU,A,B,X,Y)/A
      IF (N .EQ. 1 .OR. BUP .EQ. 0.0) RETURN
      NM1 = N - 1
      W = D
C
C          LET K BE THE INDEX OF THE MAXIMUM TERM 
C
      K = 0
      IF (B .LE. 1.0) GO TO 40
      IF (Y .GT. 1.E-4) GO TO 20
         K = NM1
         GO TO 30
   20 R = (B - 1.0)*X/Y - A
      IF (R .LT. 1.0) GO TO 40
      K = NM1
      T = NM1
      IF (R .LT. T) K = R
C
C          ADD THE INCREASING TERMS OF THE SERIES 
C
   30 DO 31 I = 1,K 
         L = I - 1
         D = ((APB + L)/(AP1 + L))*X*D
         W = W + D
   31 CONTINUE
      IF (K .EQ. NM1) GO TO 50
C
C          ADD THE REMAINING TERMS OF THE SERIES
C
   40 KP1 = K + 1
      DO 41 I = KP1,NM1
         L = I - 1
         D = ((APB + L)/(AP1 + L))*X*D
         W = W + D
         IF (D .LE. EPS*W) GO TO 50
   41 CONTINUE
C
C               TERMINATE THE PROCEDURE 
C
   50 BUP = BUP*W
      RETURN
      END 
      REAL FUNCTION BFRAC(A, B, X, Y, LAMBDA, EPS)
C-----------------------------------------------------------------------
C     CONTINUED FRACTION EXPANSION FOR IX(A,B) WHEN A,B .GT. 1.
C     IT IS ASSUMED THAT  LAMBDA = (A + B)*Y - B. 
C-----------------------------------------------------------------------
      REAL LAMBDA, N
C--------------------
      BFRAC = BRCOMP(A,B,X,Y) 
      IF (BFRAC .EQ. 0.0) RETURN
C
      C = 1.0 + LAMBDA
      C0 = B/A
      C1 = 1.0 + 1.0/A
      YP1 = Y + 1.0 
C
      N = 0.0
      P = 1.0
      S = A + 1.0
      AN = 0.0
      BN = 1.0
      ANP1 = 1.0
      BNP1 = C/C1
      R = C1/C
C
C        CONTINUED FRACTION CALCULATION 
C
   10    N = N + 1.0
         T = N/A
         W = N*(B - N)*X
         E = A/S
         ALPHA = (P*(P + C0)*E*E)*(W*X) 
         E = (1.0 + T)/(C1 + T + T)
         BETA = N + W/S + E*(C + N*YP1) 
         P = 1.0 + T
         S = S + 2.0
C
C        UPDATE AN, BN, ANP1, AND BNP1
C
         T = ALPHA*AN + BETA*ANP1
         AN = ANP1
         ANP1 = T
         T = ALPHA*BN + BETA*BNP1
         BN = BNP1
         BNP1 = T
C
         R0 = R
         R = ANP1/BNP1
         IF (ABS(R - R0) .LE. EPS*R) GO TO 20
C
C        RESCALE AN, BN, ANP1, AND BNP1 
C
         AN = AN/BNP1
         BN = BN/BNP1
         ANP1 = R
         BNP1 = 1.0 
         GO TO 10
C
C                 TERMINATION 
C
   20 BFRAC = BFRAC*R
      RETURN
      END 
      REAL FUNCTION BRCOMP (A, B, X, Y) 
C-----------------------------------------------------------------------
C               EVALUATION OF X**A*Y**B/BETA(A,B) 
C-----------------------------------------------------------------------
      REAL LAMBDA, LNX, LNY
C-----------------
C     CONST = 1/SQRT(2*PI)
C-----------------
      DATA CONST/.398942280401433/
C
      BRCOMP = 0.0
      IF (X .EQ. 0.0 .OR. Y .EQ. 0.0) RETURN
      A0 = AMIN1(A,B)
      IF (A0 .GE. 8.0) GO TO 100
C
      IF (X .GT. 0.375) GO TO 10
         LNX = ALOG(X)
         LNY = ALNREL(-X)
         GO TO 20
   10 IF (Y .GT. 0.375) GO TO 11
         LNX = ALNREL(-Y)
         LNY = ALOG(Y)
         GO TO 20
   11 LNX = ALOG(X) 
      LNY = ALOG(Y) 
C
   20 Z = A*LNX + B*LNY
      IF (A0 .LT. 1.0) GO TO 30
      Z = Z - BETALN(A,B)
      BRCOMP = EXP(Z)
      RETURN
C-----------------------------------------------------------------------
C              PROCEDURE FOR A .LT. 1 OR B .LT. 1 
C-----------------------------------------------------------------------
   30 B0 = AMAX1(A,B)
      IF (B0 .GE. 8.0) GO TO 80
      IF (B0 .GT. 1.0) GO TO 60
C
C                   ALGORITHM FOR B0 .LE. 1
C
      BRCOMP = EXP(Z)
      IF (BRCOMP .EQ. 0.0) RETURN
C
      APB = A + B
      IF (APB .GT. 1.0) GO TO 40
         Z = 1.0 + GAM1(APB)
         GO TO 50
   40 U = DBLE(A) + DBLE(B) - 1.D0
      Z = (1.0 + GAM1(U))/APB 
C
   50 C = (1.0 + GAM1(A))*(1.0 + GAM1(B))/Z
      BRCOMP = BRCOMP*(A0*C)/(1.0 + A0/B0)
      RETURN
C
C                ALGORITHM FOR 1 .LT. B0 .LT. 8
C
   60 U = GAMLN1(A0)
      N = B0 - 1.0
      IF (N .LT. 1) GO TO 70
      C = 1.0
      DO 61 I = 1,N 
         B0 = B0 - 1.0
         C = C*(B0/(A0 + B0)) 
   61 CONTINUE
      U = ALOG(C) + U
C
   70 Z = Z - U
      B0 = B0 - 1.0 
      APB = A0 + B0 
      IF (APB .GT. 1.0) GO TO 71
         T = 1.0 + GAM1(APB)
         GO TO 72
   71 U = DBLE(A0) + DBLE(B0) - 1.D0
      T = (1.0 + GAM1(U))/APB 
   72 BRCOMP = A0*EXP(Z)*(1.0 + GAM1(B0))/T
      RETURN
C
C                   ALGORITHM FOR B0 .GE. 8
C
   80 U = GAMLN1(A0) + ALGDIV(A0,B0)
      BRCOMP = A0*EXP(Z - U)
      RETURN
C-----------------------------------------------------------------------
C              PROCEDURE FOR A .GE. 8 AND B .GE. 8
C-----------------------------------------------------------------------
  100 IF (A .GT. B) GO TO 101 
         H = A/B
         X0 = H/(1.0 + H)
         Y0 = 1.0/(1.0 + H)
         LAMBDA = A - (A + B)*X
         GO TO 110
  101 H = B/A
      X0 = 1.0/(1.0 + H)
      Y0 = H/(1.0 + H)
      LAMBDA = (A + B)*Y - B
C
  110 E = -LAMBDA/A 
      IF (ABS(E) .GT. 0.6) GO TO 111
         U = RLOG1(E)
         GO TO 120
  111 U = E - ALOG(X/X0)
C
  120 E = LAMBDA/B
      IF (ABS(E) .GT. 0.6) GO TO 121
         V = RLOG1(E)
         GO TO 130
  121 V = E - ALOG(Y/Y0)
C
  130 Z = EXP(-(A*U + B*V))
      BRCOMP = CONST*SQRT(B*X0)*Z*EXP(-BCORR(A,B))
      RETURN
      END 
      REAL FUNCTION BRCMP1 (MU, A, B, X, Y)
C-----------------------------------------------------------------------
C          EVALUATION OF  EXP(MU) * (X**A*Y**B/BETA(A,B))
C-----------------------------------------------------------------------
      REAL LAMBDA, LNX, LNY
C-----------------
C     CONST = 1/SQRT(2*PI)
C-----------------
      DATA CONST/.398942280401433/
C
      A0 = AMIN1(A,B)
      IF (A0 .GE. 8.0) GO TO 100
C
      IF (X .GT. 0.375) GO TO 10
         LNX = ALOG(X)
         LNY = ALNREL(-X)
         GO TO 20
   10 IF (Y .GT. 0.375) GO TO 11
         LNX = ALNREL(-Y)
         LNY = ALOG(Y)
         GO TO 20
   11 LNX = ALOG(X) 
      LNY = ALOG(Y) 
C
   20 Z = A*LNX + B*LNY
      IF (A0 .LT. 1.0) GO TO 30
      Z = Z - BETALN(A,B)
      BRCMP1 = ESUM(MU,Z)
      RETURN
C-----------------------------------------------------------------------
C              PROCEDURE FOR A .LT. 1 OR B .LT. 1 
C-----------------------------------------------------------------------
   30 B0 = AMAX1(A,B)
      IF (B0 .GE. 8.0) GO TO 80
      IF (B0 .GT. 1.0) GO TO 60
C
C                   ALGORITHM FOR B0 .LE. 1
C
      BRCMP1 = ESUM(MU,Z)
      IF (BRCMP1 .EQ. 0.0) RETURN
C
      APB = A + B
      IF (APB .GT. 1.0) GO TO 40
         Z = 1.0 + GAM1(APB)
         GO TO 50
   40 U = DBLE(A) + DBLE(B) - 1.D0
      Z = (1.0 + GAM1(U))/APB 
C
   50 C = (1.0 + GAM1(A))*(1.0 + GAM1(B))/Z
      BRCMP1 = BRCMP1*(A0*C)/(1.0 + A0/B0)
      RETURN
C
C                ALGORITHM FOR 1 .LT. B0 .LT. 8
C
   60 U = GAMLN1(A0)
      N = B0 - 1.0
      IF (N .LT. 1) GO TO 70
      C = 1.0
      DO 61 I = 1,N 
         B0 = B0 - 1.0
         C = C*(B0/(A0 + B0)) 
   61 CONTINUE
      U = ALOG(C) + U
C
   70 Z = Z - U
      B0 = B0 - 1.0 
      APB = A0 + B0 
      IF (APB .GT. 1.0) GO TO 71
         T = 1.0 + GAM1(APB)
         GO TO 72
   71 U = DBLE(A0) + DBLE(B0) - 1.D0
      T = (1.0 + GAM1(U))/APB 
   72 BRCMP1 = A0*ESUM(MU,Z)*(1.0 + GAM1(B0))/T
      RETURN
C
C                   ALGORITHM FOR B0 .GE. 8
C
   80 U = GAMLN1(A0) + ALGDIV(A0,B0)
      BRCMP1 = A0*ESUM(MU,Z - U)
      RETURN
C-----------------------------------------------------------------------
C              PROCEDURE FOR A .GE. 8 AND B .GE. 8
C-----------------------------------------------------------------------
  100 IF (A .GT. B) GO TO 101 
         H = A/B
         X0 = H/(1.0 + H)
         Y0 = 1.0/(1.0 + H)
         LAMBDA = A - (A + B)*X
         GO TO 110
  101 H = B/A
      X0 = 1.0/(1.0 + H)
      Y0 = H/(1.0 + H)
      LAMBDA = (A + B)*Y - B
C
  110 E = -LAMBDA/A 
      IF (ABS(E) .GT. 0.6) GO TO 111
         U = RLOG1(E)
         GO TO 120
  111 U = E - ALOG(X/X0)
C
  120 E = LAMBDA/B
      IF (ABS(E) .GT. 0.6) GO TO 121
         V = RLOG1(E)
         GO TO 130
  121 V = E - ALOG(Y/Y0)
C
  130 Z = ESUM(MU,-(A*U + B*V))
      BRCMP1 = CONST*SQRT(B*X0)*Z*EXP(-BCORR(A,B))
      RETURN
      END 
      SUBROUTINE BGRAT(A, B, X, Y, W, EPS, IERR)
C-----------------------------------------------------------------------
C     ASYMPTOTIC EXPANSION FOR IX(A,B) WHEN A IS LARGER THAN B.
C     THE RESULT OF THE EXPANSION IS ADDED TO W. IT IS ASSUMED
C     THAT A .GE. 15 AND B .LE. 1.  EPS IS THE TOLERANCE USED.
C     IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
C-----------------------------------------------------------------------
      REAL J, L, LNX, NU, N2
      REAL C(30), D(30)
C
      BM1 = (B - 0.5) - 0.5
      NU = A + 0.5*BM1
      IF (Y .GT. 0.375) GO TO 10
         LNX = ALNREL(-Y)
         GO TO 11
   10 LNX = ALOG(X) 
   11 Z = -NU*LNX
      IF (B*Z .EQ. 0.0) GO TO 100
C
C                 COMPUTATION OF THE EXPANSION
C                 SET R = EXP(-Z)*Z**B/GAMMA(B)
C
      R = B*(1.0 + GAM1(B))*EXP(B*ALOG(Z))
      R = R*EXP(A*LNX)*EXP(0.5*BM1*LNX) 
      U = ALGDIV(B,A) + B*ALOG(NU)
      U = R*EXP(-U) 
      IF (U .EQ. 0.0) GO TO 100
      CALL GRAT1(B,Z,R,P,Q,EPS)
C
      V = 0.25*(1.0/NU)**2
      T2 = 0.25*LNX*LNX
      L = W/U
      J = Q/R
      SUM = J
      T = 1.0
      CN = 1.0
      N2 = 0.0
      DO 22 N = 1,30
         BP2N = B + N2
         J = (BP2N*(BP2N + 1.0)*J + (Z + BP2N + 1.0)*T)*V
         N2 = N2 + 2.0
         T = T*T2
         CN = CN/(N2*(N2 + 1.0))
         C(N) = CN
         S = 0.0
         IF (N .EQ. 1) GO TO 21
            NM1 = N - 1
            COEF = B - N
            DO 20 I = 1,NM1
               S = S + COEF*C(I)*D(N-I) 
   20          COEF = COEF + B
   21    D(N) = BM1*CN + S/N
         DJ = D(N)*J
         SUM = SUM + DJ
         IF (SUM .LE. 0.0) GO TO 100
         IF (ABS(DJ) .LE. EPS*(SUM + L)) GO TO 30 
   22 CONTINUE
C
C                    ADD THE RESULTS TO W
C
   30 IERR = 0
      W = W + U*SUM 
      RETURN
C
C               THE EXPANSION CANNOT BE COMPUTED
C
  100 IERR = 1
      RETURN
      END 
      SUBROUTINE GRAT1 (A,X,R,P,Q,EPS)
      REAL J, L
C-----------------------------------------------------------------------
C        EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS 
C                      P(A,X) AND Q(A,X)
C
C     IT IS ASSUMED THAT A .LE. 1.  EPS IS THE TOLERANCE TO BE USED.
C     THE INPUT ARGUMENT R HAS THE VALUE E**(-X)*X**A/GAMMA(A).
C-----------------------------------------------------------------------
      IF (A*X .EQ. 0.0) GO TO 130
      IF (A .EQ. 0.5) GO TO 120
      IF (X .LT. 1.1) GO TO 10
      GO TO 50
C
C             TAYLOR SERIES FOR P(A,X)/X**A
C
   10 AN = 3.0
      C = X
      SUM = X/(A + 3.0)
      TOL = 0.1*EPS/(A + 1.0) 
   11    AN = AN + 1.0
         C = -C*(X/AN)
         T = C/(A + AN)
         SUM = SUM + T
         IF (ABS(T) .GT. TOL) GO TO 11
      J = A*X*((SUM/6.0 - 0.5/(A + 2.0))*X + 1.0/(A + 1.0)) 
C
      Z = A*ALOG(X) 
      H = GAM1(A)
      G = 1.0 + H
      IF (X .LT. 0.25) GO TO 20
         IF (A .LT. X/2.59) GO TO 40
         GO TO 30
   20 IF (Z .GT. -.13394) GO TO 40
C
   30 W = EXP(Z)
      P = W*G*(0.5 + (0.5 - J))
      Q = 0.5 + (0.5 - P)
      RETURN
C
   40 L = REXP(Z)
      W = 0.5 + (0.5 + L)
      Q = (W*J - L)*G - H
      IF (Q .LT. 0.0) GO TO 110
      P = 0.5 + (0.5 - Q)
      RETURN
C
C              CONTINUED FRACTION EXPANSION
C
   50 A2NM1 = 1.0
      A2N = 1.0
      B2NM1 = X
      B2N = X + (1.0 - A)
      C = 1.0
   51    A2NM1 = X*A2N + C*A2NM1
         B2NM1 = X*B2N + C*B2NM1
         AM0 = A2NM1/B2NM1
         C = C + 1.0
         CMA = C - A
         A2N = A2NM1 + CMA*A2N
         B2N = B2NM1 + CMA*B2N
         AN0 = A2N/B2N
         IF (ABS(AN0 - AM0) .GE. EPS*AN0) GO TO 51
      Q = R*AN0
      P = 0.5 + (0.5 - Q)
      RETURN
C
C                SPECIAL CASES
C
  100 P = 0.0
      Q = 1.0
      RETURN
C
  110 P = 1.0
      Q = 0.0
      RETURN
C
  120 IF (X .GE. 0.25) GO TO 121
      P = ERF(SQRT(X))
      Q = 0.5 + (0.5 - P)
      RETURN
  121 Q = ERFC1(0,SQRT(X))
      P = 0.5 + (0.5 - Q)
      RETURN
C
  130 IF (X .LE. A) GO TO 100 
      GO TO 110
      END 
      REAL FUNCTION BASYM(A, B, LAMBDA, EPS)
C-----------------------------------------------------------------------
C     ASYMPTOTIC EXPANSION FOR IX(A,B) FOR LARGE A AND B.
C     LAMBDA = (A + B)*Y - B  AND EPS IS THE TOLERANCE USED.
C     IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT
C     A AND B ARE GREATER THAN OR EQUAL TO 15.
C-----------------------------------------------------------------------
      REAL J0, J1, LAMBDA
      REAL A0(21), B0(21), C(21), D(21) 
C------------------------
C     ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP
C            ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN. 
C            THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1.
C
                      DATA NUM/20/
C------------------------
C     E0 = 2/SQRT(PI)
C     E1 = 2**(-3/2)
C------------------------
      DATA E0/1.12837916709551/, E1/.353553390593274/
C------------------------
      BASYM = 0.0
      IF (A .GE. B) GO TO 10
         H = A/B
         R0 = 1.0/(1.0 + H)
         R1 = (B - A)/B
         W0 = 1.0/SQRT(A*(1.0 + H))
         GO TO 20
   10 H = B/A
      R0 = 1.0/(1.0 + H)
      R1 = (B - A)/A
      W0 = 1.0/SQRT(B*(1.0 + H))
C
   20 F = A*RLOG1(-LAMBDA/A) + B*RLOG1(LAMBDA/B)
      T = EXP(-F)
      IF (T .EQ. 0.0) RETURN
      Z0 = SQRT(F)
      Z = 0.5*(Z0/E1)
      Z2 = F + F
C
      A0(1) = (2.0/3.0)*R1
      C(1) = - 0.5*A0(1)
      D(1) = - C(1) 
      J0 = (0.5/E0)*ERFC1(1,Z0)
      J1 = E1
      SUM = J0 + D(1)*W0*J1
C
      S = 1.0
      H2 = H*H
      HN = 1.0
      W = W0
      ZNM1 = Z
      ZN = Z2
      DO 50 N = 2, NUM, 2
         HN = H2*HN 
         A0(N) = 2.0*R0*(1.0 + H*HN)/(N + 2.0)
         NP1 = N + 1
         S = S + HN 
         A0(NP1) = 2.0*R1*S/(N + 3.0)
C
         DO 41 I = N, NP1
         R = -0.5*(I + 1.0)
         B0(1) = R*A0(1)
         DO 31 M = 2, I
            BSUM = 0.0
            MM1 = M - 1
            DO 30 J = 1, MM1
               MMJ = M - J
   30          BSUM = BSUM + (J*R - MMJ)*A0(J)*B0(MMJ)
   31       B0(M) = R*A0(M) + BSUM/M
         C(I) = B0(I)/(I + 1.0)
C
         DSUM = 0.0 
         IM1 = I - 1
         DO 40 J = 1, IM1
            IMJ = I - J
   40       DSUM = DSUM + D(IMJ)*C(J)
   41    D(I) = -(DSUM + C(I))
C
         J0 = E1*ZNM1 + (N - 1.0)*J0
         J1 = E1*ZN + N*J1
         ZNM1 = Z2*ZNM1
         ZN = Z2*ZN 
         W = W0*W
         T0 = D(N)*W*J0
         W = W0*W
         T1 = D(NP1)*W*J1
         SUM = SUM + (T0 + T1)
         IF ((ABS(T0) + ABS(T1)) .LE. EPS*SUM) GO TO 60
   50    CONTINUE
C
   60 U = EXP(-BCORR(A,B))
      BASYM = E0*T*U*SUM
      RETURN
      END 
      REAL FUNCTION SPMPAR (I)
C-----------------------------------------------------------------------
C
C     SPMPAR PROVIDES THE SINGLE PRECISION MACHINE CONSTANTS FOR
C     THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT
C     I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE
C     SINGLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND
C     ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN
C
C        SPMPAR(1) = B**(1 - M), THE MACHINE PRECISION,
C
C        SPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE, 
C
C        SPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE.
C
C-----------------------------------------------------------------------
C     WRITTEN BY
C        ALFRED H. MORRIS, JR.
C        NAVAL SURFACE WARFARE CENTER
C        DAHLGREN VIRGINIA
C-----------------------------------------------------------------------
      INTEGER EMIN, EMAX
C
      IF (I .GT. 1) GO TO 10
         B = IPMPAR(4)
         M = IPMPAR(5)
         SPMPAR = B**(1 - M)
         RETURN
C
   10 IF (I .GT. 2) GO TO 20
         B = IPMPAR(4)
         EMIN = IPMPAR(6)
         ONE = FLOAT(1)
         BINV = ONE/B
         W = B**(EMIN + 2)
         SPMPAR = ((W * BINV) * BINV) * BINV
         RETURN
C
   20 IBETA = IPMPAR(4)
      M = IPMPAR(5) 
      EMAX = IPMPAR(7)
C
      B = IBETA
      BM1 = IBETA - 1
      ONE = FLOAT(1)
      Z = B**(M - 1)
      W = ((Z - ONE)*B + BM1)/(B*Z)
C
      Z = B**(EMAX - 2)
      SPMPAR = ((W * Z) * B) * B
      RETURN
      END 
      REAL FUNCTION EXPARG (L)
C-------------------------------------------------------------------- 
C     IF L = 0 THEN  EXPARG(L) = THE LARGEST POSITIVE W FOR WHICH
C     EXP(W) CAN BE COMPUTED. 
C
C     IF L IS NONZERO THEN  EXPARG(L) = THE LARGEST NEGATIVE W FOR
C     WHICH THE COMPUTED VALUE OF EXP(W) IS NONZERO.
C
C     NOTE... ONLY AN APPROXIMATE VALUE FOR EXPARG(L) IS NEEDED.
C-------------------------------------------------------------------- 
      INTEGER B
      REAL LNB
C
      B = IPMPAR(4) 
      IF (B .NE. 2) GO TO 10
         LNB = .69314718055995
         GO TO 50
   10 IF (B .NE. 8) GO TO 20
         LNB = 2.0794415416798
         GO TO 50
   20 IF (B .NE. 16) GO TO 30 
         LNB = 2.7725887222398
         GO TO 50
   30 LNB = ALOG(FLOAT(B))
C
   50 IF (L .EQ. 0) GO TO 60
         M = IPMPAR(6) - 1
         EXPARG = 0.99999 * (M * LNB)
         RETURN
   60 M = IPMPAR(7) 
      EXPARG = 0.99999 * (M * LNB)
      RETURN
      END 
      REAL FUNCTION ESUM (MU, X)
C-----------------------------------------------------------------------
C                    EVALUATION OF EXP(MU + X)
C-----------------------------------------------------------------------
      IF (X .GT. 0.0) GO TO 10
C
      IF (MU .LT. 0) GO TO 20 
         W = MU + X 
         IF (W .GT. 0.0) GO TO 20
         ESUM = EXP(W)
         RETURN
C
   10 IF (MU .GT. 0) GO TO 20 
         W = MU + X 
         IF (W .LT. 0.0) GO TO 20
         ESUM = EXP(W)
         RETURN
C
   20 W = MU
      ESUM = EXP(W)*EXP(X)
      RETURN
      END 
      REAL FUNCTION REXP (X)
C-----------------------------------------------------------------------
C            EVALUATION OF THE FUNCTION EXP(X) - 1
C-----------------------------------------------------------------------
      DATA P1/ .914041914819518E-09/, P2/ .238082361044469E-01/,
     *     Q1/-.499999999085958E+00/, Q2/ .107141568980644E+00/,
     *     Q3/-.119041179760821E-01/, Q4/ .595130811860248E-03/
C-----------------------
      IF (ABS(X) .GT. 0.15) GO TO 10
      REXP = X*(((P2*X + P1)*X + 1.0)/((((Q4*X + Q3)*X + Q2)*X
     *                 + Q1)*X + 1.0))
      RETURN
C
   10 W = EXP(X)
      IF (X .GT. 0.0) GO TO 20
         REXP = (W - 0.5) - 0.5
         RETURN
   20 REXP = W*(0.5 + (0.5 - 1.0/W))
      RETURN
      END 
      REAL FUNCTION ALNREL(A) 
C-----------------------------------------------------------------------
C            EVALUATION OF THE FUNCTION LN(1 + A) 
C-----------------------------------------------------------------------
      DATA P1/-.129418923021993E+01/, P2/.405303492862024E+00/,
     *     P3/-.178874546012214E-01/
      DATA Q1/-.162752256355323E+01/, Q2/.747811014037616E+00/,
     *     Q3/-.845104217945565E-01/
C--------------------------
      IF (ABS(A) .GT. 0.375) GO TO 10
      T = A/(A + 2.0)
      T2 = T*T
      W = (((P3*T2 + P2)*T2 + P1)*T2 + 1.0)/
     *    (((Q3*T2 + Q2)*T2 + Q1)*T2 + 1.0)
      ALNREL = 2.0*T*W
      RETURN
C
   10 X = 1.D0 + DBLE(A)
      ALNREL = ALOG(X)
      RETURN
      END 
      REAL FUNCTION RLOG1(X)
C-----------------------------------------------------------------------
C             EVALUATION OF THE FUNCTION X - LN(1 + X)
C-----------------------------------------------------------------------
      DATA A/.566749439387324E-01/
      DATA B/.456512608815524E-01/
C------------------------
      DATA P0/ .333333333333333E+00/, P1/-.224696413112536E+00/,
     *     P2/ .620886815375787E-02/
      DATA Q1/-.127408923933623E+01/, Q2/ .354508718369557E+00/
C------------------------
      IF (X .LT. -0.39 .OR. X .GT. 0.57) GO TO 100
      IF (X .LT. -0.18) GO TO 10
      IF (X .GT.  0.18) GO TO 20
C
C              ARGUMENT REDUCTION
C
      H = X
      W1 = 0.0
      GO TO 30
C
   10 H = DBLE(X) + 0.3D0
      H = H/0.7
      W1 = A - H*0.3
      GO TO 30
C
   20 H = 0.75D0*DBLE(X) - 0.25D0
      W1 = B + H/3.0
C
C               SERIES EXPANSION
C
   30 R = H/(H + 2.0)
      T = R*R
      W = ((P2*T + P1)*T + P0)/((Q2*T + Q1)*T + 1.0)
      RLOG1 = 2.0*T*(1.0/(1.0 - R) - R*W) + W1
      RETURN
C
C
  100 W = (X + 0.5) + 0.5
      RLOG1 = X - ALOG(W)
      RETURN
      END 
      REAL FUNCTION ERF (X)
C-----------------------------------------------------------------------
C             EVALUATION OF THE REAL ERROR FUNCTION
C-----------------------------------------------------------------------
      REAL A(5),B(3),P(8),Q(8),R(5),S(4)
C-------------------------
      DATA C /.564189583547756/
C-------------------------
      DATA A(1) /.771058495001320E-04/, A(2)/-.133733772997339E-02/,
     *     A(3) /.323076579225834E-01/, A(4) /.479137145607681E-01/,
     *     A(5) /.128379167095513E+00/
      DATA B(1) /.301048631703895E-02/, B(2) /.538971687740286E-01/,
     *     B(3) /.375795757275549E+00/
C-------------------------
      DATA P(1)/-1.36864857382717E-07/, P(2) /5.64195517478974E-01/,
     *     P(3) /7.21175825088309E+00/, P(4) /4.31622272220567E+01/,
     *     P(5) /1.52989285046940E+02/, P(6) /3.39320816734344E+02/,
     *     P(7) /4.51918953711873E+02/, P(8) /3.00459261020162E+02/
      DATA Q(1) /1.00000000000000E+00/, Q(2) /1.27827273196294E+01/,
     *     Q(3) /7.70001529352295E+01/, Q(4) /2.77585444743988E+02/,
     *     Q(5) /6.38980264465631E+02/, Q(6) /9.31354094850610E+02/,
     *     Q(7) /7.90950925327898E+02/, Q(8) /3.00459260956983E+02/
C-------------------------
      DATA R(1) /2.10144126479064E+00/, R(2) /2.62370141675169E+01/,
     *     R(3) /2.13688200555087E+01/, R(4) /4.65807828718470E+00/,
     *     R(5) /2.82094791773523E-01/
      DATA S(1) /9.41537750555460E+01/, S(2) /1.87114811799590E+02/,
     *     S(3) /9.90191814623914E+01/, S(4) /1.80124575948747E+01/
C-------------------------
      AX = ABS(X)
      IF (AX .GT. 0.5) GO TO 10
      T = X*X
      TOP = ((((A(1)*T + A(2))*T + A(3))*T + A(4))*T + A(5)) + 1.0
      BOT = ((B(1)*T + B(2))*T + B(3))*T + 1.0
      ERF = X*(TOP/BOT)
      RETURN
C
   10 IF (AX .GT. 4.0) GO TO 20
      TOP = ((((((P(1)*AX + P(2))*AX + P(3))*AX + P(4))*AX + P(5))*AX 
     *                    + P(6))*AX + P(7))*AX + P(8)
      BOT = ((((((Q(1)*AX + Q(2))*AX + Q(3))*AX + Q(4))*AX + Q(5))*AX 
     *                    + Q(6))*AX + Q(7))*AX + Q(8)
      ERF = 0.5 + (0.5 - EXP(-X*X)*TOP/BOT)
      IF (X .LT. 0.0) ERF = -ERF
      RETURN
C
   20 IF (AX .GE. 5.8) GO TO 30
      X2 = X*X
      T = 1.0/X2
      TOP = (((R(1)*T + R(2))*T + R(3))*T + R(4))*T + R(5)
      BOT = (((S(1)*T + S(2))*T + S(3))*T + S(4))*T + 1.0
      ERF = (C - TOP/(X2*BOT)) / AX
      ERF = 0.5 + (0.5 - EXP(-X2)*ERF)
      IF (X .LT. 0.0) ERF = -ERF
      RETURN
C
   30 ERF = SIGN(1.0,X)
      RETURN
      END 
      REAL FUNCTION ERFC1 (IND, X)
C-----------------------------------------------------------------------
C         EVALUATION OF THE COMPLEMENTARY ERROR FUNCTION
C
C          ERFC1(IND,X) = ERFC(X)            IF IND = 0
C          ERFC1(IND,X) = EXP(X*X)*ERFC(X)   OTHERWISE
C-----------------------------------------------------------------------
      REAL A(5),B(3),P(8),Q(8),R(5),S(4)
      DOUBLE PRECISION W
C-------------------------
      DATA C /.564189583547756/
C-------------------------
      DATA A(1) /.771058495001320E-04/, A(2)/-.133733772997339E-02/,
     *     A(3) /.323076579225834E-01/, A(4) /.479137145607681E-01/,
     *     A(5) /.128379167095513E+00/
      DATA B(1) /.301048631703895E-02/, B(2) /.538971687740286E-01/,
     *     B(3) /.375795757275549E+00/
C-------------------------
      DATA P(1)/-1.36864857382717E-07/, P(2) /5.64195517478974E-01/,
     *     P(3) /7.21175825088309E+00/, P(4) /4.31622272220567E+01/,
     *     P(5) /1.52989285046940E+02/, P(6) /3.39320816734344E+02/,
     *     P(7) /4.51918953711873E+02/, P(8) /3.00459261020162E+02/
      DATA Q(1) /1.00000000000000E+00/, Q(2) /1.27827273196294E+01/,
     *     Q(3) /7.70001529352295E+01/, Q(4) /2.77585444743988E+02/,
     *     Q(5) /6.38980264465631E+02/, Q(6) /9.31354094850610E+02/,
     *     Q(7) /7.90950925327898E+02/, Q(8) /3.00459260956983E+02/
C-------------------------
      DATA R(1) /2.10144126479064E+00/, R(2) /2.62370141675169E+01/,
     *     R(3) /2.13688200555087E+01/, R(4) /4.65807828718470E+00/,
     *     R(5) /2.82094791773523E-01/
      DATA S(1) /9.41537750555460E+01/, S(2) /1.87114811799590E+02/,
     *     S(3) /9.90191814623914E+01/, S(4) /1.80124575948747E+01/
C-------------------------
C
C                     ABS(X) .LE. 0.5
C
      AX = ABS(X)
      IF (AX .GT. 0.5) GO TO 10
      T = X*X
      TOP = ((((A(1)*T + A(2))*T + A(3))*T + A(4))*T + A(5)) + 1.0
      BOT = ((B(1)*T + B(2))*T + B(3))*T + 1.0
      ERFC1 = 0.5 + (0.5 - X*(TOP/BOT)) 
      IF (IND .NE. 0) ERFC1 = EXP(T) * ERFC1
      RETURN
C
C                  0.5 .LT. ABS(X) .LE. 4
C
   10 IF (AX .GT. 4.0) GO TO 20
      TOP = ((((((P(1)*AX + P(2))*AX + P(3))*AX + P(4))*AX + P(5))*AX 
     *                    + P(6))*AX + P(7))*AX + P(8)
      BOT = ((((((Q(1)*AX + Q(2))*AX + Q(3))*AX + Q(4))*AX + Q(5))*AX 
     *                    + Q(6))*AX + Q(7))*AX + Q(8)
      ERFC1 = TOP/BOT
      GO TO 40
C
C                      ABS(X) .GT. 4
C
   20 IF (X .LE. -5.6) GO TO 50
      IF (IND .NE. 0) GO TO 30
      IF (X .GT. 100.0) GO TO 60
      IF (X*X .GT. -EXPARG(1)) GO TO 60 
C
   30 T = (1.0/X)**2
      TOP = (((R(1)*T + R(2))*T + R(3))*T + R(4))*T + R(5)
      BOT = (((S(1)*T + S(2))*T + S(3))*T + S(4))*T + 1.0
      ERFC1 = (C - T*TOP/BOT)/AX
C
C                      FINAL ASSEMBLY
C
   40 IF (IND .EQ. 0) GO TO 41
         IF (X .LT. 0.0) ERFC1 = 2.0*EXP(X*X) - ERFC1
         RETURN
   41 W = DBLE(X)*DBLE(X)
      T = W
      E = W - DBLE(T)
      ERFC1 = ((0.5 + (0.5 - E)) * EXP(-T)) * ERFC1
      IF (X .LT. 0.0) ERFC1 = 2.0 - ERFC1
      RETURN
C
C             LIMIT VALUE FOR LARGE NEGATIVE X
C
   50 ERFC1 = 2.0
      IF (IND .NE. 0) ERFC1 = 2.0*EXP(X*X)
      RETURN
C
C             LIMIT VALUE FOR LARGE POSITIVE X
C                       WHEN IND = 0
C
   60 ERFC1 = 0.0
      RETURN
      END 
      REAL FUNCTION GAM1(A)
C     ------------------------------------------------------------------
C     COMPUTATION OF 1/GAMMA(A+1) - 1  FOR -0.5 .LE. A .LE. 1.5
C     ------------------------------------------------------------------
      REAL P(7), Q(5), R(9)
C     -------------------
      DATA P(1)/ .577215664901533E+00/, P(2)/-.409078193005776E+00/,
     *     P(3)/-.230975380857675E+00/, P(4)/ .597275330452234E-01/,
     *     P(5)/ .766968181649490E-02/, P(6)/-.514889771323592E-02/,
     *     P(7)/ .589597428611429E-03/
C     -------------------
      DATA Q(1)/ .100000000000000E+01/, Q(2)/ .427569613095214E+00/,
     *     Q(3)/ .158451672430138E+00/, Q(4)/ .261132021441447E-01/,
     *     Q(5)/ .423244297896961E-02/
C     -------------------
      DATA R(1)/-.422784335098468E+00/, R(2)/-.771330383816272E+00/,
     *     R(3)/-.244757765222226E+00/, R(4)/ .118378989872749E+00/,
     *     R(5)/ .930357293360349E-03/, R(6)/-.118290993445146E-01/,
     *     R(7)/ .223047661158249E-02/, R(8)/ .266505979058923E-03/,
     *     R(9)/-.132674909766242E-03/
C     -------------------
      DATA S1  / .273076135303957E+00/, S2  / .559398236957378E-01/
C     -------------------
      T = A
      D = A - 0.5
      IF (D .GT. 0.0) T = D - 0.5
      IF (T) 30,10,20
C
   10 GAM1 = 0.0
      RETURN
C
   20 TOP = (((((P(7)*T + P(6))*T + P(5))*T + P(4))*T + P(3))*T
     *                  + P(2))*T + P(1)
      BOT = (((Q(5)*T + Q(4))*T + Q(3))*T + Q(2))*T + 1.0
      W = TOP/BOT
      IF (D .GT. 0.0) GO TO 21
         GAM1 = A*W 
         RETURN
   21 GAM1 = (T/A)*((W - 0.5) - 0.5)
      RETURN
C
   30 TOP = (((((((R(9)*T + R(8))*T + R(7))*T + R(6))*T + R(5))*T
     *                    + R(4))*T + R(3))*T + R(2))*T + R(1)
      BOT = (S2*T + S1)*T + 1.0
      W = TOP/BOT
      IF (D .GT. 0.0) GO TO 31
         GAM1 = A*((W + 0.5) + 0.5)
         RETURN
   31 GAM1 = T*W/A
      RETURN
      END 
      REAL FUNCTION GAMLN1 (A)
C-----------------------------------------------------------------------
C     EVALUATION OF LN(GAMMA(1 + A)) FOR -0.2 .LE. A .LE. 1.25
C-----------------------------------------------------------------------
      DATA P0/ .577215664901533E+00/, P1/ .844203922187225E+00/,
     *     P2/-.168860593646662E+00/, P3/-.780427615533591E+00/,
     *     P4/-.402055799310489E+00/, P5/-.673562214325671E-01/,
     *     P6/-.271935708322958E-02/
      DATA Q1/ .288743195473681E+01/, Q2/ .312755088914843E+01/,
     *     Q3/ .156875193295039E+01/, Q4/ .361951990101499E+00/,
     *     Q5/ .325038868253937E-01/, Q6/ .667465618796164E-03/
C----------------------
      DATA R0/.422784335098467E+00/,  R1/.848044614534529E+00/,
     *     R2/.565221050691933E+00/,  R3/.156513060486551E+00/,
     *     R4/.170502484022650E-01/,  R5/.497958207639485E-03/
      DATA S1/.124313399877507E+01/,  S2/.548042109832463E+00/,
     *     S3/.101552187439830E+00/,  S4/.713309612391000E-02/,
     *     S5/.116165475989616E-03/
C----------------------
      IF (A .GE. 0.6) GO TO 10
      W = ((((((P6*A + P5)*A + P4)*A + P3)*A + P2)*A + P1)*A + P0)/
     *    ((((((Q6*A + Q5)*A + Q4)*A + Q3)*A + Q2)*A + Q1)*A + 1.0)
      GAMLN1 = -A*W 
      RETURN
C
   10 X = (A - 0.5) - 0.5
      W = (((((R5*X + R4)*X + R3)*X + R2)*X + R1)*X + R0)/
     *    (((((S5*X + S4)*X + S3)*X + S2)*X + S1)*X + 1.0)
      GAMLN1 = X*W
      RETURN
      END 
      REAL FUNCTION PSI(XX)
C---------------------------------------------------------------------
C
C                 EVALUATION OF THE DIGAMMA FUNCTION
C
C                           ----------- 
C
C     PSI(XX) IS ASSIGNED THE VALUE 0 WHEN THE DIGAMMA FUNCTION CANNOT
C     BE COMPUTED.
C
C     THE MAIN COMPUTATION INVOLVES EVALUATION OF RATIONAL CHEBYSHEV
C     APPROXIMATIONS PUBLISHED IN MATH. COMP. 27, 123-127(1973) BY
C     CODY, STRECOK AND THACHER.
C
C---------------------------------------------------------------------
C     PSI WAS WRITTEN AT ARGONNE NATIONAL LABORATORY FOR THE FUNPACK
C     PACKAGE OF SPECIAL FUNCTION SUBROUTINES. PSI WAS MODIFIED BY
C     A.H. MORRIS (NSWC).
C---------------------------------------------------------------------
      REAL P1(7), P2(4), Q1(6), Q2(4)
      DOUBLE PRECISION DX0
C---------------------------------------------------------------------
C
C     PIOV4 = PI/4
C     DX0 = ZERO OF PSI TO EXTENDED PRECISION
C
C---------------------------------------------------------------------
      DATA PIOV4/.785398163397448E0/
      DATA DX0/1.461632144968362341262659542325721325D0/
C---------------------------------------------------------------------
C
C     COEFFICIENTS FOR RATIONAL APPROXIMATION OF
C     PSI(X) / (X - X0),  0.5 .LE. X .LE. 3.0
C
C---------------------------------------------------------------------
      DATA P1(1)/.895385022981970E-02/,  P1(2)/.477762828042627E+01/, 
     *     P1(3)/.142441585084029E+03/,  P1(4)/.118645200713425E+04/, 
     *     P1(5)/.363351846806499E+04/,  P1(6)/.413810161269013E+04/, 
     *     P1(7)/.130560269827897E+04/
      DATA Q1(1)/.448452573429826E+02/,  Q1(2)/.520752771467162E+03/, 
     *     Q1(3)/.221000799247830E+04/,  Q1(4)/.364127349079381E+04/, 
     *     Q1(5)/.190831076596300E+04/,  Q1(6)/.691091682714533E-05/
C---------------------------------------------------------------------
C
C     COEFFICIENTS FOR RATIONAL APPROXIMATION OF
C     PSI(X) - LN(X) + 1 / (2*X),  X .GT. 3.0
C
C---------------------------------------------------------------------
      DATA P2(1)/-.212940445131011E+01/, P2(2)/-.701677227766759E+01/,
     *     P2(3)/-.448616543918019E+01/, P2(4)/-.648157123766197E+00/ 
      DATA Q2(1)/ .322703493791143E+02/, Q2(2)/ .892920700481861E+02/,
     *     Q2(3)/ .546117738103215E+02/, Q2(4)/ .777788548522962E+01/ 
C---------------------------------------------------------------------
C
C     MACHINE DEPENDENT CONSTANTS ...
C
C        XMAX1  = THE SMALLEST POSITIVE FLOATING POINT CONSTANT
C                 WITH ENTIRELY INTEGER REPRESENTATION.  ALSO USED
C                 AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE
C                 ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH 
C                 PSI MAY BE REPRESENTED AS ALOG(X).
C
C        XSMALL = ABSOLUTE ARGUMENT BELOW WHICH PI*COTAN(PI*X)
C                 MAY BE REPRESENTED BY 1/X.
C
C---------------------------------------------------------------------
      XMAX1 = IPMPAR(3)
      XMAX1 = AMIN1(XMAX1, 1.0/SPMPAR(1))
      XSMALL = 1.E-9
C---------------------------------------------------------------------
      X = XX
      AUG = 0.0E0
      IF (X .GE. 0.5E0) GO TO 200
C---------------------------------------------------------------------
C     X .LT. 0.5,  USE REFLECTION FORMULA
C     PSI(1-X) = PSI(X) + PI * COTAN(PI*X)
C---------------------------------------------------------------------
      IF (ABS(X) .GT. XSMALL) GO TO 100 
      IF (X .EQ. 0.0E0) GO TO 400
C---------------------------------------------------------------------
C     0 .LT. ABS(X) .LE. XSMALL.  USE 1/X AS A SUBSTITUTE
C     FOR  PI*COTAN(PI*X)
C---------------------------------------------------------------------
      AUG = -1.0E0 / X
      GO TO 150
C---------------------------------------------------------------------
C     REDUCTION OF ARGUMENT FOR COTAN
C---------------------------------------------------------------------
  100 W = - X
      SGN = PIOV4
      IF (W .GT. 0.0E0) GO TO 120
      W = - W
      SGN = -SGN
C---------------------------------------------------------------------
C     MAKE AN ERROR EXIT IF X .LE. -XMAX1
C---------------------------------------------------------------------
  120 IF (W .GE. XMAX1) GO TO 400
      NQ = INT(W)
      W = W - FLOAT(NQ)
      NQ = INT(W*4.0E0)
      W = 4.0E0 * (W - FLOAT(NQ) * .25E0)
C---------------------------------------------------------------------
C     W IS NOW RELATED TO THE FRACTIONAL PART OF  4.0 * X.
C     ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST
C     QUADRANT AND DETERMINE SIGN
C---------------------------------------------------------------------
      N = NQ / 2
      IF ((N+N) .NE. NQ) W = 1.0E0 - W
      Z = PIOV4 * W 
      M = N / 2
      IF ((M+M) .NE. N) SGN = - SGN
C---------------------------------------------------------------------
C     DETERMINE FINAL VALUE FOR  -PI*COTAN(PI*X)
C---------------------------------------------------------------------
      N = (NQ + 1) / 2
      M = N / 2
      M = M + M
      IF (M .NE. N) GO TO 140 
C---------------------------------------------------------------------
C     CHECK FOR SINGULARITY
C---------------------------------------------------------------------
      IF (Z .EQ. 0.0E0) GO TO 400
C---------------------------------------------------------------------
C     USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND
C     SIN/COS AS A SUBSTITUTE FOR TAN
C---------------------------------------------------------------------
      AUG = SGN * ((COS(Z) / SIN(Z)) * 4.0E0)
      GO TO 150
  140 AUG = SGN * ((SIN(Z) / COS(Z)) * 4.0E0)
  150 X = 1.0E0 - X 
  200 IF (X .GT. 3.0E0) GO TO 300
C---------------------------------------------------------------------
C     0.5 .LE. X .LE. 3.0
C---------------------------------------------------------------------
      DEN = X
      UPPER = P1(1) * X
C
      DO 210 I = 1, 5
         DEN = (DEN + Q1(I)) * X
         UPPER = (UPPER + P1(I+1)) * X
  210 CONTINUE
C
      DEN = (UPPER + P1(7)) / (DEN + Q1(6))
      XMX0 = DBLE(X) - DX0
      PSI = DEN * XMX0 + AUG
      RETURN
C---------------------------------------------------------------------
C     IF X .GE. XMAX1, PSI = LN(X)
C---------------------------------------------------------------------
  300 IF (X .GE. XMAX1) GO TO 350
C---------------------------------------------------------------------
C     3.0 .LT. X .LT. XMAX1
C---------------------------------------------------------------------
      W = 1.0E0 / (X * X)
      DEN = W
      UPPER = P2(1) * W
C
      DO 310 I = 1, 3
         DEN = (DEN + Q2(I)) * W
         UPPER = (UPPER + P2(I+1)) * W
  310 CONTINUE
C
      AUG = UPPER / (DEN + Q2(4)) - 0.5E0 / X + AUG
  350 PSI = AUG + ALOG(X)
      RETURN
C---------------------------------------------------------------------
C     ERROR RETURN
C---------------------------------------------------------------------
  400 PSI = 0.0E0
      RETURN
      END 
      REAL FUNCTION BETALN (A0, B0)
C-----------------------------------------------------------------------
C     EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION
C-----------------------------------------------------------------------
C     E = 0.5*LN(2*PI)
C--------------------------
      DATA E /.918938533204673/
C--------------------------
      A = AMIN1(A0,B0)
      B = AMAX1(A0,B0)
      IF (A .GE. 8.0) GO TO 60
      IF (A .GE. 1.0) GO TO 20
C-----------------------------------------------------------------------
C                   PROCEDURE WHEN A .LT. 1
C-----------------------------------------------------------------------
      IF (B .GE. 8.0) GO TO 10
         BETALN = GAMLN(A) + (GAMLN(B) - GAMLN(A + B))
         RETURN
   10 BETALN = GAMLN(A) + ALGDIV(A,B)
      RETURN
C-----------------------------------------------------------------------
C                PROCEDURE WHEN 1 .LE. A .LT. 8
C-----------------------------------------------------------------------
   20 IF (A .GT. 2.0) GO TO 30
      IF (B .GT. 2.0) GO TO 21
         BETALN = GAMLN(A) + GAMLN(B) - GSUMLN(A,B)
         RETURN
   21 W = 0.0
      IF (B .LT. 8.0) GO TO 40
         BETALN = GAMLN(A) + ALGDIV(A,B)
         RETURN
C
C                REDUCTION OF A WHEN B .LE. 1000
C
   30 IF (B .GT. 1000.0) GO TO 50
      N = A - 1.0
      W = 1.0
      DO 31 I = 1,N 
         A = A - 1.0
         H = A/B
         W = W * (H/(1.0 + H))
   31 CONTINUE
      W = ALOG(W)
      IF (B .LT. 8.0) GO TO 40
      BETALN = W + GAMLN(A) + ALGDIV(A,B)
      RETURN
C
C                 REDUCTION OF B WHEN B .LT. 8
C
   40 N = B - 1.0
      Z = 1.0
      DO 41 I = 1,N 
         B = B - 1.0
         Z = Z * (B/(A + B))
   41 CONTINUE
      BETALN = W + ALOG(Z) + (GAMLN(A) + (GAMLN(B) - GSUMLN(A,B)))
      RETURN
C
C                REDUCTION OF A WHEN B .GT. 1000
C
   50 N = A - 1.0
      W = 1.0
      DO 51 I = 1,N 
         A = A - 1.0
         W = W * (A/(1.0 + A/B))
   51 CONTINUE
      BETALN = (ALOG(W) - N*ALOG(B)) + (GAMLN(A) + ALGDIV(A,B))
      RETURN
C-----------------------------------------------------------------------
C                   PROCEDURE WHEN A .GE. 8
C-----------------------------------------------------------------------
   60 W = BCORR(A,B)
      H = A/B
      C = H/(1.0 + H)
      U = -(A - 0.5)*ALOG(C)
      V = B*ALNREL(H)
      IF (U .LE. V) GO TO 61
         BETALN = (((-0.5*ALOG(B) + E) + W) - V) - U
         RETURN
   61 BETALN = (((-0.5*ALOG(B) + E) + W) - U) - V 
      RETURN
      END 
      REAL FUNCTION GSUMLN (A, B)
C-----------------------------------------------------------------------
C          EVALUATION OF THE FUNCTION LN(GAMMA(A + B))
C          FOR 1 .LE. A .LE. 2  AND  1 .LE. B .LE. 2
C-----------------------------------------------------------------------
      X = DBLE(A) + DBLE(B) - 2.D0
      IF (X .GT. 0.25) GO TO 10
         GSUMLN = GAMLN1(1.0 + X)
         RETURN
   10 IF (X .GT. 1.25) GO TO 20
         GSUMLN = GAMLN1(X) + ALNREL(X) 
         RETURN
   20 GSUMLN = GAMLN1(X - 1.0) + ALOG(X*(1.0 + X))
      RETURN
      END 
      REAL FUNCTION BCORR (A0, B0)
C-----------------------------------------------------------------------
C
C     EVALUATION OF  DEL(A0) + DEL(B0) - DEL(A0 + B0)  WHERE
C     LN(GAMMA(A)) = (A - 0.5)*LN(A) - A + 0.5*LN(2*PI) + DEL(A).
C     IT IS ASSUMED THAT A0 .GE. 8 AND B0 .GE. 8. 
C
C-----------------------------------------------------------------------
      DATA C0/.833333333333333E-01/, C1/-.277777777760991E-02/,
     *     C2/.793650666825390E-03/, C3/-.595202931351870E-03/,
     *     C4/.837308034031215E-03/, C5/-.165322962780713E-02/
C------------------------
      A = AMIN1(A0, B0)
      B = AMAX1(A0, B0)
C
      H = A/B
      C = H/(1.0 + H)
      X = 1.0/(1.0 + H)
      X2 = X*X
C
C                SET SN = (1 - X**N)/(1 - X)
C
      S3 = 1.0 + (X + X2)
      S5 = 1.0 + (X + X2*S3)
      S7 = 1.0 + (X + X2*S5)
      S9 = 1.0 + (X + X2*S7)
      S11 = 1.0 + (X + X2*S9) 
C
C                SET W = DEL(B) - DEL(A + B)
C
      T = (1.0/B)**2
      W = ((((C5*S11*T + C4*S9)*T + C3*S7)*T + C2*S5)*T + C1*S3)*T + C0
      W = W*(C/B)
C
C                   COMPUTE  DEL(A) + W 
C
      T = (1.0/A)**2
      BCORR = (((((C5*T + C4)*T + C3)*T + C2)*T + C1)*T + C0)/A + W
      RETURN
      END 
      REAL FUNCTION ALGDIV (A, B)
C-----------------------------------------------------------------------
C
C     COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B .GE. 8
C
C                         --------
C
C     IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY
C     LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X).
C
C-----------------------------------------------------------------------
      DATA C0/.833333333333333E-01/, C1/-.277777777760991E-02/,
     *     C2/.793650666825390E-03/, C3/-.595202931351870E-03/,
     *     C4/.837308034031215E-03/, C5/-.165322962780713E-02/
C------------------------
      IF (A .LE. B) GO TO 10
         H = B/A
         C = 1.0/(1.0 + H)
         X = H/(1.0 + H)
         D = A + (B - 0.5)
         GO TO 20
   10 H = A/B
      C = H/(1.0 + H)
      X = 1.0/(1.0 + H)
      D = B + (A - 0.5)
C
C                SET SN = (1 - X**N)/(1 - X)
C
   20 X2 = X*X
      S3 = 1.0 + (X + X2)
      S5 = 1.0 + (X + X2*S3)
      S7 = 1.0 + (X + X2*S5)
      S9 = 1.0 + (X + X2*S7)
      S11 = 1.0 + (X + X2*S9) 
C
C                SET W = DEL(B) - DEL(A + B)
C
      T = (1.0/B)**2
      W = ((((C5*S11*T + C4*S9)*T + C3*S7)*T + C2*S5)*T + C1*S3)*T + C0
      W = W*(C/B)
C
C                    COMBINE THE RESULTS
C
      U = D*ALNREL(A/B)
      V = A*(ALOG(B) - 1.0)
      IF (U .LE. V) GO TO 30
         ALGDIV = (W - V) - U 
         RETURN
   30 ALGDIV = (W - U) - V
      RETURN
      END 
      REAL FUNCTION GAMLN (A) 
C-----------------------------------------------------------------------
C            EVALUATION OF LN(GAMMA(A)) FOR POSITIVE A
C-----------------------------------------------------------------------
C     WRITTEN BY ALFRED H. MORRIS
C          NAVAL SURFACE WARFARE CENTER 
C          DAHLGREN, VIRGINIA 
C--------------------------
C     D = 0.5*(LN(2*PI) - 1)
C--------------------------
      DATA D/.418938533204673/
C--------------------------
      DATA C0/.833333333333333E-01/, C1/-.277777777760991E-02/,
     *     C2/.793650666825390E-03/, C3/-.595202931351870E-03/,
     *     C4/.837308034031215E-03/, C5/-.165322962780713E-02/
C-----------------------------------------------------------------------
      IF (A .GT. 0.8) GO TO 10
         GAMLN = GAMLN1(A) - ALOG(A)
         RETURN
   10 IF (A .GT. 2.25) GO TO 20
         T = (A - 0.5) - 0.5
         GAMLN = GAMLN1(T)
         RETURN
C
   20 IF (A .GE. 10.0) GO TO 30
      N = A - 1.25
      T = A
      W = 1.0
      DO 21 I = 1,N 
         T = T - 1.0
   21    W = T*W
      GAMLN = GAMLN1(T - 1.0) + ALOG(W) 
      RETURN
C
   30 T = (1.0/A)**2
      W = (((((C5*T + C4)*T + C3)*T + C2)*T + C1)*T + C0)/A 
      GAMLN = (D + W) + (A - 0.5)*(ALOG(A) - 1.0) 
      END 
      INTEGER FUNCTION IPMPAR (I)
C-----------------------------------------------------------------------
C
C     IPMPAR PROVIDES THE INTEGER MACHINE CONSTANTS FOR THE COMPUTER
C     THAT IS USED. IT IS ASSUMED THAT THE ARGUMENT I IS AN INTEGER
C     HAVING ONE OF THE VALUES 1-10. IPMPAR(I) HAS THE VALUE ...
C
C  INTEGERS.
C
C     ASSUME INTEGERS ARE REPRESENTED IN THE N-DIGIT, BASE-A FORM
C
C               SIGN ( X(N-1)*A**(N-1) + ... + X(1)*A + X(0) )
C
C               WHERE 0 .LE. X(I) .LT. A FOR I=0,...,N-1.
C
C     IPMPAR(1) = A, THE BASE.
C
C     IPMPAR(2) = N, THE NUMBER OF BASE-A DIGITS. 
C
C     IPMPAR(3) = A**N - 1, THE LARGEST MAGNITUDE.
C
C  FLOATING-POINT NUMBERS.
C
C     IT IS ASSUMED THAT THE SINGLE AND DOUBLE PRECISION FLOATING
C     POINT ARITHMETICS HAVE THE SAME BASE, SAY B, AND THAT THE
C     NONZERO NUMBERS ARE REPRESENTED IN THE FORM 
C
C               SIGN (B**E) * (X(1)/B + ... + X(M)/B**M)
C
C               WHERE X(I) = 0,1,...,B-1 FOR I=1,...,M,
C               X(1) .GE. 1, AND EMIN .LE. E .LE. EMAX.
C
C     IPMPAR(4) = B, THE BASE.
C
C  SINGLE-PRECISION 
C
C     IPMPAR(5) = M, THE NUMBER OF BASE-B DIGITS. 
C
C     IPMPAR(6) = EMIN, THE SMALLEST EXPONENT E.
C
C     IPMPAR(7) = EMAX, THE LARGEST EXPONENT E.
C
C  DOUBLE-PRECISION 
C
C     IPMPAR(8) = M, THE NUMBER OF BASE-B DIGITS. 
C
C     IPMPAR(9) = EMIN, THE SMALLEST EXPONENT E.
C
C     IPMPAR(10) = EMAX, THE LARGEST EXPONENT E.
C
C-----------------------------------------------------------------------
C
C     TO DEFINE THIS FUNCTION FOR THE COMPUTER BEING USED, ACTIVATE
C     THE DATA STATMENTS FOR THE COMPUTER BY REMOVING THE C FROM
C     COLUMN 1. (ALL THE OTHER DATA STATEMENTS SHOULD HAVE C IN
C     COLUMN 1.)
C
C-----------------------------------------------------------------------
C
C     IPMPAR IS AN ADAPTATION OF THE FUNCTION I1MACH, WRITTEN BY
C     P.A. FOX, A.D. HALL, AND N.L. SCHRYER (BELL LABORATORIES).
C     IPMPAR WAS FORMED BY A.H. MORRIS (NSWC). THE CONSTANTS ARE
C     FROM BELL LABORATORIES, NSWC, AND OTHER SOURCES.
C
C-----------------------------------------------------------------------
      INTEGER IMACH(10)
C
C     MACHINE CONSTANTS FOR AMDAHL MACHINES.
C
C     DATA IMACH( 1) /   2 /
C     DATA IMACH( 2) /  31 /
C     DATA IMACH( 3) / 2147483647 /
C     DATA IMACH( 4) /  16 /
C     DATA IMACH( 5) /   6 /
C     DATA IMACH( 6) / -64 /
C     DATA IMACH( 7) /  63 /
C     DATA IMACH( 8) /  14 /
C     DATA IMACH( 9) / -64 /
C     DATA IMACH(10) /  63 /
C
C     MACHINE CONSTANTS FOR THE AT&T 3B SERIES, AT&T
C     PC 7300, AND AT&T 6300. 
C
C     DATA IMACH( 1) /     2 /
C     DATA IMACH( 2) /    31 /
C     DATA IMACH( 3) / 2147483647 /
C     DATA IMACH( 4) /     2 /
C     DATA IMACH( 5) /    24 /
C     DATA IMACH( 6) /  -125 /
C     DATA IMACH( 7) /   128 /
C     DATA IMACH( 8) /    53 /
C     DATA IMACH( 9) / -1021 /
C     DATA IMACH(10) /  1024 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C     DATA IMACH( 1) /    2 / 
C     DATA IMACH( 2) /   33 / 
C     DATA IMACH( 3) / 8589934591 /
C     DATA IMACH( 4) /    2 / 
C     DATA IMACH( 5) /   24 / 
C     DATA IMACH( 6) / -256 / 
C     DATA IMACH( 7) /  255 / 
C     DATA IMACH( 8) /   60 / 
C     DATA IMACH( 9) / -256 / 
C     DATA IMACH(10) /  255 / 
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
C
C     DATA IMACH( 1) /    2 / 
C     DATA IMACH( 2) /   39 / 
C     DATA IMACH( 3) / 549755813887 /
C     DATA IMACH( 4) /    8 / 
C     DATA IMACH( 5) /   13 / 
C     DATA IMACH( 6) /  -50 / 
C     DATA IMACH( 7) /   76 / 
C     DATA IMACH( 8) /   26 / 
C     DATA IMACH( 9) /  -50 / 
C     DATA IMACH(10) /   76 / 
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
C
C     DATA IMACH( 1) /      2 /
C     DATA IMACH( 2) /     39 /
C     DATA IMACH( 3) / 549755813887 /
C     DATA IMACH( 4) /      8 /
C     DATA IMACH( 5) /     13 /
C     DATA IMACH( 6) /    -50 /
C     DATA IMACH( 7) /     76 /
C     DATA IMACH( 8) /     26 /
C     DATA IMACH( 9) / -32754 /
C     DATA IMACH(10) /  32780 /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES
C     60 BIT ARITHMETIC, AND THE CDC CYBER 995 64 BIT
C     ARITHMETIC (NOS OPERATING SYSTEM).
C
      DATA IMACH( 1) /    2 / 
      DATA IMACH( 2) /   48 / 
      DATA IMACH( 3) / 281474976710655 /
      DATA IMACH( 4) /    2 / 
      DATA IMACH( 5) /   48 / 
      DATA IMACH( 6) / -974 / 
      DATA IMACH( 7) / 1070 / 
      DATA IMACH( 8) /   95 / 
      DATA IMACH( 9) / -926 / 
      DATA IMACH(10) / 1070 / 
C
C     MACHINE CONSTANTS FOR THE CDC CYBER 995 64 BIT
C     ARITHMETIC (NOS/VE OPERATING SYSTEM).
C
C     DATA IMACH( 1) /     2 /
C     DATA IMACH( 2) /    63 /
C     DATA IMACH( 3) / 9223372036854775807 /
C     DATA IMACH( 4) /     2 /
C     DATA IMACH( 5) /    48 /
C     DATA IMACH( 6) / -4096 /
C     DATA IMACH( 7) /  4095 /
C     DATA IMACH( 8) /    96 /
C     DATA IMACH( 9) / -4096 /
C     DATA IMACH(10) /  4095 /
C
C     MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3.
C
C     DATA IMACH( 1) /     2 /
C     DATA IMACH( 2) /    63 /
C     DATA IMACH( 3) / 9223372036854775807 /
C     DATA IMACH( 4) /     2 /
C     DATA IMACH( 5) /    47 /
C     DATA IMACH( 6) / -8189 /
C     DATA IMACH( 7) /  8190 /
C     DATA IMACH( 8) /    94 /
C     DATA IMACH( 9) / -8099 /
C     DATA IMACH(10) /  8190 /
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200. 
C
C     DATA IMACH( 1) /    2 / 
C     DATA IMACH( 2) /   15 / 
C     DATA IMACH( 3) / 32767 /
C     DATA IMACH( 4) /   16 / 
C     DATA IMACH( 5) /    6 / 
C     DATA IMACH( 6) /  -64 / 
C     DATA IMACH( 7) /   63 / 
C     DATA IMACH( 8) /   14 / 
C     DATA IMACH( 9) /  -64 / 
C     DATA IMACH(10) /   63 / 
C
C     MACHINE CONSTANTS FOR THE HARRIS 220.
C
C     DATA IMACH( 1) /    2 / 
C     DATA IMACH( 2) /   23 / 
C     DATA IMACH( 3) / 8388607 /
C     DATA IMACH( 4) /    2 / 
C     DATA IMACH( 5) /   23 / 
C     DATA IMACH( 6) / -127 / 
C     DATA IMACH( 7) /  127 / 
C     DATA IMACH( 8) /   38 / 
C     DATA IMACH( 9) / -127 / 
C     DATA IMACH(10) /  127 / 
C
C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000
C     AND DPS 8/70 SERIES.
C
C     DATA IMACH( 1) /    2 / 
C     DATA IMACH( 2) /   35 / 
C     DATA IMACH( 3) / 34359738367 /
C     DATA IMACH( 4) /    2 / 
C     DATA IMACH( 5) /   27 / 
C     DATA IMACH( 6) / -127 / 
C     DATA IMACH( 7) /  127 / 
C     DATA IMACH( 8) /   63 / 
C     DATA IMACH( 9) / -127 / 
C     DATA IMACH(10) /  127 / 
C
C     MACHINE CONSTANTS FOR THE HP 2100 
C     3 WORD DOUBLE PRECISION OPTION WITH FTN4
C
C     DATA IMACH( 1) /    2 / 
C     DATA IMACH( 2) /   15 / 
C     DATA IMACH( 3) / 32767 /
C     DATA IMACH( 4) /    2 / 
C     DATA IMACH( 5) /   23 / 
C     DATA IMACH( 6) / -128 / 
C     DATA IMACH( 7) /  127 / 
C     DATA IMACH( 8) /   39 / 
C     DATA IMACH( 9) / -128 / 
C     DATA IMACH(10) /  127 / 
C
C     MACHINE CONSTANTS FOR THE HP 2100 
C     4 WORD DOUBLE PRECISION OPTION WITH FTN4
C
C     DATA IMACH( 1) /    2 / 
C     DATA IMACH( 2) /   15 / 
C     DATA IMACH( 3) / 32767 /
C     DATA IMACH( 4) /    2 / 
C     DATA IMACH( 5) /   23 / 
C     DATA IMACH( 6) / -128 / 
C     DATA IMACH( 7) /  127 / 
C     DATA IMACH( 8) /   55 / 
C     DATA IMACH( 9) / -128 / 
C     DATA IMACH(10) /  127 / 
C
C     MACHINE CONSTANTS FOR THE HP 9000.
C
C     DATA IMACH( 1) /     2 /
C     DATA IMACH( 2) /    31 /
C     DATA IMACH( 3) / 2147483647 /
C     DATA IMACH( 4) /     2 /
C     DATA IMACH( 5) /    24 /
C     DATA IMACH( 6) /  -126 /
C     DATA IMACH( 7) /   128 /
C     DATA IMACH( 8) /    53 /
C     DATA IMACH( 9) / -1021 /
C     DATA IMACH(10) /  1024 /
C
C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C     THE ICL 2900, THE ITEL AS/6, THE XEROX SIGMA
C     5/7/9 AND THE SEL SYSTEMS 85/86.
C
C     DATA IMACH( 1) /    2 / 
C     DATA IMACH( 2) /   31 / 
C     DATA IMACH( 3) / 2147483647 /
C     DATA IMACH( 4) /   16 / 
C     DATA IMACH( 5) /    6 / 
C     DATA IMACH( 6) /  -64 / 
C     DATA IMACH( 7) /   63 / 
C     DATA IMACH( 8) /   14 / 
C     DATA IMACH( 9) /  -64 / 
C     DATA IMACH(10) /   63 / 
C
C     MACHINE CONSTANTS FOR THE IBM PC. 
C
C     DATA IMACH( 1) /     2 /
C     DATA IMACH( 2) /    31 /
C     DATA IMACH( 3) / 2147483647 /
C     DATA IMACH( 4) /     2 /
C     DATA IMACH( 5) /    24 /
C     DATA IMACH( 6) /  -125 /
C     DATA IMACH( 7) /   128 /
C     DATA IMACH( 8) /    53 /
C     DATA IMACH( 9) / -1021 /
C     DATA IMACH(10) /  1024 /
C
C     MACHINE CONSTANTS FOR THE MACINTOSH II - ABSOFT
C     MACFORTRAN II.
C
C     DATA IMACH( 1) /     2 /
C     DATA IMACH( 2) /    31 /
C     DATA IMACH( 3) / 2147483647 /
C     DATA IMACH( 4) /     2 /
C     DATA IMACH( 5) /    24 /
C     DATA IMACH( 6) /  -125 /
C     DATA IMACH( 7) /   128 /
C     DATA IMACH( 8) /    53 /
C     DATA IMACH( 9) / -1021 /
C     DATA IMACH(10) /  1024 /
C
C     MACHINE CONSTANTS FOR THE MICROVAX - VMS FORTRAN.
C
C     DATA IMACH( 1) /    2 / 
C     DATA IMACH( 2) /   31 / 
C     DATA IMACH( 3) / 2147483647 /
C     DATA IMACH( 4) /    2 / 
C     DATA IMACH( 5) /   24 / 
C     DATA IMACH( 6) / -127 / 
C     DATA IMACH( 7) /  127 / 
C     DATA IMACH( 8) /   56 / 
C     DATA IMACH( 9) / -127 / 
C     DATA IMACH(10) /  127 / 
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
C
C     DATA IMACH( 1) /    2 / 
C     DATA IMACH( 2) /   35 / 
C     DATA IMACH( 3) / 34359738367 /
C     DATA IMACH( 4) /    2 / 
C     DATA IMACH( 5) /   27 / 
C     DATA IMACH( 6) / -128 / 
C     DATA IMACH( 7) /  127 / 
C     DATA IMACH( 8) /   54 / 
C     DATA IMACH( 9) / -101 / 
C     DATA IMACH(10) /  127 / 
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
C
C     DATA IMACH( 1) /    2 / 
C     DATA IMACH( 2) /   35 / 
C     DATA IMACH( 3) / 34359738367 /
C     DATA IMACH( 4) /    2 / 
C     DATA IMACH( 5) /   27 / 
C     DATA IMACH( 6) / -128 / 
C     DATA IMACH( 7) /  127 / 
C     DATA IMACH( 8) /   62 / 
C     DATA IMACH( 9) / -128 / 
C     DATA IMACH(10) /  127 / 
C
C     MACHINE CONSTANTS FOR THE PDP-11 FORTRAN SUPPORTING
C     32-BIT INTEGER ARITHMETIC.
C
C     DATA IMACH( 1) /    2 / 
C     DATA IMACH( 2) /   31 / 
C     DATA IMACH( 3) / 2147483647 /
C     DATA IMACH( 4) /    2 / 
C     DATA IMACH( 5) /   24 / 
C     DATA IMACH( 6) / -127 / 
C     DATA IMACH( 7) /  127 / 
C     DATA IMACH( 8) /   56 / 
C     DATA IMACH( 9) / -127 / 
C     DATA IMACH(10) /  127 / 
C
C     MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000.
C
C     DATA IMACH( 1) /     2 /
C     DATA IMACH( 2) /    31 /
C     DATA IMACH( 3) / 2147483647 /
C     DATA IMACH( 4) /     2 /
C     DATA IMACH( 5) /    24 /
C     DATA IMACH( 6) /  -125 /
C     DATA IMACH( 7) /   128 /
C     DATA IMACH( 8) /    53 /
C     DATA IMACH( 9) / -1021 /
C     DATA IMACH(10) /  1024 /
C
C     MACHINE CONSTANTS FOR THE SILICON GRAPHICS IRIS-4D
C     SERIES (MIPS R3000 PROCESSOR).
C
C     DATA IMACH( 1) /     2 /
C     DATA IMACH( 2) /    31 /
C     DATA IMACH( 3) / 2147483647 /
C     DATA IMACH( 4) /     2 /
C     DATA IMACH( 5) /    24 /
C     DATA IMACH( 6) /  -125 /
C     DATA IMACH( 7) /   128 /
C     DATA IMACH( 8) /    53 /
C     DATA IMACH( 9) / -1021 /
C     DATA IMACH(10) /  1024 /
C
C     MACHINE CONSTANTS FOR THE SUN 3.
C
C     DATA IMACH( 1) /     2 /
C     DATA IMACH( 2) /    31 /
C     DATA IMACH( 3) / 2147483647 /
C     DATA IMACH( 4) /     2 /
C     DATA IMACH( 5) /    24 /
C     DATA IMACH( 6) /  -125 /
C     DATA IMACH( 7) /   128 /
C     DATA IMACH( 8) /    53 /
C     DATA IMACH( 9) / -1021 /
C     DATA IMACH(10) /  1024 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
C
C     DATA IMACH( 1) /    2 / 
C     DATA IMACH( 2) /   35 / 
C     DATA IMACH( 3) / 34359738367 /
C     DATA IMACH( 4) /    2 / 
C     DATA IMACH( 5) /   27 / 
C     DATA IMACH( 6) / -128 / 
C     DATA IMACH( 7) /  127 / 
C     DATA IMACH( 8) /   60 / 
C     DATA IMACH( 9) /-1024 / 
C     DATA IMACH(10) / 1023 / 
C
C     MACHINE CONSTANTS FOR THE VAX 11/780.
C
C     DATA IMACH( 1) /    2 / 
C     DATA IMACH( 2) /   31 / 
C     DATA IMACH( 3) / 2147483647 /
C     DATA IMACH( 4) /    2 / 
C     DATA IMACH( 5) /   24 / 
C     DATA IMACH( 6) / -127 / 
C     DATA IMACH( 7) /  127 / 
C     DATA IMACH( 8) /   56 / 
C     DATA IMACH( 9) / -127 / 
C     DATA IMACH(10) /  127 / 
C
      IPMPAR = IMACH(I)
      RETURN
      END 

		
.

NEW PAGES:

[ODDNUGGET]

[GOPHER]