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