C**** SUBROUTINE DMRRRR(NRA,NCA,A,LDA,NRB,NCB,B,LDB,NRC,NCC,C,LDC) C=============================================================================== C C SIMPLE MATRIX MULTIPLICATION C C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DIMENSION A(LDA,NCA),B(LDB,NCB),C(LDC,NCC) IF( 1 ( NCA .NE. NRB ) .OR. 2 ( NRA .NE. NRC ) .OR. 3 ( NCB .NE. NCC ) ) THEN WRITE(*,*) ' NON-CONFORMING MATRICES (DMRRRR)' STOP ENDIF DO 10 I=1,NRC DO 10 J=1,NCC C(I,J) = 0.D0 DO 10 K=1,NCA C(I,J) = C(I,J) + A(I,K)*B(K,J) 10 CONTINUE RETURN END C**** SUBROUTINE DMXTYF(NRA,NCA,A,LDA,NRB,NCB,B,LDB,NRC,NCC,C,LDC) C=============================================================================== C C MATRIX MULTIPLICATION TRANS(A)*B C C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DIMENSION A(LDA,NCA),B(LDB,NCB),C(LDC,NCC) IF( 1 ( NRA .NE. NRB ) .OR. 2 ( NCA .NE. NRC ) .OR. 3 ( NCB .NE. NCC ) ) THEN WRITE(*,*) ' NON-CONFORMING MATRICES (DMXTYF)' STOP ENDIF DO 10 I=1,NRC DO 10 J=1,NCC C(I,J) = 0.D0 DO 10 K=1,NRA C(I,J) = C(I,J) + A(K,I)*B(K,J) 10 CONTINUE RETURN END C**** SUBROUTINE DMXYTF(NRA,NCA,A,LDA,NRB,NCB,B,LDB,NRC,NCC,C,LDC) C=============================================================================== C C MATRIX MULTIPLICATION A*TRANS(B) C C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DIMENSION A(LDA,NCA),B(LDB,NCB),C(LDC,NCC) IF( 1 ( NRA .NE. NRB ) .OR. 2 ( NRA .NE. NRC ) .OR. 3 ( NRB .NE. NCC ) ) THEN WRITE(*,*) ' NON-CONFORMING MATRICES (DMXYTF)' STOP ENDIF DO 10 I=1,NRC DO 10 J=1,NCC C(I,J) = 0.D0 DO 10 K=1,NRA C(I,J) = C(I,J) + A(I,K)*B(J,K) 10 CONTINUE RETURN END C**** SUBROUTINE DMSUM(NRA,NCA,A,LDA,NRB,NCB,B,LDB,NRC,NCC,C,LDC) C=============================================================================== C C MATRIX ADDITION C = A + B C C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DIMENSION A(LDA,NCA),B(LDB,NCB),C(LDC,NCC) IF( 1 ( NRA .NE. NRB ) .OR. 2 ( NCA .NE. NCB ) .OR. 3 ( NRB .NE. NRC ) .OR. 3 ( NCB .NE. NCC ) ) THEN WRITE(*,*) ' NON-CONFORMING MATRICES (DMSUM)' STOP ENDIF DO 10 I=1,NRC DO 10 J=1,NCC C(I,J) = A(I,J) + B(I,J) 10 CONTINUE RETURN END C**** SUBROUTINE DMMULT(VAL,NRA,NCA,A,LDA,NRC,NCC,C,LDC) C=============================================================================== C C MATRIX ADDITION C = VAL * A C C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DIMENSION A(LDA,NCA),C(LDC,NCC) IF( 1 ( NRA .NE. NRC ) .OR. 2 ( NCA .NE. NCC ) ) THEN WRITE(*,*) ' NON-CONFORMING MATRICES (DMSUM)' STOP ENDIF DO 10 I=1,NRC DO 10 J=1,NCC C(I,J) = VAL*A(I,J) 10 CONTINUE RETURN END C**** SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) C=============================================================================== C C LET Y(I) = DA * X(I) + Y(I) C C C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DIMENSION DX(*),DY(*) IF(INCX .GE. 0) THEN INDX = 1 ELSE INDX = (1-N)*INCX+1 ENDIF IF(INCY .GE. 0) THEN INDY = 1 ELSE INDY = (1-N)*INCY+1 ENDIF DO 10 I=1,N DY(INDY) = DY(INDY) + DA*DX(INDX) INDX = INDX + INCX INDY = INDY + INCY 10 CONTINUE RETURN END C**** FUNCTION DLNGAM(XX) C=============================================================================== C C LOGARITHM OF THE GAMMA FUNCTION C C N.B. VALID ONLY FOR XX>0 C THIS IS DIFFERENT FROM IMSL ROUTINE C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DIMENSION COF(6) DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0, * -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/ DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/ X=XX-ONE TMP=X+FPF TMP=(X+HALF)*LOG(TMP)-TMP SER=ONE DO 11 J=1,6 X=X+ONE SER=SER+COF(J)/X 11 CONTINUE DLNGAM=TMP+LOG(STP*SER) RETURN END C**** SUBROUTINE DLINRG(N,A,LDA,AINV,LDAINV) C=============================================================================== C C SIMPLE MATRIX INVERSION C C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DIMENSION A(LDA,LDA),AINV(LDAINV,N) DO 10 I=1,N DO 10 J=1,N AINV(I,J) = A(I,J) 10 CONTINUE CALL GAUSSJ(AINV,N,LDAINV,B,0,1) RETURN END C**** SUBROUTINE GAUSSJ(A,N,NP,B,M,MP) C C PERFORMS GAUSS-JORDAN ELIMINATION C C ON RETURN A CONTAINS INVERSE OF A, AND, IF B CONTAINS RHS C OF EQUATION C AX = B C THEN ON RETURN B WILL CONTAIN SOLUTION X C C TIM COHN 8/22/97 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE PARAMETER (NMAX=50) DIMENSION A(NP,NP),B(NP,MP),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX) DO 11 J=1,N IPIV(J)=0 11 CONTINUE DO 22 I=1,N BIG=0. DO 13 J=1,N IF(IPIV(J).NE.1)THEN DO 12 K=1,N IF (IPIV(K).EQ.0) THEN IF (ABS(A(J,K)).GE.BIG)THEN BIG=ABS(A(J,K)) IROW=J ICOL=K ENDIF ELSE IF (IPIV(K).GT.1) THEN PAUSE 'SINGULAR MATRIX' ENDIF 12 CONTINUE ENDIF 13 CONTINUE IPIV(ICOL)=IPIV(ICOL)+1 IF (IROW.NE.ICOL) THEN DO 14 L=1,N DUM=A(IROW,L) A(IROW,L)=A(ICOL,L) A(ICOL,L)=DUM 14 CONTINUE DO 15 L=1,M DUM=B(IROW,L) B(IROW,L)=B(ICOL,L) B(ICOL,L)=DUM 15 CONTINUE ENDIF INDXR(I)=IROW INDXC(I)=ICOL IF (A(ICOL,ICOL).EQ.0.) PAUSE 'SINGULAR MATRIX.' PIVINV=1./A(ICOL,ICOL) A(ICOL,ICOL)=1. DO 16 L=1,N A(ICOL,L)=A(ICOL,L)*PIVINV 16 CONTINUE DO 17 L=1,M B(ICOL,L)=B(ICOL,L)*PIVINV 17 CONTINUE DO 21 LL=1,N IF(LL.NE.ICOL)THEN DUM=A(LL,ICOL) A(LL,ICOL)=0. DO 18 L=1,N A(LL,L)=A(LL,L)-A(ICOL,L)*DUM 18 CONTINUE DO 19 L=1,M B(LL,L)=B(LL,L)-B(ICOL,L)*DUM 19 CONTINUE ENDIF 21 CONTINUE 22 CONTINUE DO 24 L=N,1,-1 IF(INDXR(L).NE.INDXC(L))THEN DO 23 K=1,N DUM=A(K,INDXR(L)) A(K,INDXR(L))=A(K,INDXC(L)) A(K,INDXC(L))=DUM 23 CONTINUE ENDIF 24 CONTINUE RETURN END C**** SUBROUTINE ISET(N,IVALUE,IVAR,ISPACE) C=============================================================================== C C SIMPLE MATRIX INVERSION C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE INTEGER IVAR(*) DO 10 I=1,N,ISPACE IVAR(I) = IVALUE 10 CONTINUE RETURN END C**** SUBROUTINE DSET(N,FVALUE,FVAR,ISPACE) C=============================================================================== C C SIMPLE MATRIX INVERSION C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DIMENSION FVAR(*) DO 10 I=1,N,ISPACE FVAR(I) = FVALUE 10 CONTINUE RETURN END C**** DOUBLE PRECISION FUNCTION DNORDF(X) C=============================================================================== C C NORMAL DISTRIBUTION CDF C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE X2 = X * 0.7071068 DNORDF = 1.D0 - ERFCC(X2)/2.D0 RETURN END C**** DOUBLE PRECISION FUNCTION ERFCC(X) C=============================================================================== C C ERROR FUNCTION FOUND IN NUMERICAL RECIPES C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE Z=ABS(X) T=1./(1.+0.5*Z) IF( Z .LT. 15.D0 ) THEN ERFCC=T*EXP(-Z*Z-1.26551223+T*(1.00002368+T*(.37409196+ * T*(.09678418+T*(-.18628806+T*(.27886807+T*(-1.13520398+ * T*(1.48851587+T*(-.82215223+T*.17087277))))))))) ELSE ERFCC = 0.D0 ENDIF IF (X.LT.0.) ERFCC=2.-ERFCC RETURN END C**** DOUBLE PRECISION FUNCTION DNORIN(P) C=============================================================================== C C NORMAL DISTRIBUTION CDF C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE REAL*4 P2,PPF P2 = P CALL NORPPF(P2,PPF) DNORIN = PPF RETURN END C**** * ====================================================================== * NIST GUIDE TO AVAILABLE MATH SOFTWARE. * FULLSOURCE FOR MODULE NORPPF FROM PACKAGE DATAPAC. * RETRIEVED FROM CAMSUN ON FRI AUG 29 09:02:15 1997. * ====================================================================== * NORPPF SUBROUTINE NORPPF(P,PPF) C C PURPOSE--THIS SUBROUTINE COMPUTES THE PERCENT POINT C FUNCTION VALUE FOR THE NORMAL (GAUSSIAN) C DISTRIBUTION WITH MEAN = 0 AND STANDARD DEVIATION = 1. C THIS DISTRIBUTION IS DEFINED FOR ALL X AND HAS C THE PROBABILITY DENSITY FUNCTION C F(X) = (1/SQRT(2*PI))*EXP(-X*X/2). C NOTE THAT THE PERCENT POINT FUNCTION OF A DISTRIBUTION C IS IDENTICALLY THE SAME AS THE INVERSE CUMULATIVE C DISTRIBUTION FUNCTION OF THE DISTRIBUTION. C INPUT ARGUMENTS--P = THE SINGLE PRECISION VALUE C (BETWEEN 0.0 AND 1.0) C AT WHICH THE PERCENT POINT C FUNCTION IS TO BE EVALUATED. C OUTPUT ARGUMENTS--PPF = THE SINGLE PRECISION PERCENT C POINT FUNCTION VALUE. C OUTPUT--THE SINGLE PRECISION PERCENT POINT C FUNCTION VALUE PPF. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--P SHOULD BE BETWEEN 0.0 AND 1.0, EXCLUSIVELY. C OTHER DATAPAC SUBROUTINES NEEDED--NONE. C FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, ALOG. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN. C REFERENCES--ODEH AND EVANS, THE PERCENTAGE POINTS C OF THE NORMAL DISTRIBUTION, ALGORTIHM 70, C APPLIED STATISTICS, 1974, PAGES 96-97. C --EVANS, ALGORITHMS FOR MINIMAL DEGREE C POLYNOMIAL AND RATIONAL APPROXIMATION, C M. SC. THESIS, 1972, UNIVERSITY C OF VICTORIA, B. C., CANADA. C --HASTINGS, APPROXIMATIONS FOR DIGITAL C COMPUTERS, 1955, PAGES 113, 191, 192. C --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS C SERIES 55, 1964, PAGE 933, FORMULA 26.2.23. C --FILLIBEN, SIMPLE AND ROBUST LINEAR ESTIMATION C OF THE LOCATION PARAMETER OF A SYMMETRIC C DISTRIBUTION (UNPUBLISHED PH.D. DISSERTATION, C PRINCETON UNIVERSITY), 1969, PAGES 21-44, 229-231. C --FILLIBEN, 'THE PERCENT POINT FUNCTION', C (UNPUBLISHED MANUSCRIPT), 1970, PAGES 28-31. C --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE C DISTRIBUTIONS--1, 1970, PAGES 40-111. C --THE KELLEY STATISTICAL TABLES, 1948. C --OWEN, HANDBOOK OF STATISTICAL TABLES, C 1962, PAGES 3-16. C --PEARSON AND HARTLEY, BIOMETRIKA TABLES C FOR STATISTICIANS, VOLUME 1, 1954, C PAGES 104-113. C COMMENTS--THE CODING AS PRESENTED BELOW C IS ESSENTIALLY IDENTICAL TO THAT C PRESENTED BY ODEH AND EVANS C AS ALGORTIHM 70 OF APPLIED STATISTICS. C THE PRESENT AUTHOR HAS MODIFIED THE C ORIGINAL ODEH AND EVANS CODE WITH ONLY C MINOR STYLISTIC CHANGES. C --AS POINTED OUT BY ODEH AND EVANS C IN APPLIED STATISTICS, C THEIR ALGORITHM REPRESENTES A C SUBSTANTIAL IMPROVEMENT OVER THE C PREVIOUSLY EMPLOYED C HASTINGS APPROXIMATION FOR THE C NORMAL PERCENT POINT FUNCTION-- C THE ACCURACY OF APPROXIMATION C BEING IMPROVED FROM 4.5*(10**-4) C TO 1.5*(10**-8). C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING LABORATORY (205.03) C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE: 301-921-2315 C ORIGINAL VERSION--JUNE 1972. C UPDATED --SEPTEMBER 1975. C UPDATED --NOVEMBER 1975. C UPDATED --OCTOBER 1976. C C--------------------------------------------------------------------- C DATA P0,P1,P2,P3,P4 1/-.322232431088,-1.0, 1 -.342242088547,-.204231210245E-1, 1 -.453642210148E-4/ DATA Q0,Q1,Q2,Q3,Q4 1/.993484626060E-1,.588581570495, 1 .531103462366,.103537752850, 1 .38560700634E-2/ C IPR=6 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(P.LE.0.0.OR.P.GE.1.0)GOTO50 GOTO90 50 WRITE(IPR,1) WRITE(IPR,46)P RETURN 90 CONTINUE 1 FORMAT(1H ,115H***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO THE 1 NORPPF SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****) 46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****) C C-----START POINT----------------------------------------------------- C IF(P.NE.0.5)GOTO150 PPF=0.0 RETURN C 150 R=P IF(P.GT.0.5)R=1.0-R T=SQRT(-2.0*ALOG(R)) ANUM=((((T*P4+P3)*T+P2)*T+P1)*T+P0) ADEN=((((T*Q4+Q3)*T+Q2)*T+Q1)*T+Q0) PPF=T+(ANUM/ADEN) IF(P.LT.0.5)PPF=-PPF RETURN C END C**** DOUBLE PRECISION FUNCTION DGAMDF(X,A) C=============================================================================== C C GAMMA DISTRIBUTION CDF C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE X2 = MAX(0.D0,X) IF(A .LE. 0.D0) THEN WRITE(*,*) 'DGAMDF: A = ',A DGAMDF = 1.D0 RETURN ENDIF IF (A .LE. 2000.) THEN ARG1 = GAMMP(A,X2) ENDIF IF (A .GE. 1500.D0) THEN C C THIS IS A WILSON-HILFERTY SOLUTION C Z = -(-1. + 9.*A - 9.*A**(2./3.)*X2**(1./3.))/(3.*SQRT(A)) ARG2 = DNORDF(Z) ENDIF W = MAX(0.D0,MIN((A-1500.)/500.,1.D0)) DGAMDF = (1.D0-W)*ARG1 + W*ARG2 RETURN END C**** DOUBLE PRECISION FUNCTION DCHIDF(X,D) C=============================================================================== C C CHI-SQUARE DISTRIBUTION CDF C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DCHIDF = DGAMDF(X/2.D0,D/2.D0) RETURN END C**** DOUBLE PRECISION FUNCTION DCHIIN(P,NU) C=============================================================================== C C INVERSE CHI-SQUARED DISTRIBUTION CDF C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DOUBLE PRECISION P,NU,X REAL*4 GAMMA,P2 IF( (P .GT. 0.0) .AND. (P .LT. 1.D0) ) THEN P2 = P GAMMA = NU/2.0 CALL MGAMINV(NU/2.D0,P,X1,IER) X = 2.D0*X1 ELSE IF (P .LE. 0.D0) THEN X = 0.0 ELSE X = 1.D99 ENDIF DCHIIN = X RETURN END C**** DOUBLE PRECISION FUNCTION DRNGAM(N,ALPHA,X) C=============================================================================== C C GENERATES RANDOM 1P GAMMA VARIATES C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DOUBLE PRECISION X(*) INTEGER N,ISEED REAL*4 GAMMA,X2 COMMON /ZZZ889/ISEED GAMMA = ALPHA DO 10 I=1,N CALL GAMRAN(1,GAMMA,ISEED,X2) X(I) = X2 10 CONTINUE RETURN END C C**** NORRAN(N,ISEED,X) C**** DOUBLE PRECISION FUNCTION DRNORMS() C=============================================================================== C C GENERATES STANDARD NORMAL VARIATE C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE INTEGER ISEED REAL*4 X2 COMMON /ZZZ889/ISEED CALL NORRAN(1,ISEED,X2) DRNORMS = X2 RETURN END C C**** SUBROUTINE RNSET(ISEED2) C=============================================================================== C C INITIALIZED RANDOM NUMBER GENERATOR C C=============================================================================== INTEGER ISEED, ISEED2 COMMON /ZZZ889/ISEED ISEED = ISEED2 RETURN END C**** * ====================================================================== * NIST GUIDE TO AVAILABLE MATH SOFTWARE. * FULLSOURCE FOR MODULE GAMCDF FROM PACKAGE DATAPAC. * RETRIEVED FROM CAMSUN ON TUE AUG 26 16:39:05 1997. * ====================================================================== * GAMCDF SUBROUTINE GAMCDF(X,GAMMA,CDF) C C PURPOSE--THIS SUBROUTINE COMPUTES THE CUMULATIVE DISTRIBUTION C FUNCTION VALUE FOR THE GAMMA C DISTRIBUTION WITH SINGLE PRECISION C TAIL LENGTH PARAMETER = GAMMA. C THE GAMMA DISTRIBUTION USED C HEREIN HAS MEAN = GAMMA C AND STANDARD DEVIATION = SQRT(GAMMA). C THIS DISTRIBUTION IS DEFINED FOR ALL POSITIVE X, C AND HAS THE PROBABILITY DENSITY FUNCTION C F(X) = (1/CONSTANT) * (X**(GAMMA-1)) * EXP(-X) C WHERE THE CONSTANT = THE GAMMA FUNCTION EVALUATED C AT THE VALUE GAMMA. C INPUT ARGUMENTS--X = THE SINGLE PRECISION VALUE C AT WHICH THE CUMULATIVE DISTRIBUTION C FUNCTION IS TO BE EVALUATED. C X SHOULD BE POSITIVE. C --GAMMA = THE SINGLE PRECISION VALUE C OF THE TAIL LENGTH PARAMETER. C GAMMA SHOULD BE POSITIVE. C OUTPUT ARGUMENTS--CDF = THE SINGLE PRECISION CUMULATIVE C DISTRIBUTION FUNCTION VALUE. C OUTPUT--THE SINGLE PRECISION CUMULATIVE DISTRIBUTION C FUNCTION VALUE CDF FOR THE GAMMA DISTRIBUTION C WITH TAIL LENGTH PARAMETER VALUE = GAMMA. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--GAMMA SHOULD BE POSITIVE. C --X SHOULD BE POSITIVE. C OTHER DATAPAC SUBROUTINES NEEDED--NONE. C FORTRAN LIBRARY SUBROUTINES NEEDED--DEXP, DLOG. C MODE OF INTERNAL OPERATIONS--DOUBLE PRECISION. C LANGUAGE--ANSI FORTRAN. C ACCURACY--(ON THE UNIVAC 1108, EXEC 8 SYSTEM AT NBS) C COMPARED TO THE KNOWN GAMMA = 1 (EXPONENTIAL) C RESULTS, AGREEMENT WAS HAD OUT TO 7 SIGNIFICANT C DIGITS FOR ALL TESTED X. C THE TESTED X VALUES COVERED THE ENTIRE C RANGE OF THE DISTRIBUTION--FROM THE 0.00001 C PERCENT POINT UP TO THE 99.99999 PERCENT POINT C OF THE DISTRIBUTION. C REFERENCES--WILK, GNANADESIKAN, AND HUYETT, 'PROBABILITY C PLOTS FOR THE GAMMA DISTRIBUTION', C TECHNOMETRICS, 1962, PAGES 1-15, C ESPECIALLY PAGES 3-5. C --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS C SERIES 55, 1964, PAGE 257, FORMULA 6.1.41. C --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE C DISTRIBUTIONS--1, 1970, PAGES 166-206. C --HASTINGS AND PEACOCK, STATISTICAL C DISTRIBUTIONS--A HANDBOOK FOR C STUDENTS AND PRACTITIONERS, 1975, C PAGES 68-73. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING LABORATORY (205.03) C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE: 301-921-2315 C ORIGINAL VERSION--NOVEMBER 1975. C C--------------------------------------------------------------------- C DOUBLE PRECISION DX,DGAMMA,AI,TERM,SUM,CUT1,CUT2,CUTOFF,T DOUBLE PRECISION Z,Z2,Z3,Z4,Z5,DEN,A,B,C,D,G DOUBLE PRECISION DEXP,DLOG DIMENSION D(10) DATA C/ .918938533204672741D0/ DATA D(1),D(2),D(3),D(4),D(5) 1 /+.833333333333333333D-1,-.277777777777777778D-2, 1+.793650793650793651D-3,-.595238095238095238D-3,+.8417508417508417 151D-3/ DATA D(6),D(7),D(8),D(9),D(10) 1 /-.191752691752691753D-2,+.641025641025641025D-2,-.2955065359 147712418D-1,+.179644372368830573D0,-.139243221690590111D1/ C C IPR=6 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(X.LE.0.0)GOTO50 IF(GAMMA.LE.0.0)GOTO55 GOTO90 50 WRITE(IPR,4) WRITE(IPR,46)X CDF=0.0 RETURN 55 WRITE(IPR,15) WRITE(IPR,46)GAMMA CDF=0.0 RETURN 90 CONTINUE 4 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUME 1NT TO THE GAMCDF SUBROUTINE IS NON-POSITIVE *****) 15 FORMAT(1H , 91H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 GAMCDF SUBROUTINE IS NON-POSITIVE *****) 46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****) C C-----START POINT----------------------------------------------------- C DX=X DGAMMA=GAMMA MAXIT=10000 C C COMPUTE THE GAMMA FUNCTION USING THE ALGORITHM IN THE C NBS APPLIED MATHEMATICS SERIES REFERENCE. C Z=DGAMMA DEN=1.0D0 300 IF(Z.GE.10.0D0)GOTO400 DEN=DEN*Z Z=Z+1 GOTO300 400 Z2=Z*Z Z3=Z*Z2 Z4=Z2*Z2 Z5=Z2*Z3 A=(Z-0.5D0)*DLOG(Z)-Z+C B=D(1)/Z+D(2)/Z3+D(3)/Z5+D(4)/(Z2*Z5)+D(5)/(Z4*Z5)+ 1D(6)/(Z*Z5*Z5)+D(7)/(Z3*Z5*Z5)+D(8)/(Z5*Z5*Z5)+D(9)/(Z2*Z5*Z5*Z5) G=DEXP(A+B)/DEN C C COMPUTE T-SUB-Q AS DEFINED ON PAGE 4 OF THE WILK, GNANADESIKAN, C AND HUYETT REFERENCE C SUM=1.0D0/DGAMMA TERM=1.0D0/DGAMMA CUT1=DX-DGAMMA CUT2=DX*10000000000.0D0 DO200I=1,MAXIT AI=I TERM=DX*TERM/(DGAMMA+AI) SUM=SUM+TERM CUTOFF=CUT1+(CUT2*TERM/SUM) IF(AI.GT.CUTOFF)GOTO250 200 CONTINUE WRITE(IPR,205)MAXIT WRITE(IPR,206)X WRITE(IPR,207)GAMMA WRITE(IPR,208) CDF=1.0 RETURN C 250 T=SUM CDF=(DX**DGAMMA)*(DEXP(-DX))*T/G C 205 FORMAT(1H ,48H*****ERROR IN INTERNAL OPERATIONS IN THE GAMCDF , 1 45HSUBROUTINE--THE NUMBER OF ITERATIONS EXCEEDS ,I7) 206 FORMAT(1H ,33H THE INPUT VALUE OF X IS ,E15.8) 207 FORMAT(1H ,33H THE INPUT VALUE OF GAMMA IS ,E15.8) 208 FORMAT(1H ,48H THE OUTPUT VALUE OF CDF HAS BEEN SET TO 1.0) C RETURN END C**** SUBROUTINE DSVRGN(N,RB,RA) C=============================================================================== C C HEAP SORT C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DIMENSION RA(N),RB(N) DO 5 I=1,N RA(I) = RB(I) 5 CONTINUE L=N/2+1 IR=N 10 CONTINUE IF(L.GT.1)THEN L=L-1 RRA=RA(L) ELSE RRA=RA(IR) RA(IR)=RA(1) IR=IR-1 IF(IR.EQ.1)THEN RA(1)=RRA RETURN ENDIF ENDIF I=L J=L+L 20 IF(J.LE.IR)THEN IF(J.LT.IR)THEN IF(RA(J).LT.RA(J+1))J=J+1 ENDIF IF(RRA.LT.RA(J))THEN RA(I)=RA(J) I=J J=J+J ELSE J=IR+1 ENDIF GO TO 20 ENDIF RA(I)=RRA GO TO 10 END C**** INTEGER FUNCTION NDAYS(IDAY, IMONTH, IYEAR) C=============================================================================== C C COMPUTE THE NUMBER OF DAYS SINCE 1 JAN 1900 (=0) C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE INTEGER IDAY, IMONTH, IYEAR C NDAYS0(1,1,1900) = 693961 NDAYS = NDAYS0(IDAY,IMONTH,IYEAR) - 693961 RETURN END C**** INTEGER FUNCTION NDAYS0(IDAY, IMONTH, IYEAR) C=============================================================================== C C COMPUTE THE NUMBER OF DAYS SINCE 1 JAN 0000 (=0) C C=============================================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE INTEGER IDAY, IMONTH, IYEAR, DM(12) DATA DM/0,31,59,90,120,151,181,212,243,273,304,334/ C C LEAP YEAR CORRECTIONS C 1. EVERY 4 YEARS C 2. EXCEPT EVERY 100 YEARS C 3. EXCEPT EVERY 400 YEARS C (I.E. LY'S IN 1896, 1904 1908, ..., 1996, 2000, C IF( IYEAR .LT. 1583) THEN WRITE(*,*) ' NDAYS FAILS FOR YEARS PRIOR TO 1583' STOP ENDIF NDAYS0 = 1 365*IYEAR + DM(IMONTH) + (IDAY - 1) + 2 IYEAR/4 - IYEAR/100 + IYEAR/400 + 1 IF(MOD(IYEAR,4) .EQ. 0 .AND. IMONTH .LT. 3) THEN IF(MOD(IYEAR,100) .NE. 0 .OR. MOD(IYEAR,400) .EQ. 0) THEN NDAYS0 = NDAYS0-1 ENDIF ENDIF RETURN END C**** C=============================================================================== C C DUMMY GRAPHICS ROUTINES, OUTPUT, ETC. C C=============================================================================== SUBROUTINE UMACH RETURN END SUBROUTINE DBOXP RETURN END SUBROUTINE PAGE RETURN END SUBROUTINE DPLOTP RETURN END C C=============================================================================== * ====================================================================== * NIST Guide to Available Math Software. * Fullsource for module GAMPPF from package DATAPAC. * Retrieved from CAMSUN on Fri Jul 10 07:59:09 1998. * ====================================================================== * GAMPPF SUBROUTINE GAMPPF(P,GAMMA,PPF) C C PURPOSE--THIS SUBROUTINE COMPUTES THE PERCENT POINT C FUNCTION VALUE FOR THE GAMMA DISTRIBUTION C WITH SINGLE PRECISION C TAIL LENGTH PARAMETER = GAMMA. C THE GAMMA DISTRIBUTION USED C HEREIN HAS MEAN = GAMMA C AND STANDARD DEVIATION = SQRT(GAMMA). C THIS DISTRIBUTION IS DEFINED FOR ALL POSITIVE X, C AND HAS THE PROBABILITY DENSITY FUNCTION C F(X) = (1/CONSTANT) * (X**(GAMMA-1)) * EXP(-X) C WHERE THE CONSTANT = THE GAMMA FUNCTION EVALUATED C AT THE VALUE GAMMA. C NOTE THAT THE PERCENT POINT FUNCTION OF A DISTRIBUTION C IS IDENTICALLY THE SAME AS THE INVERSE CUMULATIVE C DISTRIBUTION FUNCTION OF THE DISTRIBUTION. C INPUT ARGUMENTS--P = THE SINGLE PRECISION VALUE C (BETWEEN 0.0 (EXCLUSIVELY) C AND 1.0 (EXCLUSIVELY)) C AT WHICH THE PERCENT POINT C FUNCTION IS TO BE EVALUATED. C --GAMMA = THE SINGLE PRECISION VALUE OF THE C TAIL LENGTH PARAMETER. C GAMMA SHOULD BE POSITIVE. C OUTPUT ARGUMENTS--PPF = THE SINGLE PRECISION PERCENT C POINT FUNCTION VALUE. C OUTPUT--THE SINGLE PRECISION PERCENT POINT FUNCTION . C VALUE PPF FOR THE GAMMA DISTRIBUTION C WITH TAIL LENGTH PARAMETER VALUE = GAMMA. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--GAMMA SHOULD BE POSITIVE. C --P SHOULD BE BETWEEN 0.0 (EXCLUSIVELY) C AND 1.0 (EXCLUSIVELY). C OTHER DATAPAC SUBROUTINES NEEDED--NONE. C FORTRAN LIBRARY SUBROUTINES NEEDED--DEXP, DLOG. C MODE OF INTERNAL OPERATIONS--DOUBLE PRECISION. C LANGUAGE--ANSI FORTRAN. C ACCURACY--(ON THE UNIVAC 1108, EXEC 8 SYSTEM AT NBS) C COMPARED TO THE KNOWN GAMMA = 1 (EXPONENTIAL) C RESULTS, AGREEMENT WAS HAD OUT TO 6 SIGNIFICANT C DIGITS FOR ALL TESTED P IN THE RANGE P = .001 TO C P = .999. FOR P = .95 AND SMALLER, THE AGREEMENT C WAS EVEN BETTER--7 SIGNIFICANT DIGITS. C (NOTE THAT THE TABULATED VALUES GIVEN IN THE WILK, C GNANADESIKAN, AND HUYETT REFERENCE BELOW, PAGE 20, C ARE IN ERROR FOR AT LEAST THE GAMMA = 1 CASE-- C THE WORST DETECTED ERROR WAS AGREEMENT TO ONLY 3 C SIGNIFICANT DIGITS (IN THEIR 8 SIGNIFICANT DIGIT TABLE) C FOR P = .999.) C REFERENCES--WILK, GNANADESIKAN, AND HUYETT, 'PROBABILITY C PLOTS FOR THE GAMMA DISTRIBUTION', C TECHNOMETRICS, 1962, PAGES 1-15, C ESPECIALLY PAGES 3-5. C --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS C SERIES 55, 1964, PAGE 257, FORMULA 6.1.41. C --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE C DISTRIBUTIONS--1, 1970, PAGES 166-206. C --HASTINGS AND PEACOCK, STATISTICAL C DISTRIBUTIONS--A HANDBOOK FOR C STUDENTS AND PRACTITIONERS, 1975, C PAGES 68-73. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING LABORATORY (205.03) C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE: 301-921-2315 C ORIGINAL VERSION--NOVEMBER 1974. C UPDATED --SEPTEMBER 1975. C UPDATED --NOVEMBER 1975. C C--------------------------------------------------------------------- C DOUBLE PRECISION DP,DGAMMA DOUBLE PRECISION Z,Z2,Z3,Z4,Z5,DEN,A,B,C,D,G DOUBLE PRECISION XMIN0,XMIN,AI,XMAX,DX,PCALC,XMID DOUBLE PRECISION XLOWER,XUPPER,XDEL DOUBLE PRECISION SUM,TERM,CUT1,CUT2,AJ,CUTOFF,T DOUBLE PRECISION DEXP,DLOG DIMENSION D(10) DATA C/ .918938533204672741D0/ DATA D(1),D(2),D(3),D(4),D(5) 1 /+.833333333333333333D-1,-.277777777777777778D-2, 1+.793650793650793651D-3,-.595238095238095238D-3,+.8417508417508417 151D-3/ DATA D(6),D(7),D(8),D(9),D(10) 1 /-.191752691752691753D-2,+.641025641025641025D-2,-.2955065359 147712418D-1,+.179644372368830573D0,-.139243221690590111D1/ C IPR=6 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(P.LE.0.0.OR.P.GE.1.0)GOTO50 IF(GAMMA.LE.0.0)GOTO55 GOTO90 50 WRITE(IPR,1) WRITE(IPR,46)P PPF=0.0 RETURN 55 WRITE(IPR,15) WRITE(IPR,46)GAMMA PPF=0.0 RETURN 90 CONTINUE 1 FORMAT(1H ,115H***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO THE 1 GAMPPF SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****) 15 FORMAT(1H , 91H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 GAMPPF SUBROUTINE IS NON-POSITIVE *****) 46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****) C C-----START POINT----------------------------------------------------- C DP=P DGAMMA=GAMMA MAXIT=10000 C C COMPUTE THE GAMMA FUNCTION USING THE ALGORITHM IN THE C NBS APPLIED MATHEMATICS SERIES REFERENCE. C THIS GAMMA FUNCTION NEED BE CALCULATED ONLY ONCE. C IT IS USED IN THE CALCULATION OF THE CDF BASED ON C THE TENTATIVE VALUE OF THE PPF IN THE ITERATION. C Z=DGAMMA DEN=1.0D0 150 IF(Z.GE.10.0D0)GOTO160 DEN=DEN*Z Z=Z+1.0D0 GOTO150 160 Z2=Z*Z Z3=Z*Z2 Z4=Z2*Z2 Z5=Z2*Z3 A=(Z-0.5D0)*DLOG(Z)-Z+C B=D(1)/Z+D(2)/Z3+D(3)/Z5+D(4)/(Z2*Z5)+D(5)/(Z4*Z5)+ 1D(6)/(Z*Z5*Z5)+D(7)/(Z3*Z5*Z5)+D(8)/(Z5*Z5*Z5)+D(9)/(Z2*Z5*Z5*Z5) G=DEXP(A+B)/DEN C C DETERMINE LOWER AND UPPER LIMITS ON THE DESIRED 100P C PERCENT POINT. C ILOOP=1 XMIN0=(DP*DGAMMA*G)**(1.0D0/DGAMMA) XMIN=XMIN0 ICOUNT=1 350 AI=ICOUNT XMAX=AI*XMIN0 DX=XMAX GOTO1000 360 IF(PCALC.GE.DP)GOTO370 XMIN=XMAX ICOUNT=ICOUNT+1 IF(ICOUNT.LE.30000)GOTO350 370 XMID=(XMIN+XMAX)/2.0D0 C C NOW ITERATE BY BISECTION UNTIL THE DESIRED ACCURACY IS ACHIEVED. C ILOOP=2 XLOWER=XMIN XUPPER=XMAX ICOUNT=0 550 DX=XMID GOTO1000 560 IF(PCALC.EQ.DP)GOTO570 IF(PCALC.GT.DP)GOTO580 XLOWER=XMID XMID=(XMID+XUPPER)/2.0D0 GOTO590 580 XUPPER=XMID XMID=(XMID+XLOWER)/2.0D0 590 XDEL=XMID-XLOWER IF(XDEL.LT.0.0D0)XDEL=-XDEL ICOUNT=ICOUNT+1 IF(XDEL.LT.0.0000000001D0.OR.ICOUNT.GT.100)GOTO570 GOTO550 570 PPF=XMID RETURN C C******************************************************************** C THIS SECTION BELOW IS LOGICALLY SEPARATE FROM THE ABOVE. C THIS SECTION COMPUTES A CDF VALUE FOR ANY GIVEN TENTATIVE C PERCENT POINT X VALUE AS DEFINED IN EITHER OF THE 2 C ITERATION LOOPS IN THE ABOVE CODE. C C COMPUTE T-SUB-Q AS DEFINED ON PAGE 4 OF THE WILK, GNANADESIKAN, C AND HUYETT REFERENCE C 1000 SUM=1.0D0/DGAMMA TERM=1.0D0/DGAMMA CUT1=DX-DGAMMA CUT2=DX*10000000000.0D0 DO700J=1,MAXIT AJ=J TERM=DX*TERM/(DGAMMA+AJ) SUM=SUM+TERM CUTOFF=CUT1+(CUT2*TERM/SUM) IF(AJ.GT.CUTOFF)GOTO750 700 CONTINUE WRITE(IPR,705)MAXIT WRITE(IPR,706)P WRITE(IPR,707)GAMMA WRITE(IPR,708) PPF=0.0 RETURN C 750 T=SUM PCALC=(DX**DGAMMA)*(DEXP(-DX))*T/G IF(ILOOP.EQ.1)GOTO360 GOTO560 C 705 FORMAT(1H ,48H*****ERROR IN INTERNAL OPERATIONS IN THE GAMPPF , 1 45HSUBROUTINE--THE NUMBER OF ITERATIONS EXCEEDS ,I7) 706 FORMAT(1H ,33H THE INPUT VALUE OF P IS ,E15.8) 707 FORMAT(1H ,33H THE INPUT VALUE OF GAMMA IS ,E15.8) 708 FORMAT(1H ,48H THE OUTPUT VALUE OF PPF HAS BEEN SET TO 0.0) C END * ====================================================================== * NIST Guide to Available Math Software. * Fullsource for module DGAMMA from package CMLIB. * Retrieved from CAMSUN on Thu Jul 9 21:57:50 1998. * ====================================================================== DOUBLE PRECISION FUNCTION DGAMMA(X) C***BEGIN PROLOGUE DGAMMA C***DATE WRITTEN 770601 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. C7A C***KEYWORDS COMPLETE GAMMA FUNCTION,DOUBLE PRECISION,GAMMA FUNCTION, C SPECIAL FUNCTION C***AUTHOR FULLERTON, W., (LANL) C***PURPOSE Computes the d.p. complete Gamma function. C***DESCRIPTION C C DGAMMA(X) calculates the double precision complete gamma function C for double precision argument X. C C Series for GAM on the interval 0. to 1.00000E+00 C with weighted error 5.79E-32 C log weighted error 31.24 C significant figures required 30.00 C decimal places required 32.05 C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH,D9LGMC,DCSEVL,DGAMLM,DINT,INITDS,XERROR C***END PROLOGUE DGAMMA DOUBLE PRECISION X, GAMCS(42), DXREL, PI, SINPIY, SQ2PIL, XMAX, 1 XMIN, Y, DINT, D9LGMC, DCSEVL, D1MACH C DATA GAM CS( 1) / +.8571195590 9893314219 2006239994 2 D-2 / DATA GAM CS( 2) / +.4415381324 8410067571 9131577165 2 D-2 / DATA GAM CS( 3) / +.5685043681 5993633786 3266458878 9 D-1 / DATA GAM CS( 4) / -.4219835396 4185605010 1250018662 4 D-2 / DATA GAM CS( 5) / +.1326808181 2124602205 8400679635 2 D-2 / DATA GAM CS( 6) / -.1893024529 7988804325 2394702388 6 D-3 / DATA GAM CS( 7) / +.3606925327 4412452565 7808221722 5 D-4 / DATA GAM CS( 8) / -.6056761904 4608642184 8554829036 5 D-5 / DATA GAM CS( 9) / +.1055829546 3022833447 3182350909 3 D-5 / DATA GAM CS( 10) / -.1811967365 5423840482 9185589116 6 D-6 / DATA GAM CS( 11) / +.3117724964 7153222777 9025459316 9 D-7 / DATA GAM CS( 12) / -.5354219639 0196871408 7408102434 7 D-8 / DATA GAM CS( 13) / +.9193275519 8595889468 8778682594 0 D-9 / DATA GAM CS( 14) / -.1577941280 2883397617 6742327395 3 D-9 / DATA GAM CS( 15) / +.2707980622 9349545432 6654043308 9 D-10 / DATA GAM CS( 16) / -.4646818653 8257301440 8166105893 3 D-11 / DATA GAM CS( 17) / +.7973350192 0074196564 6076717535 9 D-12 / DATA GAM CS( 18) / -.1368078209 8309160257 9949917230 9 D-12 / DATA GAM CS( 19) / +.2347319486 5638006572 3347177168 8 D-13 / DATA GAM CS( 20) / -.4027432614 9490669327 6657053469 9 D-14 / DATA GAM CS( 21) / +.6910051747 3721009121 3833697525 7 D-15 / DATA GAM CS( 22) / -.1185584500 2219929070 5238712619 2 D-15 / DATA GAM CS( 23) / +.2034148542 4963739552 0102605193 2 D-16 / DATA GAM CS( 24) / -.3490054341 7174058492 7401294910 8 D-17 / DATA GAM CS( 25) / +.5987993856 4853055671 3505106602 6 D-18 / DATA GAM CS( 26) / -.1027378057 8722280744 9006977843 1 D-18 / DATA GAM CS( 27) / +.1762702816 0605298249 4275966074 8 D-19 / DATA GAM CS( 28) / -.3024320653 7353062609 5877211204 2 D-20 / DATA GAM CS( 29) / +.5188914660 2183978397 1783355050 6 D-21 / DATA GAM CS( 30) / -.8902770842 4565766924 4925160106 6 D-22 / DATA GAM CS( 31) / +.1527474068 4933426022 7459689130 6 D-22 / DATA GAM CS( 32) / -.2620731256 1873629002 5732833279 9 D-23 / DATA GAM CS( 33) / +.4496464047 8305386703 3104657066 6 D-24 / DATA GAM CS( 34) / -.7714712731 3368779117 0390152533 3 D-25 / DATA GAM CS( 35) / +.1323635453 1260440364 8657271466 6 D-25 / DATA GAM CS( 36) / -.2270999412 9429288167 0231381333 3 D-26 / DATA GAM CS( 37) / +.3896418998 0039914493 2081663999 9 D-27 / DATA GAM CS( 38) / -.6685198115 1259533277 9212799999 9 D-28 / DATA GAM CS( 39) / +.1146998663 1400243843 4761386666 6 D-28 / DATA GAM CS( 40) / -.1967938586 3451346772 9510399999 9 D-29 / DATA GAM CS( 41) / +.3376448816 5853380903 3489066666 6 D-30 / DATA GAM CS( 42) / -.5793070335 7821357846 2549333333 3 D-31 / DATA PI / 3.1415926535 8979323846 2643383279 50 D0 / DATA SQ2PIL / 0.9189385332 0467274178 0329736405 62 D0 / DATA NGAM, XMIN, XMAX, DXREL / 0, 3*0.D0 / C***FIRST EXECUTABLE STATEMENT DGAMMA IF (NGAM.NE.0) GO TO 10 NGAM = INITDS (GAMCS, 42, 0.1*SNGL(D1MACH(3)) ) C CALL DGAMLM (XMIN, XMAX) DXREL = DSQRT (D1MACH(4)) C 10 Y = DABS(X) IF (Y.GT.10.D0) GO TO 50 C C COMPUTE GAMMA(X) FOR -XBND .LE. X .LE. XBND. REDUCE INTERVAL AND FIND C GAMMA(1+Y) FOR 0.0 .LE. Y .LT. 1.0 FIRST OF ALL. C N = X IF (X.LT.0.D0) N = N - 1 Y = X - DBLE(FLOAT(N)) N = N - 1 DGAMMA = 0.9375D0 + DCSEVL (2.D0*Y-1.D0, GAMCS, NGAM) IF (N.EQ.0) RETURN C IF (N.GT.0) GO TO 30 C C COMPUTE GAMMA(X) FOR X .LT. 1.0 C N = -N IF (X.EQ.0.D0) CALL XERROR ( 'DGAMMA X IS 0', 14, 4, 2) IF (X.LT.0.0 .AND. X+DBLE(FLOAT(N-2)).EQ.0.D0) CALL XERROR ( 'DGAM 1MA X IS A NEGATIVE INTEGER', 31, 4, 2) IF (X.LT.(-0.5D0) .AND. DABS((X-DINT(X-0.5D0))/X).LT.DXREL) CALL 1 XERROR ( 'DGAMMA ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NE 2GATIVE INTEGER', 68, 1, 1) C DO 20 I=1,N DGAMMA = DGAMMA/(X+DBLE(FLOAT(I-1)) ) 20 CONTINUE RETURN C C GAMMA(X) FOR X .GE. 2.0 AND X .LE. 10.0 C 30 DO 40 I=1,N DGAMMA = (Y+DBLE(FLOAT(I))) * DGAMMA 40 CONTINUE RETURN C C GAMMA(X) FOR DABS(X) .GT. 10.0. RECALL Y = DABS(X). C 50 IF (X.GT.XMAX) CALL XERROR ( 'DGAMMA X SO BIG GAMMA OVERFLOWS', 1 32, 3, 2) C DGAMMA = 0.D0 IF (X.LT.XMIN) CALL XERROR ( 'DGAMMA X SO SMALL GAMMA UNDERFLOWS' 1 , 35, 2, 1) IF (X.LT.XMIN) RETURN C DGAMMA = DEXP ((Y-0.5D0)*DLOG(Y) - Y + SQ2PIL + D9LGMC(Y) ) IF (X.GT.0.D0) RETURN C IF (DABS((X-DINT(X-0.5D0))/X).LT.DXREL) CALL XERROR ( 'DGAMMA ANS 1WER LT HALF PRECISION, X TOO NEAR NEGATIVE INTEGER' , 61, 1, 1) C SINPIY = DSIN (PI*Y) IF (SINPIY.EQ.0.D0) CALL XERROR ( 'DGAMMA X IS A NEGATIVE INTEGER 1', 31, 4, 2) C DGAMMA = -PI/(Y*SINPIY*DGAMMA) C RETURN END FUNCTION INITDS(DOS,NOS,ETA) C***BEGIN PROLOGUE INITDS C***DATE WRITTEN 770601 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. C3A2 C***KEYWORDS CHEBYSHEV,DOUBLE PRECISION,INITIALIZE, C ORTHOGONAL POLYNOMIAL,SERIES,SPECIAL FUNCTION C***AUTHOR FULLERTON, W., (LANL) C***PURPOSE Initializes the d.p. properly normalized orthogonal C polynomial series to determine the number of terms needed C for specific accuracy. C***DESCRIPTION C C Initialize the double precision orthogonal series DOS so that INITDS C is the number of terms needed to insure the error is no larger than C ETA. Ordinarily ETA will be chosen to be one-tenth machine precision C C Input Arguments -- C DOS dble prec array of NOS coefficients in an orthogonal series. C NOS number of coefficients in DOS. C ETA requested accuracy of series. C***REFERENCES (NONE) C***ROUTINES CALLED XERROR C***END PROLOGUE INITDS C DOUBLE PRECISION DOS(NOS) C***FIRST EXECUTABLE STATEMENT INITDS IF (NOS.LT.1) CALL XERROR ( 'INITDS NUMBER OF COEFFICIENTS LT 1', 1 35, 2, 2) C ERR = 0. DO 10 II=1,NOS I = NOS + 1 - II ERR = ERR + ABS(SNGL(DOS(I))) IF (ERR.GT.ETA) GO TO 20 10 CONTINUE C 20 IF (I.EQ.NOS) CALL XERROR ( 'INITDS ETA MAY BE TOO SMALL', 28, 1 1, 2) INITDS = I C RETURN END SUBROUTINE XERROR(MESSG,NMESSG,NERR,LEVEL) C***BEGIN PROLOGUE XERROR C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Processes an error (diagnostic) message. C***DESCRIPTION C Abstract C XERROR processes a diagnostic message, in a manner C determined by the value of LEVEL and the current value C of the library error control flag, KONTRL. C (See subroutine XSETF for details.) C C Description of Parameters C --Input-- C MESSG - the Hollerith message to be processed, containing C no more than 72 characters. C NMESSG- the actual number of characters in MESSG. C NERR - the error number associated with this message. C NERR must not be zero. C LEVEL - error category. C =2 means this is an unconditionally fatal error. C =1 means this is a recoverable error. (I.e., it is C non-fatal if XSETF has been appropriately called.) C =0 means this is a warning message only. C =-1 means this is a warning message which is to be C printed at most once, regardless of how many C times this call is executed. C C Examples C CALL XERROR('SMOOTH -- NUM WAS ZERO.',23,1,2) C CALL XERROR('INTEG -- LESS THAN FULL ACCURACY ACHIEVED.', C 43,2,1) C CALL XERROR('ROOTER -- ACTUAL ZERO OF F FOUND BEFORE INTERVAL F C 1ULLY COLLAPSED.',65,3,0) C CALL XERROR('EXP -- UNDERFLOWS BEING SET TO ZERO.',39,1,-1) C C Latest revision --- 19 MAR 1980 C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED XERRWV C***END PROLOGUE XERROR CHARACTER*(*) MESSG C***FIRST EXECUTABLE STATEMENT XERROR CALL XERRWV(MESSG,NMESSG,NERR,LEVEL,0,0,0,0,0.,0.) RETURN END SUBROUTINE XERRWV(MESSG,NMESSG,NERR,LEVEL,NI,I1,I2,NR,R1,R2) C***BEGIN PROLOGUE XERRWV C***DATE WRITTEN 800319 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Processes error message allowing 2 integer and two real C values to be included in the message. C***DESCRIPTION C Abstract C XERRWV processes a diagnostic message, in a manner C determined by the value of LEVEL and the current value C of the library error control flag, KONTRL. C (See subroutine XSETF for details.) C In addition, up to two integer values and two real C values may be printed along with the message. C C Description of Parameters C --Input-- C MESSG - the Hollerith message to be processed. C NMESSG- the actual number of characters in MESSG. C NERR - the error number associated with this message. C NERR must not be zero. C LEVEL - error category. C =2 means this is an unconditionally fatal error. C =1 means this is a recoverable error. (I.e., it is C non-fatal if XSETF has been appropriately called.) C =0 means this is a warning message only. C =-1 means this is a warning message which is to be C printed at most once, regardless of how many C times this call is executed. C NI - number of integer values to be printed. (0 to 2) C I1 - first integer value. C I2 - second integer value. C NR - number of real values to be printed. (0 to 2) C R1 - first real value. C R2 - second real value. C C Examples C CALL XERRWV('SMOOTH -- NUM (=I1) WAS ZERO.',29,1,2, C 1 1,NUM,0,0,0.,0.) C CALL XERRWV('QUADXY -- REQUESTED ERROR (R1) LESS THAN MINIMUM ( C 1R2).,54,77,1,0,0,0,2,ERRREQ,ERRMIN) C C Latest revision --- 19 MAR 1980 C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED FDUMP,I1MACH,J4SAVE,XERABT,XERCTL,XERPRT,XERSAV, C XGETUA C***END PROLOGUE XERRWV CHARACTER*(*) MESSG CHARACTER*20 LFIRST CHARACTER*37 FORM DIMENSION LUN(5) C GET FLAGS C***FIRST EXECUTABLE STATEMENT XERRWV LKNTRL = J4SAVE(2,0,.FALSE.) MAXMES = J4SAVE(4,0,.FALSE.) C CHECK FOR VALID INPUT IF ((NMESSG.GT.0).AND.(NERR.NE.0).AND. 1 (LEVEL.GE.(-1)).AND.(LEVEL.LE.2)) GO TO 10 IF (LKNTRL.GT.0) CALL XERPRT('FATAL ERROR IN...',17) CALL XERPRT('XERROR -- INVALID INPUT',23) IF (LKNTRL.GT.0) CALL FDUMP IF (LKNTRL.GT.0) CALL XERPRT('JOB ABORT DUE TO FATAL ERROR.', 1 29) IF (LKNTRL.GT.0) CALL XERSAV(' ',0,0,0,KDUMMY) CALL XERABT('XERROR -- INVALID INPUT',23) RETURN 10 CONTINUE C RECORD MESSAGE JUNK = J4SAVE(1,NERR,.TRUE.) CALL XERSAV(MESSG,NMESSG,NERR,LEVEL,KOUNT) C LET USER OVERRIDE LFIRST = MESSG LMESSG = NMESSG LERR = NERR LLEVEL = LEVEL CALL XERCTL(LFIRST,LMESSG,LERR,LLEVEL,LKNTRL) C RESET TO ORIGINAL VALUES LMESSG = NMESSG LERR = NERR LLEVEL = LEVEL LKNTRL = MAX0(-2,MIN0(2,LKNTRL)) MKNTRL = IABS(LKNTRL) C DECIDE WHETHER TO PRINT MESSAGE IF ((LLEVEL.LT.2).AND.(LKNTRL.EQ.0)) GO TO 100 IF (((LLEVEL.EQ.(-1)).AND.(KOUNT.GT.MIN0(1,MAXMES))) 1.OR.((LLEVEL.EQ.0) .AND.(KOUNT.GT.MAXMES)) 2.OR.((LLEVEL.EQ.1) .AND.(KOUNT.GT.MAXMES).AND.(MKNTRL.EQ.1)) 3.OR.((LLEVEL.EQ.2) .AND.(KOUNT.GT.MAX0(1,MAXMES)))) GO TO 100 IF (LKNTRL.LE.0) GO TO 20 CALL XERPRT(' ',1) C INTRODUCTION IF (LLEVEL.EQ.(-1)) CALL XERPRT 1('WARNING MESSAGE...THIS MESSAGE WILL ONLY BE PRINTED ONCE.',57) IF (LLEVEL.EQ.0) CALL XERPRT('WARNING IN...',13) IF (LLEVEL.EQ.1) CALL XERPRT 1 ('RECOVERABLE ERROR IN...',23) IF (LLEVEL.EQ.2) CALL XERPRT('FATAL ERROR IN...',17) 20 CONTINUE C MESSAGE CALL XERPRT(MESSG,LMESSG) CALL XGETUA(LUN,NUNIT) ISIZEI = LOG10(FLOAT(I1MACH(9))) + 1.0 ISIZEF = LOG10(FLOAT(I1MACH(10))**I1MACH(11)) + 1.0 DO 50 KUNIT=1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) DO 22 I=1,MIN(NI,2) WRITE (FORM,21) I,ISIZEI 21 FORMAT ('(11X,21HIN ABOVE MESSAGE, I',I1,'=,I',I2,') ') IF (I.EQ.1) WRITE (IUNIT,FORM) I1 IF (I.EQ.2) WRITE (IUNIT,FORM) I2 22 CONTINUE DO 24 I=1,MIN(NR,2) WRITE (FORM,23) I,ISIZEF+10,ISIZEF 23 FORMAT ('(11X,21HIN ABOVE MESSAGE, R',I1,'=,E', 1 I2,'.',I2,')') IF (I.EQ.1) WRITE (IUNIT,FORM) R1 IF (I.EQ.2) WRITE (IUNIT,FORM) R2 24 CONTINUE IF (LKNTRL.LE.0) GO TO 40 C ERROR NUMBER WRITE (IUNIT,30) LERR 30 FORMAT (15H ERROR NUMBER =,I10) 40 CONTINUE 50 CONTINUE C TRACE-BACK IF (LKNTRL.GT.0) CALL FDUMP 100 CONTINUE IFATAL = 0 IF ((LLEVEL.EQ.2).OR.((LLEVEL.EQ.1).AND.(MKNTRL.EQ.2))) 1IFATAL = 1 C QUIT HERE IF MESSAGE IS NOT FATAL IF (IFATAL.LE.0) RETURN IF ((LKNTRL.LE.0).OR.(KOUNT.GT.MAX0(1,MAXMES))) GO TO 120 C PRINT REASON FOR ABORT IF (LLEVEL.EQ.1) CALL XERPRT 1 ('JOB ABORT DUE TO UNRECOVERED ERROR.',35) IF (LLEVEL.EQ.2) CALL XERPRT 1 ('JOB ABORT DUE TO FATAL ERROR.',29) C PRINT ERROR SUMMARY CALL XERSAV(' ',-1,0,0,KDUMMY) 120 CONTINUE C ABORT IF ((LLEVEL.EQ.2).AND.(KOUNT.GT.MAX0(1,MAXMES))) LMESSG = 0 CALL XERABT(MESSG,LMESSG) RETURN END SUBROUTINE XERSAV(MESSG,NMESSG,NERR,LEVEL,ICOUNT) C***BEGIN PROLOGUE XERSAV C***DATE WRITTEN 800319 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. Z C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Records that an error occurred. C***DESCRIPTION C Abstract C Record that this error occurred. C C Description of Parameters C --Input-- C MESSG, NMESSG, NERR, LEVEL are as in XERROR, C except that when NMESSG=0 the tables will be C dumped and cleared, and when NMESSG is less than zero the C tables will be dumped and not cleared. C --Output-- C ICOUNT will be the number of times this message has C been seen, or zero if the table has overflowed and C does not contain this message specifically. C When NMESSG=0, ICOUNT will not be altered. C C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C Latest revision --- 19 Mar 1980 C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED I1MACH,S88FMT,XGETUA C***END PROLOGUE XERSAV INTEGER LUN(5) CHARACTER*(*) MESSG CHARACTER*20 MESTAB(10),MES DIMENSION NERTAB(10),LEVTAB(10),KOUNT(10) SAVE MESTAB,NERTAB,LEVTAB,KOUNT,KOUNTX C NEXT TWO DATA STATEMENTS ARE NECESSARY TO PROVIDE A BLANK C ERROR TABLE INITIALLY DATA KOUNT(1),KOUNT(2),KOUNT(3),KOUNT(4),KOUNT(5), 1 KOUNT(6),KOUNT(7),KOUNT(8),KOUNT(9),KOUNT(10) 2 /0,0,0,0,0,0,0,0,0,0/ DATA KOUNTX/0/ C***FIRST EXECUTABLE STATEMENT XERSAV IF (NMESSG.GT.0) GO TO 80 C DUMP THE TABLE IF (KOUNT(1).EQ.0) RETURN C PRINT TO EACH UNIT CALL XGETUA(LUN,NUNIT) DO 60 KUNIT=1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) C PRINT TABLE HEADER WRITE (IUNIT,10) 10 FORMAT (32H0 ERROR MESSAGE SUMMARY/ 1 51H MESSAGE START NERR LEVEL COUNT) C PRINT BODY OF TABLE DO 20 I=1,10 IF (KOUNT(I).EQ.0) GO TO 30 WRITE (IUNIT,15) MESTAB(I),NERTAB(I),LEVTAB(I),KOUNT(I) 15 FORMAT (1X,A20,3I10) 20 CONTINUE 30 CONTINUE C PRINT NUMBER OF OTHER ERRORS IF (KOUNTX.NE.0) WRITE (IUNIT,40) KOUNTX 40 FORMAT (41H0OTHER ERRORS NOT INDIVIDUALLY TABULATED=,I10) WRITE (IUNIT,50) 50 FORMAT (1X) 60 CONTINUE IF (NMESSG.LT.0) RETURN C CLEAR THE ERROR TABLES DO 70 I=1,10 70 KOUNT(I) = 0 KOUNTX = 0 RETURN 80 CONTINUE C PROCESS A MESSAGE... C SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG, C OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL. MES = MESSG DO 90 I=1,10 II = I IF (KOUNT(I).EQ.0) GO TO 110 IF (MES.NE.MESTAB(I)) GO TO 90 IF (NERR.NE.NERTAB(I)) GO TO 90 IF (LEVEL.NE.LEVTAB(I)) GO TO 90 GO TO 100 90 CONTINUE C THREE POSSIBLE CASES... C TABLE IS FULL KOUNTX = KOUNTX+1 ICOUNT = 1 RETURN C MESSAGE FOUND IN TABLE 100 KOUNT(II) = KOUNT(II) + 1 ICOUNT = KOUNT(II) RETURN C EMPTY SLOT FOUND FOR NEW MESSAGE 110 MESTAB(II) = MES NERTAB(II) = NERR LEVTAB(II) = LEVEL KOUNT(II) = 1 ICOUNT = 1 RETURN END SUBROUTINE XGETUA(IUNITA,N) C***BEGIN PROLOGUE XGETUA C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Returns unit number(s) to which error messages are being C sent. C***DESCRIPTION C Abstract C XGETUA may be called to determine the unit number or numbers C to which error messages are being sent. C These unit numbers may have been set by a call to XSETUN, C or a call to XSETUA, or may be a default value. C C Description of Parameters C --Output-- C IUNIT - an array of one to five unit numbers, depending C on the value of N. A value of zero refers to the C default unit, as defined by the I1MACH machine C constant routine. Only IUNIT(1),...,IUNIT(N) are C defined by XGETUA. The values of IUNIT(N+1),..., C IUNIT(5) are not defined (for N .LT. 5) or altered C in any way by XGETUA. C N - the number of units to which copies of the C error messages are being sent. N will be in the C range from 1 to 5. C C Latest revision --- 19 MAR 1980 C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED J4SAVE C***END PROLOGUE XGETUA DIMENSION IUNITA(5) C***FIRST EXECUTABLE STATEMENT XGETUA N = J4SAVE(5,0,.FALSE.) DO 30 I=1,N INDEX = I+4 IF (I.EQ.1) INDEX = 3 IUNITA(I) = J4SAVE(INDEX,0,.FALSE.) 30 CONTINUE RETURN END DOUBLE PRECISION FUNCTION D1MACH(I) C***BEGIN PROLOGUE D1MACH C***DATE WRITTEN 750101 (YYMMDD) C***REVISION DATE 910131 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE Returns double precision machine dependent constants C***DESCRIPTION C C This is the CMLIB version of D1MACH, the double precision machine C constants subroutine originally developed for the PORT library. C C D1MACH can be used to obtain machine-dependent parameters C for the local machine environment. It is a function C subprogram with one (input) argument, and can be called C as follows, for example C C D = D1MACH(I) C C where I=1,...,5. The (output) value of D above is C determined by the (input) value of I. The results for C various values of I are discussed below. C C Double-precision machine constants C D1MACH( 1) = B**(EMIN-1), the smallest positive magnitude. C D1MACH( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude. C D1MACH( 3) = B**(-T), the smallest relative spacing. C D1MACH( 4) = B**(1-T), the largest relative spacing. C D1MACH( 5) = LOG10(B) C***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A C PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. C***ROUTINES CALLED XERROR C***END PROLOGUE D1MACH C INTEGER SMALL(4) INTEGER LARGE(4) INTEGER RIGHT(4) INTEGER DIVER(4) INTEGER LOG10(4) C DOUBLE PRECISION DMACH(5) C EQUIVALENCE (DMACH(1),SMALL(1)) EQUIVALENCE (DMACH(2),LARGE(1)) EQUIVALENCE (DMACH(3),RIGHT(1)) EQUIVALENCE (DMACH(4),DIVER(1)) EQUIVALENCE (DMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T C 3B SERIES AND MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T C PC 7300), IN WHICH THE MOST SIGNIFICANT BYTE IS STORED FIRST. C C === MACHINE = IEEE.MOST-SIG-BYTE-FIRST C === MACHINE = SUN C === MACHINE = 68000 C === MACHINE = ATT.3B C === MACHINE = ATT.7300 DATA SMALL(1),SMALL(2) / 1048576, 0 / DATA LARGE(1),LARGE(2) / 2146435071, -1 / DATA RIGHT(1),RIGHT(2) / 1017118720, 0 / DATA DIVER(1),DIVER(2) / 1018167296, 0 / DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 / C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES AND 8087-BASED C MICROS, SUCH AS THE IBM PC AND AT&T 6300, IN WHICH THE LEAST C SIGNIFICANT BYTE IS STORED FIRST. C C === MACHINE = IEEE.LEAST-SIG-BYTE-FIRST C === MACHINE = 8087 C === MACHINE = IBM.PC C === MACHINE = ATT.6300 C DATA SMALL(1),SMALL(2) / 0, 1048576 / C DATA LARGE(1),LARGE(2) / -1, 2146435071 / C DATA RIGHT(1),RIGHT(2) / 0, 1017118720 / C DATA DIVER(1),DIVER(2) / 0, 1018167296 / C DATA LOG10(1),LOG10(2) / 1352628735, 1070810131 / C C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C === MACHINE = AMDAHL C DATA SMALL(1),SMALL(2) / 1048576, 0 / C DATA LARGE(1),LARGE(2) / 2147483647, -1 / C DATA RIGHT(1),RIGHT(2) / 856686592, 0 / C DATA DIVER(1),DIVER(2) / 873463808, 0 / C DATA LOG10(1),LOG10(2) / 1091781651, 1352628735 / C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C === MACHINE = BURROUGHS.1700 C DATA SMALL(1) / ZC00800000 / C DATA SMALL(2) / Z000000000 / C DATA LARGE(1) / ZDFFFFFFFF / C DATA LARGE(2) / ZFFFFFFFFF / C DATA RIGHT(1) / ZCC5800000 / C DATA RIGHT(2) / Z000000000 / C DATA DIVER(1) / ZCC6800000 / C DATA DIVER(2) / Z000000000 / C DATA LOG10(1) / ZD00E730E7 / C DATA LOG10(2) / ZC77800DC0 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C === MACHINE = BURROUGHS.5700 C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O0000000000000000 / C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O0007777777777777 / C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C === MACHINE = BURROUGHS.6700 C === MACHINE = BURROUGHS.7700 C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O7770000000000000 / C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O7777777777777777 / C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE CONVEX C-120 (NATIVE MODE) C WITH OR WITHOUT -R8 OPTION C C === MACHINE = CONVEX.C1 C === MACHINE = CONVEX.C1.R8 C DATA DMACH(1) / 5.562684646268007D-309 / C DATA DMACH(2) / 8.988465674311577D+307 / C DATA DMACH(3) / 1.110223024625157D-016 / C DATA DMACH(4) / 2.220446049250313D-016 / C DATA DMACH(5) / 3.010299956639812D-001 / C C MACHINE CONSTANTS FOR THE CONVEX C-120 (IEEE MODE) C WITH OR WITHOUT -R8 OPTION C C === MACHINE = CONVEX.C1.IEEE C === MACHINE = CONVEX.C1.IEEE.R8 C DATA DMACH(1) / 2.225073858507202D-308 / C DATA DMACH(2) / 1.797693134862315D+308 / C DATA DMACH(3) / 1.110223024625157D-016 / C DATA DMACH(4) / 2.220446049250313D-016 / C DATA DMACH(5) / 3.010299956639812D-001 / C C MACHINE CONSTANTS FOR THE CYBER 170/180 SERIES USING NOS (FTN5). C C === MACHINE = CYBER.170.NOS C === MACHINE = CYBER.180.NOS C DATA SMALL(1) / O"00604000000000000000" / C DATA SMALL(2) / O"00000000000000000000" / C DATA LARGE(1) / O"37767777777777777777" / C DATA LARGE(2) / O"37167777777777777777" / C DATA RIGHT(1) / O"15604000000000000000" / C DATA RIGHT(2) / O"15000000000000000000" / C DATA DIVER(1) / O"15614000000000000000" / C DATA DIVER(2) / O"15010000000000000000" / C DATA LOG10(1) / O"17164642023241175717" / C DATA LOG10(2) / O"16367571421742254654" / C C MACHINE CONSTANTS FOR THE CDC 180 SERIES USING NOS/VE C C === MACHINE = CYBER.180.NOS/VE C DATA SMALL(1) / Z"3001800000000000" / C DATA SMALL(2) / Z"3001000000000000" / C DATA LARGE(1) / Z"4FFEFFFFFFFFFFFE" / C DATA LARGE(2) / Z"4FFE000000000000" / C DATA RIGHT(1) / Z"3FD2800000000000" / C DATA RIGHT(2) / Z"3FD2000000000000" / C DATA DIVER(1) / Z"3FD3800000000000" / C DATA DIVER(2) / Z"3FD3000000000000" / C DATA LOG10(1) / Z"3FFF9A209A84FBCF" / C DATA LOG10(2) / Z"3FFFF7988F8959AC" / C C MACHINE CONSTANTS FOR THE CYBER 205 C C === MACHINE = CYBER.205 C DATA SMALL(1) / X'9000400000000000' / C DATA SMALL(2) / X'8FD1000000000000' / C DATA LARGE(1) / X'6FFF7FFFFFFFFFFF' / C DATA LARGE(2) / X'6FD07FFFFFFFFFFF' / C DATA RIGHT(1) / X'FF74400000000000' / C DATA RIGHT(2) / X'FF45000000000000' / C DATA DIVER(1) / X'FF75400000000000' / C DATA DIVER(2) / X'FF46000000000000' / C DATA LOG10(1) / X'FFD04D104D427DE7' / C DATA LOG10(2) / X'FFA17DE623E2566A' / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C === MACHINE = CDC.6000 C === MACHINE = CDC.7000 C DATA SMALL(1) / 00604000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C DATA LARGE(1) / 37767777777777777777B / C DATA LARGE(2) / 37167777777777777777B / C DATA RIGHT(1) / 15604000000000000000B / C DATA RIGHT(2) / 15000000000000000000B / C DATA DIVER(1) / 15614000000000000000B / C DATA DIVER(2) / 15010000000000000000B / C DATA LOG10(1) / 17164642023241175717B / C DATA LOG10(2) / 16367571421742254654B / C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C C === MACHINE = CRAY.46-BIT-INTEGER C === MACHINE = CRAY.64-BIT-INTEGER C DATA SMALL(1) / 201354000000000000000B / C DATA SMALL(2) / 000000000000000000000B / C DATA LARGE(1) / 577767777777777777777B / C DATA LARGE(2) / 000007777777777777776B / C DATA RIGHT(1) / 376434000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C DATA DIVER(1) / 376444000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE - C STATIC DMACH(5) C C === MACHINE = DATA_GENERAL.ECLIPSE.S/200 C DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/ C DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/ C DATA LOG10/40423K,42023K,50237K,74776K/ C C ELXSI 6400 C C === MACHINE = ELSXI.6400 C DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X / C DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF'X / C DATA RIGHT(1), RIGHT(2) / '3CB00000'X,'00000000'X / C DATA DIVER(1), DIVER(2) / '3CC00000'X,'00000000'X / C DATA LOG10(1), DIVER(2) / '3FD34413'X,'509F79FF'X / C C MACHINE CONSTANTS FOR THE HARRIS 220 C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7 C C === MACHINE = HARRIS.220 C === MACHINE = HARRIS.SLASH6 C === MACHINE = HARRIS.SLASH7 C DATA SMALL(1),SMALL(2) / '20000000, '00000201 / C DATA LARGE(1),LARGE(2) / '37777777, '37777577 / C DATA RIGHT(1),RIGHT(2) / '20000000, '00000333 / C DATA DIVER(1),DIVER(2) / '20000000, '00000334 / C DATA LOG10(1),LOG10(2) / '23210115, '10237777 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. C C === MACHINE = HONEYWELL.600/6000 C === MACHINE = HONEYWELL.DPS.8/70 C DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 / C C MACHINE CONSTANTS FOR THE HP 2100 C 3 WORD DOUBLE PRECISION OPTION WITH FTN4 C C === MACHINE = HP.2100.3_WORD_DP C DATA SMALL(1), SMALL(2), SMALL(3) / 40000B, 0, 1 / C DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B / C DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B, 0, 265B / C DATA DIVER(1), DIVER(2), DIVER(3) / 40000B, 0, 276B / C DATA LOG10(1), LOG10(2), LOG10(3) / 46420B, 46502B, 77777B / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION OPTION WITH FTN4 C C === MACHINE = HP.2100.4_WORD_DP C DATA SMALL(1), SMALL(2) / 40000B, 0 / C DATA SMALL(3), SMALL(4) / 0, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177777B / C DATA LARGE(3), LARGE(4) / 177777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 0 / C DATA RIGHT(3), RIGHT(4) / 0, 225B / C DATA DIVER(1), DIVER(2) / 40000B, 0 / C DATA DIVER(3), DIVER(4) / 0, 227B / C DATA LOG10(1), LOG10(2) / 46420B, 46502B / C DATA LOG10(3), LOG10(4) / 76747B, 176377B / C C HP 9000 C C D1MACH(1) = 2.8480954D-306 C D1MACH(2) = 1.40444776D+306 C D1MACH(3) = 2.22044605D-16 C D1MACH(4) = 4.44089210D-16 C D1MACH(5) = 3.01029996D-1 C C === MACHINE = HP.9000 C DATA SMALL(1), SMALL(2) / 00040000000B, 00000000000B / C DATA LARGE(1), LARGE(2) / 17737777777B, 37777777777B / C DATA RIGHT(1), RIGHT(2) / 07454000000B, 00000000000B / C DATA DIVER(1), DIVER(2) / 07460000000B, 00000000000B / C DATA LOG10(1), LOG10(2) / 07764642023B, 12047674777B / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND C THE INTERDATA 3230 AND INTERDATA 7/32. C C === MACHINE = IBM.360 C === MACHINE = IBM.370 C === MACHINE = XEROX.SIGMA.5 C === MACHINE = XEROX.SIGMA.7 C === MACHINE = XEROX.SIGMA.9 C === MACHINE = SEL.85 C === MACHINE = SEL.86 C === MACHINE = INTERDATA.3230 C === MACHINE = INTERDATA.7/32 C DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 / C DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / C DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 / C DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 / C DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF / C C MACHINE CONSTANTS FOR THE INTERDATA 8/32 C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. C C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. C C === MACHINE = INTERDATA.8/32.UNIX C DATA SMALL(1),SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1),LARGE(2) / Z'7EFFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1),RIGHT(2) / Z'33100000', Z'00000000' / C DATA DIVER(1),DIVER(2) / Z'34100000', Z'00000000' / C DATA LOG10(1),LOG10(2) / Z'41134413', Z'509F79FF' / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C === MACHINE = PDP-10.KA C DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 / C DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C === MACHINE = PDP-10.KI C DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 / C DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "047674776746 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C === MACHINE = PDP-11.32-BIT C DATA SMALL(1),SMALL(2) / 8388608, 0 / C DATA LARGE(1),LARGE(2) / 2147483647, -1 / C DATA RIGHT(1),RIGHT(2) / 612368384, 0 / C DATA DIVER(1),DIVER(2) / 620756992, 0 / C DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 / C C DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 / C DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 / C DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 / C DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 / C DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C === MACHINE = PDP-11.16-BIT C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA SMALL(3),SMALL(4) / 0, 0 / C DATA LARGE(1),LARGE(2) / 32767, -1 / C DATA LARGE(3),LARGE(4) / -1, -1 / C DATA RIGHT(1),RIGHT(2) / 9344, 0 / C DATA RIGHT(3),RIGHT(4) / 0, 0 / C DATA DIVER(1),DIVER(2) / 9472, 0 / C DATA DIVER(3),DIVER(4) / 0, 0 / C DATA LOG10(1),LOG10(2) / 16282, 8346 / C DATA LOG10(3),LOG10(4) / -31493, -12296 / C C DATA SMALL(1),SMALL(2) / O000200, O000000 / C DATA SMALL(3),SMALL(4) / O000000, O000000 / C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA LARGE(3),LARGE(4) / O177777, O177777 / C DATA RIGHT(1),RIGHT(2) / O022200, O000000 / C DATA RIGHT(3),RIGHT(4) / O000000, O000000 / C DATA DIVER(1),DIVER(2) / O022400, O000000 / C DATA DIVER(3),DIVER(4) / O000000, O000000 / C DATA LOG10(1),LOG10(2) / O037632, O020232 / C DATA LOG10(3),LOG10(4) / O102373, O147770 / C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000 C C === MACHINE = SEQUENT.BALANCE.8000 C DATA SMALL(1),SMALL(2) / $00000000, $00100000 / C DATA LARGE(1),LARGE(2) / $FFFFFFFF, $7FEFFFFF / C DATA RIGHT(1),RIGHT(2) / $00000000, $3CA00000 / C DATA DIVER(1),DIVER(2) / $00000000, $3CB00000 / C DATA LOG10(1),LOG10(2) / $509F79FF, $3FD34413 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER C C === MACHINE = UNIVAC.1100 C DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 / C C MACHINE CONSTANTS FOR VAX 11/780 C (EXPRESSED IN INTEGER AND HEXADECIMAL) C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C === MACHINE = VAX.11/780 C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 9344, 0 / C DATA DIVER(1), DIVER(2) / 9472, 0 / C DATA LOG10(1), LOG10(2) / 546979738, -805796613 / C C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS*** C DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 / C DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB / C C MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING) C (EXPRESSED IN INTEGER AND HEXADECIMAL) C *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS*** C C DATA SMALL(1), SMALL(2) / 16, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 15552, 0 / C DATA DIVER(1), DIVER(2) / 15568, 0 / C DATA LOG10(1), LOG10(2) / 1142112243, 2046775455 / C C ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS*** C DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 / C DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F / C C C***FIRST EXECUTABLE STATEMENT D1MACH IF (I .LT. 1 .OR. I .GT. 5) 1 CALL XERROR( 'D1MACH -- I OUT OF BOUNDS',25,1,2) C D1MACH = DMACH(I) RETURN C END DOUBLE PRECISION FUNCTION D9LGMC(X) C***BEGIN PROLOGUE D9LGMC C***DATE WRITTEN 770601 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. C7E C***KEYWORDS COMPLETE GAMMA FUNCTION,CORRECTION FACTOR, C DOUBLE PRECISION,GAMMA FUNCTION,LOGARITHM, C SPECIAL FUNCTION C***AUTHOR FULLERTON, W., (LANL) C***PURPOSE Computes the d.p. log Gamma correction factor for C X .GE. 10. so that DLOG(DGAMMA(X)) = DLOG(DSQRT(2*PI)) + C (X-5.)*DLOG(X) - X + D9LGMC(X) C***DESCRIPTION C C Compute the log gamma correction factor for X .GE. 10. so that C DLOG (DGAMMA(X)) = DLOG(DSQRT(2*PI)) + (X-.5)*DLOG(X) - X + D9lGMC(X) C C Series for ALGM on the interval 0. to 1.00000E-02 C with weighted error 1.28E-31 C log weighted error 30.89 C significant figures required 29.81 C decimal places required 31.48 C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH,DCSEVL,INITDS,XERROR C***END PROLOGUE D9LGMC DOUBLE PRECISION X, ALGMCS(15), XBIG, XMAX, DCSEVL, D1MACH DATA ALGMCS( 1) / +.1666389480 4518632472 0572965082 2 D+0 / DATA ALGMCS( 2) / -.1384948176 0675638407 3298605913 5 D-4 / DATA ALGMCS( 3) / +.9810825646 9247294261 5717154748 7 D-8 / DATA ALGMCS( 4) / -.1809129475 5724941942 6330626671 9 D-10 / DATA ALGMCS( 5) / +.6221098041 8926052271 2601554341 6 D-13 / DATA ALGMCS( 6) / -.3399615005 4177219443 0333059966 6 D-15 / DATA ALGMCS( 7) / +.2683181998 4826987489 5753884666 6 D-17 / DATA ALGMCS( 8) / -.2868042435 3346432841 4462239999 9 D-19 / DATA ALGMCS( 9) / +.3962837061 0464348036 7930666666 6 D-21 / DATA ALGMCS( 10) / -.6831888753 9857668701 1199999999 9 D-23 / DATA ALGMCS( 11) / +.1429227355 9424981475 7333333333 3 D-24 / DATA ALGMCS( 12) / -.3547598158 1010705471 9999999999 9 D-26 / DATA ALGMCS( 13) / +.1025680058 0104709120 0000000000 0 D-27 / DATA ALGMCS( 14) / -.3401102254 3167487999 9999999999 9 D-29 / DATA ALGMCS( 15) / +.1276642195 6300629333 3333333333 3 D-30 / DATA NALGM, XBIG, XMAX / 0, 2*0.D0 / C***FIRST EXECUTABLE STATEMENT D9LGMC IF (NALGM.NE.0) GO TO 10 NALGM = INITDS (ALGMCS, 15, SNGL(D1MACH(3)) ) XBIG = 1.0D0/DSQRT(D1MACH(3)) XMAX = DEXP (DMIN1(DLOG(D1MACH(2)/12.D0), -DLOG(12.D0*D1MACH(1)))) C 10 IF (X.LT.10.D0) CALL XERROR ( 'D9LGMC X MUST BE GE 10', 23, 1, 2) IF (X.GE.XMAX) GO TO 20 C D9LGMC = 1.D0/(12.D0*X) IF (X.LT.XBIG) D9LGMC = DCSEVL (2.0D0*(10.D0/X)**2-1.D0, ALGMCS, 1 NALGM) / X RETURN C 20 D9LGMC = 0.D0 CALL XERROR ( 'D9LGMC X SO BIG D9LGMC UNDERFLOWS', 34, 2, 1) RETURN C END DOUBLE PRECISION FUNCTION DCSEVL(X,A,N) C***BEGIN PROLOGUE DCSEVL C***DATE WRITTEN 770401 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. C3A2 C***KEYWORDS CHEBYSHEV,FNLIB,SPECIAL FUNCTION C***AUTHOR FULLERTON, W., (LANL) C***PURPOSE Evaluate the double precision N-term Chebyshev series A C at X. C***DESCRIPTION C C Evaluate the N-term Chebyshev series A at X. Adapted from C R. Broucke, Algorithm 446, C.A.C.M., 16, 254 (1973). C W. Fullerton, C-3, Los Alamos Scientific Laboratory. C C Input Arguments -- C X double precision value at which the series is to be evaluated. C A double precision array of N terms of a Chebyshev series. In C evaluating A, only half of the first coefficient is summed. C N number of terms in array A. C***REFERENCES (NONE) C***ROUTINES CALLED XERROR C***END PROLOGUE DCSEVL C DOUBLE PRECISION A(N),X,TWOX,B0,B1,B2 C***FIRST EXECUTABLE STATEMENT DCSEVL IF(N.LT.1)CALL XERROR( 'DCSEVL NUMBER OF TERMS LE 0', 28, 2,2) IF(N.GT.1000) CALL XERROR ( 'DCSEVL NUMBER OF TERMS GT 1000', 1 31, 3, 2) IF ((X.LT.-1.D0) .OR. (X.GT.1.D0)) CALL XERROR ( 'DCSEVL X OUTSI 1DE (-1,+1)', 25, 1, 1) C TWOX = 2.0D0*X B1 = 0.D0 B0=0.D0 DO 10 I=1,N B2=B1 B1=B0 NI = N - I + 1 B0 = TWOX*B1 - B2 + A(NI) 10 CONTINUE C DCSEVL = 0.5D0 * (B0-B2) C RETURN END SUBROUTINE DGAMLM(XMIN,XMAX) C***BEGIN PROLOGUE DGAMLM C***DATE WRITTEN 770601 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. C7A,R2 C***KEYWORDS COMPLETE GAMMA FUNCTION,DOUBLE PRECISION,GAMMA FUNCTION, C LIMITS,SPECIAL FUNCTION C***AUTHOR FULLERTON, W., (LANL) C***PURPOSE Computes the d.p. minimum and maximum bounds for X in C GAMMA(X). C***DESCRIPTION C C Calculate the minimum and maximum legal bounds for X in gamma(X). C XMIN and XMAX are not the only bounds, but they are the only non- C trivial ones to calculate. C C Output Arguments -- C XMIN double precision minimum legal value of X in gamma(X). Any C smaller value of X might result in underflow. C XMAX double precision maximum legal value of X in gamma(X). Any C larger value of X might cause overflow. C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH,XERROR C***END PROLOGUE DGAMLM DOUBLE PRECISION XMIN, XMAX, ALNBIG, ALNSML, XLN, XOLD, D1MACH C***FIRST EXECUTABLE STATEMENT DGAMLM ALNSML = DLOG(D1MACH(1)) XMIN = -ALNSML DO 10 I=1,10 XOLD = XMIN XLN = DLOG(XMIN) XMIN = XMIN - XMIN*((XMIN+0.5D0)*XLN - XMIN - 0.2258D0 + ALNSML) 1 / (XMIN*XLN+0.5D0) IF (DABS(XMIN-XOLD).LT.0.005D0) GO TO 20 10 CONTINUE CALL XERROR ( 'DGAMLM UNABLE TO FIND XMIN', 27, 1, 2) C 20 XMIN = -XMIN + 0.01D0 C ALNBIG = DLOG (D1MACH(2)) XMAX = ALNBIG DO 30 I=1,10 XOLD = XMAX XLN = DLOG(XMAX) XMAX = XMAX - XMAX*((XMAX-0.5D0)*XLN - XMAX + 0.9189D0 - ALNBIG) 1 / (XMAX*XLN-0.5D0) IF (DABS(XMAX-XOLD).LT.0.005D0) GO TO 40 30 CONTINUE CALL XERROR ( 'DGAMLM UNABLE TO FIND XMAX', 27, 2, 2) C 40 XMAX = XMAX - 0.01D0 XMIN = DMAX1 (XMIN, -XMAX+1.D0) C RETURN END SUBROUTINE FDUMP C***BEGIN PROLOGUE FDUMP C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. Z C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Symbolic dump (should be locally written). C***DESCRIPTION C ***Note*** Machine Dependent Routine C FDUMP is intended to be replaced by a locally written C version which produces a symbolic dump. Failing this, C it should be replaced by a version which prints the C subprogram nesting list. Note that this dump must be C printed on each of up to five files, as indicated by the C XGETUA routine. See XSETUA and XGETUA for details. C C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C Latest revision --- 23 May 1979 C***ROUTINES CALLED (NONE) C***END PROLOGUE FDUMP C***FIRST EXECUTABLE STATEMENT FDUMP RETURN END INTEGER FUNCTION I1MACH(I) C***BEGIN PROLOGUE I1MACH C***DATE WRITTEN 750101 (YYMMDD) C***REVISION DATE 910131 (YYMMDD) C***CATEGORY NO. R1 C***KEYWORDS MACHINE CONSTANTS C***AUTHOR FOX, P. A., (BELL LABS) C HALL, A. D., (BELL LABS) C SCHRYER, N. L., (BELL LABS) C***PURPOSE Returns integer machine dependent constants C***DESCRIPTION C C This is the CMLIB version of I1MACH, the integer machine C constants subroutine originally developed for the PORT library. C C I1MACH can be used to obtain machine-dependent parameters C for the local machine environment. It is a function C subroutine with one (input) argument, and can be called C as follows, for example C C K = I1MACH(I) C C where I=1,...,16. The (output) value of K above is C determined by the (input) value of I. The results for C various values of I are discussed below. C C I/O unit numbers. C I1MACH( 1) = the standard input unit. C I1MACH( 2) = the standard output unit. C I1MACH( 3) = the standard punch unit. C I1MACH( 4) = the standard error message unit. C C Words. C I1MACH( 5) = the number of bits per integer storage unit. C I1MACH( 6) = the number of characters per integer storage unit. C C Integers. C assume integers are represented in the S-digit, base-A form C C sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) ) C C where 0 .LE. X(I) .LT. A for I=0,...,S-1. C I1MACH( 7) = A, the base. C I1MACH( 8) = S, the number of base-A digits. C I1MACH( 9) = A**S - 1, the largest magnitude. C C Floating-Point Numbers. C Assume floating-point numbers are represented in the T-digit, C base-B form C sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) C C where 0 .LE. X(I) .LT. B for I=1,...,T, C 0 .LT. X(1), and EMIN .LE. E .LE. EMAX. C I1MACH(10) = B, the base. C C Single-Precision C I1MACH(11) = T, the number of base-B digits. C I1MACH(12) = EMIN, the smallest exponent E. C I1MACH(13) = EMAX, the largest exponent E. C C Double-Precision C I1MACH(14) = T, the number of base-B digits. C I1MACH(15) = EMIN, the smallest exponent E. C I1MACH(16) = EMAX, the largest exponent E. C C To alter this function for a particular environment, C the desired set of DATA statements should be activated by C removing the C from column 1. Also, the values of C I1MACH(1) - I1MACH(4) should be checked for consistency C with the local operating system. C***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A C PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. C***ROUTINES CALLED (NONE) C***END PROLOGUE I1MACH C INTEGER IMACH(16),OUTPUT EQUIVALENCE (IMACH(4),OUTPUT) C C MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T C 3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T C PC 7300), AND 8087 BASED MICROS (E.G. IBM PC AND AT&T 6300). C C === MACHINE = IEEE.MOST-SIG-BYTE-FIRST C === MACHINE = IEEE.LEAST-SIG-BYTE-FIRST C === MACHINE = SUN C === MACHINE = 68000 C === MACHINE = 8087 C === MACHINE = IBM.PC C === MACHINE = ATT.3B C === MACHINE = ATT.7300 C === MACHINE = ATT.6300 DATA IMACH( 1) / 5 / DATA IMACH( 2) / 6 / DATA IMACH( 3) / 7 / DATA IMACH( 4) / 6 / DATA IMACH( 5) / 32 / DATA IMACH( 6) / 4 / DATA IMACH( 7) / 2 / DATA IMACH( 8) / 31 / DATA IMACH( 9) / 2147483647 / DATA IMACH(10) / 2 / DATA IMACH(11) / 24 / DATA IMACH(12) / -125 / DATA IMACH(13) / 128 / DATA IMACH(14) / 53 / DATA IMACH(15) / -1021 / DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C === MACHINE = AMDAHL C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C === MACHINE = BURROUGHS.1700 C DATA IMACH( 1) / 7 / C DATA IMACH( 2) / 2 / C DATA IMACH( 3) / 2 / C DATA IMACH( 4) / 2 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 33 / C DATA IMACH( 9) / Z1FFFFFFFF / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -256 / C DATA IMACH(13) / 255 / C DATA IMACH(14) / 60 / C DATA IMACH(15) / -256 / C DATA IMACH(16) / 255 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C === MACHINE = BURROUGHS.5700 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -50 / C DATA IMACH(16) / 76 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C === MACHINE = BURROUGHS.6700 C === MACHINE = BURROUGHS.7700 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 48 / C DATA IMACH( 6) / 6 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 39 / C DATA IMACH( 9) / O0007777777777777 / C DATA IMACH(10) / 8 / C DATA IMACH(11) / 13 / C DATA IMACH(12) / -50 / C DATA IMACH(13) / 76 / C DATA IMACH(14) / 26 / C DATA IMACH(15) / -32754 / C DATA IMACH(16) / 32780 / C C MACHINE CONSTANTS FOR THE CONVEX C-120 (NATIVE MODE) C C === MACHINE = CONVEX.C1 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1023 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE CONVEX (NATIVE MODE) C WITH -R8 OPTION C C === MACHINE = CONVEX.C1.R8 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 53 / C DATA IMACH(12) / -1023 / C DATA IMACH(13) / 1023 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1023 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE CONVEX C-120 (IEEE MODE) C C === MACHINE = CONVEX.C1.IEEE C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE CONVEX (IEEE MODE) C WITH -R8 OPTION C C === MACHINE = CONVEX.C1.IEEE.R8 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 53 / C DATA IMACH(12) / -1021 / C DATA IMACH(13) / 1024 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE CYBER 170/180 SERIES USING NOS (FTN5). C C === MACHINE = CYBER.170.NOS C === MACHINE = CYBER.180.NOS C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 60 / C DATA IMACH( 6) / 10 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 48 / C DATA IMACH( 9) / O"00007777777777777777" / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 48 / C DATA IMACH(12) / -974 / C DATA IMACH(13) / 1070 / C DATA IMACH(14) / 96 / C DATA IMACH(15) / -927 / C DATA IMACH(16) / 1070 / C C MACHINE CONSTANTS FOR THE CDC 180 SERIES USING NOS/VE C C === MACHINE = CYBER.180.NOS/VE C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / 9223372036854775807 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -4095 / C DATA IMACH(13) / 4094 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -4095 / C DATA IMACH(16) / 4094 / C C MACHINE CONSTANTS FOR THE CYBER 205 C C === MACHINE = CYBER.205 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 47 / C DATA IMACH( 9) / X'00007FFFFFFFFFFF' / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -28625 / C DATA IMACH(13) / 28718 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -28625 / C DATA IMACH(16) / 28718 / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES. C C === MACHINE = CDC.6000 C === MACHINE = CDC.7000 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 60 / C DATA IMACH( 6) / 10 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 48 / C DATA IMACH( 9) / 00007777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 48 / C DATA IMACH(12) / -974 / C DATA IMACH(13) / 1070 / C DATA IMACH(14) / 96 / C DATA IMACH(15) / -927 / C DATA IMACH(16) / 1070 / C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C USING THE 46 BIT INTEGER COMPILER OPTION C C === MACHINE = CRAY.46-BIT-INTEGER C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 102 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 46 / C DATA IMACH( 9) / 777777777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -8189 / C DATA IMACH(13) / 8190 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -8099 / C DATA IMACH(16) / 8190 / C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C USING THE 64 BIT INTEGER COMPILER OPTION C C === MACHINE = CRAY.64-BIT-INTEGER C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 102 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 64 / C DATA IMACH( 6) / 8 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 63 / C DATA IMACH( 9) / 777777777777777777777B / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 47 / C DATA IMACH(12) / -8189 / C DATA IMACH(13) / 8190 / C DATA IMACH(14) / 94 / C DATA IMACH(15) / -8099 / C DATA IMACH(16) / 8190 /C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C === MACHINE = DATA_GENERAL.ECLIPSE.S/200 C DATA IMACH( 1) / 11 / C DATA IMACH( 2) / 12 / C DATA IMACH( 3) / 8 / C DATA IMACH( 4) / 10 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) /32767 / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C ELXSI 6400 C C === MACHINE = ELSXI.6400 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 32 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -126 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1022 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE HARRIS 220 C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7 C C === MACHINE = HARRIS.220 C === MACHINE = HARRIS.SLASH6 C === MACHINE = HARRIS.SLASH7 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 0 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 24 / C DATA IMACH( 6) / 3 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 23 / C DATA IMACH( 9) / 8388607 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 23 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 38 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES. C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. C C === MACHINE = HONEYWELL.600/6000 C === MACHINE = HONEYWELL.DPS.8/70 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 43 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 63 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 3 WORD DOUBLE PRECISION OPTION WITH FTN4 C C === MACHINE = HP.2100.3_WORD_DP C DATA IMACH(1) / 5/ C DATA IMACH(2) / 6 / C DATA IMACH(3) / 4 / C DATA IMACH(4) / 1 / C DATA IMACH(5) / 16 / C DATA IMACH(6) / 2 / C DATA IMACH(7) / 2 / C DATA IMACH(8) / 15 / C DATA IMACH(9) / 32767 / C DATA IMACH(10)/ 2 / C DATA IMACH(11)/ 23 / C DATA IMACH(12)/ -128 / C DATA IMACH(13)/ 127 / C DATA IMACH(14)/ 39 / C DATA IMACH(15)/ -128 / C DATA IMACH(16)/ 127 / C C MACHINE CONSTANTS FOR THE HP 2100 C 4 WORD DOUBLE PRECISION OPTION WITH FTN4 C C === MACHINE = HP.2100.4_WORD_DP C DATA IMACH(1) / 5 / C DATA IMACH(2) / 6 / C DATA IMACH(3) / 4 / C DATA IMACH(4) / 1 / C DATA IMACH(5) / 16 / C DATA IMACH(6) / 2 / C DATA IMACH(7) / 2 / C DATA IMACH(8) / 15 / C DATA IMACH(9) / 32767 / C DATA IMACH(10)/ 2 / C DATA IMACH(11)/ 23 / C DATA IMACH(12)/ -128 / C DATA IMACH(13)/ 127 / C DATA IMACH(14)/ 55 / C DATA IMACH(15)/ -128 / C DATA IMACH(16)/ 127 / C C HP 9000 C C === MACHINE = HP.9000 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 7 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 32 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -126 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1015 / C DATA IMACH(16) / 1017 / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86 AND C THE INTERDATA 3230 AND INTERDATA 7/32. C C === MACHINE = IBM.360 C === MACHINE = IBM.370 C === MACHINE = XEROX.SIGMA.5 C === MACHINE = XEROX.SIGMA.7 C === MACHINE = XEROX.SIGMA.9 C === MACHINE = SEL.85 C === MACHINE = SEL.86 C === MACHINE = INTERDATA.3230 C === MACHINE = INTERDATA.7/32 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / Z7FFFFFFF / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 63 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 63 / C C MACHINE CONSTANTS FOR THE INTERDATA 8/32 C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. C C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. C C === MACHINE = INTERDATA.8/32.UNIX C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 6 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / Z'7FFFFFFF' / C DATA IMACH(10) / 16 / C DATA IMACH(11) / 6 / C DATA IMACH(12) / -64 / C DATA IMACH(13) / 62 / C DATA IMACH(14) / 14 / C DATA IMACH(15) / -64 / C DATA IMACH(16) / 62 / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C === MACHINE = PDP-10.KA C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 54 / C DATA IMACH(15) / -101 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C === MACHINE = PDP-10.KI C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 5 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / "377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 62 / C DATA IMACH(15) / -128 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 32-BIT INTEGER ARITHMETIC. C C === MACHINE = PDP-11.32-BIT C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 16-BIT INTEGER ARITHMETIC. C C === MACHINE = PDP-11.16-BIT C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 16 / C DATA IMACH( 6) / 2 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 15 / C DATA IMACH( 9) / 32767 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -127 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 56 / C DATA IMACH(15) / -127 / C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000. C C === MACHINE = SEQUENT.BALANCE.8000 C DATA IMACH( 1) / 0 / C DATA IMACH( 2) / 0 / C DATA IMACH( 3) / 7 / C DATA IMACH( 4) / 0 / C DATA IMACH( 5) / 32 / C DATA IMACH( 6) / 1 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 31 / C DATA IMACH( 9) / 2147483647 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 24 / C DATA IMACH(12) / -125 / C DATA IMACH(13) / 128 / C DATA IMACH(14) / 53 / C DATA IMACH(15) / -1021 / C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER C C === MACHINE = UNIVAC.1100 C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / C DATA IMACH( 3) / 1 / C DATA IMACH( 4) / 6 / C DATA IMACH( 5) / 36 / C DATA IMACH( 6) / 4 / C DATA IMACH( 7) / 2 / C DATA IMACH( 8) / 35 / C DATA IMACH( 9) / O377777777777 / C DATA IMACH(10) / 2 / C DATA IMACH(11) / 27 / C DATA IMACH(12) / -128 / C DATA IMACH(13) / 127 / C DATA IMACH(14) / 60 / C DATA IMACH(15) /-1024 / C DATA IMACH(16) / 1023 / C C MACHINE CONSTANTS FOR THE VAX 11/780 C C === MACHINE = VAX.11/780 C DATA IMACH(1) / 5 / C DATA IMACH(2) / 6 / C DATA IMACH(3) / 5 / C DATA IMACH(4) / 6 / C DATA IMACH(5) / 32 / C DATA IMACH(6) / 4 / C DATA IMACH(7) / 2 / C DATA IMACH(8) / 31 / C DATA IMACH(9) /2147483647 / C DATA IMACH(10)/ 2 / C DATA IMACH(11)/ 24 / C DATA IMACH(12)/ -127 / C DATA IMACH(13)/ 127 / C DATA IMACH(14)/ 56 / C DATA IMACH(15)/ -127 / C DATA IMACH(16)/ 127 / C C C***FIRST EXECUTABLE STATEMENT I1MACH IF (I .LT. 1 .OR. I .GT. 16) 1 CALL XERROR ( 'I1MACH -- I OUT OF BOUNDS',25,1,2) C I1MACH=IMACH(I) RETURN C END FUNCTION J4SAVE(IWHICH,IVALUE,ISET) C***BEGIN PROLOGUE J4SAVE C***REFER TO XERROR C Abstract C J4SAVE saves and recalls several global variables needed C by the library error handling routines. C C Description of Parameters C --Input-- C IWHICH - Index of item desired. C = 1 Refers to current error number. C = 2 Refers to current error control flag. C = 3 Refers to current unit number to which error C messages are to be sent. (0 means use standard.) C = 4 Refers to the maximum number of times any C message is to be printed (as set by XERMAX). C = 5 Refers to the total number of units to which C each error message is to be written. C = 6 Refers to the 2nd unit for error messages C = 7 Refers to the 3rd unit for error messages C = 8 Refers to the 4th unit for error messages C = 9 Refers to the 5th unit for error messages C IVALUE - The value to be set for the IWHICH-th parameter, C if ISET is .TRUE. . C ISET - If ISET=.TRUE., the IWHICH-th parameter will BE C given the value, IVALUE. If ISET=.FALSE., the C IWHICH-th parameter will be unchanged, and IVALUE C is a dummy parameter. C --Output-- C The (old) value of the IWHICH-th parameter will be returned C in the function value, J4SAVE. C C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C Adapted from Bell Laboratories PORT Library Error Handler C Latest revision --- 23 MAY 1979 C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED (NONE) C***END PROLOGUE J4SAVE LOGICAL ISET INTEGER IPARAM(9) SAVE IPARAM DATA IPARAM(1),IPARAM(2),IPARAM(3),IPARAM(4)/0,2,0,10/ DATA IPARAM(5)/1/ DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/ C***FIRST EXECUTABLE STATEMENT J4SAVE J4SAVE = IPARAM(IWHICH) IF (ISET) IPARAM(IWHICH) = IVALUE RETURN END SUBROUTINE XERABT(MESSG,NMESSG) C***BEGIN PROLOGUE XERABT C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Aborts program execution and prints error message. C***DESCRIPTION C Abstract C ***Note*** machine dependent routine C XERABT aborts the execution of the program. C The error message causing the abort is given in the calling C sequence, in case one needs it for printing on a dayfile, C for example. C C Description of Parameters C MESSG and NMESSG are as in XERROR, except that NMESSG may C be zero, in which case no message is being supplied. C C Written by Ron Jones, with SLATEC Common Math Library Subcommittee C Latest revision --- 19 MAR 1980 C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED (NONE) C***END PROLOGUE XERABT CHARACTER*(*) MESSG C***FIRST EXECUTABLE STATEMENT XERABT STOP END SUBROUTINE XERCTL(MESSG1,NMESSG,NERR,LEVEL,KONTRL) C***BEGIN PROLOGUE XERCTL C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. R3C C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Allows user control over handling of individual errors. C***DESCRIPTION C Abstract C Allows user control over handling of individual errors. C Just after each message is recorded, but before it is C processed any further (i.e., before it is printed or C a decision to abort is made), a call is made to XERCTL. C If the user has provided his own version of XERCTL, he C can then override the value of KONTROL used in processing C this message by redefining its value. C KONTRL may be set to any value from -2 to 2. C The meanings for KONTRL are the same as in XSETF, except C that the value of KONTRL changes only for this message. C If KONTRL is set to a value outside the range from -2 to 2, C it will be moved back into that range. C C Description of Parameters C C --Input-- C MESSG1 - the first word (only) of the error message. C NMESSG - same as in the call to XERROR or XERRWV. C NERR - same as in the call to XERROR or XERRWV. C LEVEL - same as in the call to XERROR or XERRWV. C KONTRL - the current value of the control flag as set C by a call to XSETF. C C --Output-- C KONTRL - the new value of KONTRL. If KONTRL is not C defined, it will remain at its original value. C This changed value of control affects only C the current occurrence of the current message. C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED (NONE) C***END PROLOGUE XERCTL CHARACTER*20 MESSG1 C***FIRST EXECUTABLE STATEMENT XERCTL RETURN END SUBROUTINE XERPRT(MESSG,NMESSG) C***BEGIN PROLOGUE XERPRT C***DATE WRITTEN 790801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. Z C***KEYWORDS ERROR,XERROR PACKAGE C***AUTHOR JONES, R. E., (SNLA) C***PURPOSE Prints error messages. C***DESCRIPTION C Abstract C Print the Hollerith message in MESSG, of length NMESSG, C on each file indicated by XGETUA. C Latest revision --- 19 MAR 1980 C***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- C HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, C 1982. C***ROUTINES CALLED I1MACH,S88FMT,XGETUA C***END PROLOGUE XERPRT INTEGER LUN(5) CHARACTER*(*) MESSG C OBTAIN UNIT NUMBERS AND WRITE LINE TO EACH UNIT C***FIRST EXECUTABLE STATEMENT XERPRT CALL XGETUA(LUN,NUNIT) LENMES = LEN(MESSG) DO 20 KUNIT=1,NUNIT IUNIT = LUN(KUNIT) IF (IUNIT.EQ.0) IUNIT = I1MACH(4) DO 10 ICHAR=1,LENMES,72 LAST = MIN0(ICHAR+71 , LENMES) WRITE (IUNIT,'(1X,A)') MESSG(ICHAR:LAST) 10 CONTINUE 20 CONTINUE RETURN END * ====================================================================== * NIST Guide to Available Math Software. * Fullsource for module GAMRAN from package DATAPAC. * Retrieved from CAMSUN on Fri Jul 10 08:30:32 1998. * ====================================================================== * GAMRAN SUBROUTINE GAMRAN(N,GAMMA,ISEED,X) C ******STILL NEEDS ALGORITHM WORK ****** C C PURPOSE--THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF SIZE N C FROM THE GAMMA DISTRIBUTION C WITH TAIL LENGTH PARAMETER VALUE = GAMMA. C THE PROTOTYPE GAMMA DISTRIBUTION USED C HEREIN HAS MEAN = GAMMA C AND STANDARD DEVIATION = SQRT(GAMMA). C THIS DISTRIBUTION IS DEFINED FOR ALL POSITIVE X, C AND HAS THE PROBABILITY DENSITY FUNCTION C F(X) = (1/CONSTANT) * (X**(GAMMA-1)) * EXP(-X) C WHERE THE CONSTANT = THE GAMMA FUNCTION EVALUATED C AT THE VALUE GAMMA. C INPUT ARGUMENTS--N = THE DESIRED INTEGER NUMBER C OF RANDOM NUMBERS TO BE C GENERATED. C --GAMMA = THE SINGLE PRECISION VALUE OF THE C TAIL LENGTH PARAMETER. C GAMMA SHOULD BE POSITIVE. C GAMMA SHOULD BE LARGER C THAN 1/3 (ALGORITHMIC RESTRICTION). C OUTPUT ARGUMENTS--X = A SINGLE PRECISION VECTOR C (OF DIMENSION AT LEAST N) C INTO WHICH THE GENERATED C RANDOM SAMPLE WILL BE PLACED. C OUTPUT--A RANDOM SAMPLE OF SIZE N C FROM THE GAMMA DISTRIBUTION C WITH TAIL LENGTH PARAMETER VALUE = GAMMA. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE C OF N FOR THIS SUBROUTINE. C --GAMMA SHOULD BE POSITIVE. C --GAMMA SHOULD BE LARGER C THAN 1/3 (ALGORITHMIC RESTRICTION). C OTHER DATAPAC SUBROUTINES NEEDED--UNIRAN, NORRAN. C FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, EXP. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN (1977) C REFERENCES--GREENWOOD, 'A FAST GENERATOR FOR C GAMMA-DISTRIBUTED RANDOM VARIABLES', C COMPSTAT 1974, PROCEEDINGS IN C COMPUTATIONAL STATISTICS, VIENNA, C SEPTEMBER, 1974, PAGES 19-27. C --TOCHER, THE ART OF SIMULATION, C 1963, PAGES 24-27. C --HAMMERSLEY AND HANDSCOMB, MONTE CARLO METHODS, C 1964, PAGES 36-37. C --WILK, GNANADESIKAN, AND HUYETT, 'PROBABILITY C PLOTS FOR THE GAMMA DISTRIBUTION', C TECHNOMETRICS, 1962, PAGES 1-15. C --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE C DISTRIBUTIONS--1, 1970, PAGES 166-206. C --HASTINGS AND PEACOCK, STATISTICAL C DISTRIBUTIONS--A HANDBOOK FOR C STUDENTS AND PRACTITIONERS, 1975, C PAGES 68-73. C --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS C SERIES 55, 1964, PAGE 952. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING DIVISION C CENTER FOR APPLIED MATHEMATICS C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE--301-921-3651 C NOTE--DATAPLOT IS A REGISTERED TRADEMARK C OF THE NATIONAL BUREAU OF STANDARDS. C THIS SUBROUTINE MAY NOT BE COPIED, EXTRACTED, C MODIFIED, OR OTHERWISE USED IN A CONTEXT C OUTSIDE OF THE DATAPLOT LANGUAGE/SYSTEM. C LANGUAGE--ANSI FORTRAN (1966) C EXCEPTION--HOLLERITH STRINGS IN FORMAT STATEMENTS C DENOTED BY QUOTES RATHER THAN NH. C VERSION NUMBER--82/7 C ORIGINAL VERSION--NOVEMBER 1975. C UPDATED --FEBRUARY 1976. C UPDATED --JUNE 1978. C UPDATED --DECEMBER 1981. C UPDATED --MARCH 1982. C UPDATED --MAY 1982. C C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES------------------- C C--------------------------------------------------------------------- C DIMENSION X(*) C C--------------------------------------------------------------------- C CCCCC CHARACTER*4 IFEEDB CCCCC CHARACTER*4 IPRINT C CCCCC COMMON /MACH/IRD,IPR,CPUMIN,CPUMAX,NUMBPC,NUMCPW,NUMBPW CCCCC COMMON /PRINT/IFEEDB,IPRINT C C-----DATA STATEMENTS------------------------------------------------- C DATA ATHIRD/0.3333333/ DATA SQRT3 /1.73205081/ C IPR=6 C C-----START POINT----------------------------------------------------- C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(N.LT.1)GOTO50 IF(GAMMA.LE.0.0)GOTO60 IF(GAMMA.LE.0.33333333)GOTO65 GOTO90 50 WRITE(IPR, 5) WRITE(IPR,47)N RETURN 60 WRITE(IPR,15) WRITE(IPR,46)GAMMA RETURN 65 WRITE(IPR,16) WRITE(IPR,17) WRITE(IPR,46)GAMMA RETURN 90 CONTINUE 5 FORMAT(1H , 91H***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO THE 1 GAMRAN SUBROUTINE IS NON-POSITIVE *****) 15 FORMAT(1H , 91H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 GAMRAN SUBROUTINE IS NON-POSITIVE *****) 16 FORMAT(1H ,114H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 GAMRAN SUBROUTINE IS SMALLER THAN OR EQUAL TO 0.33333333 *****) 17 FORMAT(1H , 44H (ALGORITHMIC RESTIRCTION)) 46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****) 47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8 ,6H *****) C C GENERATE N GAMMA DISTRIBUTION RANDOM NUMBERS C USING GREENWOOD'S REJECTION ALGORITHM-- C 1) GENERATE A NORMAL RANDOM NUMBER; C 2) TRANSFORM THE NORMAL VARIATE TO AN APPROXIMATE C GAMMA VARIATE USING THE WILSON-HILFERTY C APPROXIMATION (SEE THE JOHNSON AND KOTZ C REFERENCE, PAGE 176); C 3) FORM THE REJECTION FUNCTION VALUE, BASED C ON THE PROBABILITY DENSITY FUNCTION VALUE C OF THE ACTUAL DISTRIBUTION OF THE PSEUDO-GAMMA C VARIATE, AND THE PROBABILITY DENSITY FUNCTION VALUE C OF A TRUE GAMMA VARIATE. C 4) GENERATE A UNIFORM RANDOM NUMBER; C 5) IF THE UNIFORM RANDOM NUMBER IS LESS THAN C THE REJECTION FUNCTION VALUE, THEN ACCEPT C THE PSEUDO-RANDOM NUMBER AS A GAMMA VARIATE; C IF THE UNIFORM RANDOM NUMBER IS LARGER THAN C THE REJECTION FUNCTION VALUE, THEN REJECT C THE PSEUDO-RANDOM NUMBER AS A GAMMA VARIATE. C A1=1.0/(9.0*GAMMA) B1=SQRT(A1) XN0=-SQRT3+B1 XG0=GAMMA*(1.0-A1+B1*XN0)**3 DO100I=1,N 150 CALL NORRAN(1,ISEED,XN) XG=GAMMA*(1.0-A1+B1*XN)**3 IF(XG.LT.0.0)GOTO150 TERM=(XG/XG0)**(GAMMA-ATHIRD) ARG=0.5*XN*XN-XG-0.5*XN0*XN0+XG0 FUNCT=TERM*EXP(ARG) CALL UNIRAN(1,ISEED,U) IF(U.LE.FUNCT)GOTO170 GOTO150 170 X(I)=XG 100 CONTINUE C RETURN END * NORRAN SUBROUTINE NORRAN(N,ISEED,X) C C PURPOSE--THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF SIZE N C FROM THE THE NORMAL (GAUSSIAN) C DISTRIBUTION WITH MEAN = 0 AND STANDARD DEVIATION = 1. C THIS DISTRIBUTION IS DEFINED FOR ALL X AND HAS C THE PROBABILITY DENSITY FUNCTION C F(X) = (1/SQRT(2*PI))*EXP(-X*X/2). C INPUT ARGUMENTS--N = THE DESIRED INTEGER NUMBER C OF RANDOM NUMBERS TO BE C GENERATED. C OUTPUT ARGUMENTS--X = A SINGLE PRECISION VECTOR C (OF DIMENSION AT LEAST N) C INTO WHICH THE GENERATED C RANDOM SAMPLE WILL BE PLACED. C OUTPUT--A RANDOM SAMPLE OF SIZE N C FROM THE NORMAL DISTRIBUTION C WITH MEAN = 0 AND STANDARD DEVIATION = 1. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE C OF N FOR THIS SUBROUTINE. C OTHER DATAPAC SUBROUTINES NEEDED--UNIRAN. C FORTRAN LIBRARY SUBROUTINES NEEDED--ALOG, SQRT, SIN, COS. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN (1977) C METHOD--BOX-MULLER ALGORITHM. C REFERENCES--BOX AND MULLER, 'A NOTE ON THE GENERATION C OF RANDOM NORMAL DEVIATES', JOURNAL OF THE C ASSOCIATION FOR COMPUTING MACHINERY, 1958, C PAGES 610-611. C --TOCHER, THE ART OF SIMULATION, C 1963, PAGES 33-34. C --HAMMERSLEY AND HANDSCOMB, MONTE CARLO METHODS, C 1964, PAGE 39. C --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE C DISTRIBUTIONS--1, 1970, PAGES 40-111. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING DIVISION C CENTER FOR APPLIED MATHEMATICS C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE--301-921-3651 C NOTE--DATAPLOT IS A REGISTERED TRADEMARK C OF THE NATIONAL BUREAU OF STANDARDS. C THIS SUBROUTINE MAY NOT BE COPIED, EXTRACTED, C MODIFIED, OR OTHERWISE USED IN A CONTEXT C OUTSIDE OF THE DATAPLOT LANGUAGE/SYSTEM. C LANGUAGE--ANSI FORTRAN (1966) C EXCEPTION--HOLLERITH STRINGS IN FORMAT STATEMENTS C DENOTED BY QUOTES RATHER THAN NH. C VERSION NUMBER--82.6 C ORIGINAL VERSION--JUNE 1972. C UPDATED --SEPTEMBER 1975. C UPDATED --NOVEMBER 1975. C UPDATED --JULY 1976. C UPDATED --DECEMBER 1981. C UPDATED --MAY 1982. C C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES------------------- C C--------------------------------------------------------------------- C DIMENSION X(*) DIMENSION Y(2) C C--------------------------------------------------------------------- C CCCCC CHARACTER*4 IFEEDB CCCCC CHARACTER*4 IPRINT C CCCCC COMMON /MACH/IRD,IPR,CPUMIN,CPUMAX,NUMBPC,NUMCPW,NUMBPW CCCCC COMMON /PRINT/IFEEDB,IPRINT C IPR=6 C C-----DATA STATEMENTS------------------------------------------------- C DATA PI/3.14159265359/ C C-----START POINT----------------------------------------------------- C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(N.LT.1)GOTO50 GOTO90 50 WRITE(IPR, 5) WRITE(IPR,47)N RETURN 90 CONTINUE 5 FORMAT(1H , 91H***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO THE 1 NORRAN SUBROUTINE IS NON-POSITIVE *****) 47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8 ,6H *****) C C GENERATE N UNIFORM (0,1) RANDOM NUMBERS; C THEN GENERATE 2 ADDITIONAL UNIFORM (0,1) RANDOM NUMBERS C (TO BE USED BELOW IN FORMING THE N-TH NORMAL C RANDOM NUMBER WHEN THE DESIRED SAMPLE SIZE N C HAPPENS TO BE ODD). C CALL UNIRAN(N,ISEED,X) CALL UNIRAN(2,ISEED,Y) C C GENERATE N NORMAL RANDOM NUMBERS C USING THE BOX-MULLER METHOD. C DO200I=1,N,2 IP1=I+1 U1=X(I) IF(I.EQ.N)GOTO210 U2=X(IP1) GOTO220 210 U2=Y(2) 220 ARG1=-2.0*ALOG(U1) ARG2=2.0*PI*U2 SQRT1=SQRT(ARG1) Z1=SQRT1*COS(ARG2) Z2=SQRT1*SIN(ARG2) X(I)=Z1 IF(I.EQ.N)GOTO200 X(IP1)=Z2 200 CONTINUE C RETURN END * UNIRAN SUBROUTINE UNIRAN(N,ISEED,X) C C PURPOSE--THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF SIZE N C FROM THE UNIFORM (RECTANGULAR) C DISTRIBUTION ON THE UNIT INTERVAL (0,1). C THIS DISTRIBUTION HAS MEAN = 0.5 C AND STANDARD DEVIATION = SQRT(1/12) = 0.28867513. C THIS DISTRIBUTION HAS THE PROBABILITY C DENSITY FUNCTION F(X) = 1. C C INPUT ARGUMENTS--N = THE DESIRED INTEGER NUMBER C OF RANDOM NUMBERS TO BE C GENERATED. C --ISEED = AN INTEGER ISEED VALUE C OUTPUT ARGUMENTS--X = A SINGLE PRECISION VECTOR C (OF DIMENSION AT LEAST N) C INTO WHICH THE GENERATED C RANDOM SAMPLE WILL BE PLACED. C OUTPUT--A RANDOM SAMPLE OF SIZE N C FROM THE RECTANGULAR DISTRIBUTION ON (0,1). C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE C OF N FOR THIS SUBROUTINE. C OTHER SUBROUTINES NEEDED--NONE. C FORTRAN LIBRARY SUBROUTINES NEEDED--NONE. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN (1977) C C ALGORITHM--FIBONACCI GENERATOR C AS DEFINED BY GEORGE MARSAGLIA. C C NOTE--THIS GENERATOR IS TRANSPORTABLE. C IT IS NOT MACHINE-INDEPENDENT C IN THE SENSE THAT FOR A GIVEN VALUE C OF THE INPUT SEED ISEED AND FOR A GIVEN VALUE C OF MDIG (TO BE DEFINED BELOW), C THE SAME SEQUENCE OF UNIRFORM RANDOM C NUMBERS WILL RESULT ON DIFFERENT COMPUTERS C (VAX, PRIME, PERKIN-ELMER, IBM, UNIVAC, HONEYWELL, ETC.) C C NOTE--IF MDIG = 32 AND IF ISEED = 305, C THEN THE OUTPUT FROM THIS GENERATOR SHOULD BE AS FOLLOWS-- C THE FIRST NUMBER TO RESULT IS .4771580... C THE SECOND NUMBER TO RESULT IS .4219293... C THE THIRD NUMBER TO RESULT IS .6646181... C ... C THE THOUSANDTH NUMBER TO RESULT IS .2036834... C C NOTE--IF MDIG = 16 AND IF ISEED = 305, C THEN THE OUTPUT FROM THIS GENERATOR SHOULD BE AS FOLLOWS-- C THE FIRST NUMBER TO RESULT IS .027832881... C THE SECOND NUMBER TO RESULT IS .56102176... C THE THIRD NUMBER TO RESULT IS .41456343... C ... C THE THOUSANDTH NUMBER TO RESULT IS .19797357... C C NOTE--IT IS RECOMMENDED THAT UPON C IMPLEMENTATION OF DATAPLOT, THE OUTPUT C FROM UNIRAN BE CHECKED FOR AGREEMENT C WITH THE ABOVE SAMPLE OUTPUT. C ALSO, THERE ARE MANY ANALYSIS AND DIAGNOSTIC C TOOLS IN DATAPLOT THAT WILL ALLOW THE C TESTING OF THE RANDOMNESS AND UNIFORMITY C OF THIS GENERATOR. C SUCH CHECKING IS ESPECIALLY IMPORTANT C IN LIGHT OF THE FACT THAT OTHER DATAPLOT RANDOM C NUMBER GENERATOR SUBROUTINES (NORRAN--NORMAL, C LOGRAN--LOGISTIC, ETC.) ALL MAKE USE OF INTERMEDIATE C OUTPUT FROM UNIRAN. C C NOTE--THE OUTPUT FROM THIS SUBROUTINE DEPENDS C ON THE INPUT SEED (ISEED) AND ON THE C VALUE OF MDIG. C MDIG MAY NOT BE SMALLER THAN 16. C MDIG MAY NOT BE LARGER THAN MAX INTEGER ON YOUR COMPUTER. C C NOTE--BECAUSE OF THE PREPONDERANCE OF MAINFRAMES C WHICH HAVE WORDS OF 32 BITS AND LARGER C (E.G, VAX (= 32 BITS), UNIVAC (= 36 BITS), CDC (= 60 BITS), ETC.) C MDIG HAS BEEN SET TO 32. C THUS THE SAME SEQUENCE OF RANDOM NUMBERS SHOULD RESULT C ON ALL OF THESE COMPUTERS. C C NOTE--FOR SMALLER WORD SIZE COMPUTERS (E.G., 24-BIT AND 16-BIT), C THE VALUE OF MDIG SHOULD BE CHANGED TO 24 OR 16. C IN SUCH CASE, THE OUTPUT WILL NOT BE IDENTICAL TO C THE OUTPUT WHEN MDIG = 32. C C NOTE--THE CYCLE OF THE RANDOM NUMBERS DEPENDS ON MDIG. C THE CYCLE FROM MDIG = 32 IS LONG ENOUGH FOR MOST C PRACTICAL APPLICATIONS. C IF A LONGER CYCLE IS DESIRED, THEN INCREASE MDIG. C C NOTE--THE SEED MAY BE ANY POSITIVE INTEGER. C NO APPRECIABLE DIFFERENCE IN THE QUALITY C OF THE RANDOM NUMBERS HAS BEEN NOTED C BY THE CHOICE OF THE SEED. THERE IS NO C NEED TO USE PRIMES, NOR TO USE EXCEPTIONALLY C LARGE NUMBERS, ETC. C C REFERENCES--MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM C NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U. C --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE C DISTRIBUTIONS--2, 1970, PAGES 57-74. C WRITTEN BY--JAMES BLUE C SCIENTIFIC COMPUTING DIVISION C CENTER FOR APPLIED MATHEMATICS C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C --DAVID KAHANER C SCIENTIFIC COMPUTING DIVISION C CENTER FOR APPLIED MATHEMATICS C NATIONAL BUREAU OF STANDARDS C --GEORGE MARSAGLIA C COMPUTER SCIENCE DEPARTMENT C WASHINGTON STATE UNIVERSITY C --JAMES J. FILLIBEN C STATISTICAL ENGINEERING DIVISION C CENTER FOR APPLIED MATHEMATICS C NATIONAL BUREAU OF STANDARDS C C LANGUAGE--ANSI FORTRAN (1977) C ORIGINAL VERSION--JUNE 1972. C UPDATED --AUGUST 1974. C UPDATED --SEPTEM