c c c -------< ACM Chi-square code from Anderson Cancer Center in Texas: >------- c SUBROUTINE cdfchi(which,p,q,x,df,status,bound) c********************************************************************** c c SUBROUTINE CDFCHI( WHICH, P, Q, X, DF, STATUS, BOUND ) c Cumulative Distribution Function c CHI-Square distribution c c c Function c c c Calculates any one parameter of the chi-square c distribution given values for the others. c c c Arguments c c c WHICH --> Integer indicating which of the next three argument c values is to be calculated from the others. c Legal range: 1..3 c iwhich = 1 : Calculate P and Q from X and DF c iwhich = 2 : Calculate X from P,Q and DF c iwhich = 3 : Calculate DF from P,Q and X c INTEGER WHICH c c P <--> The integral from 0 to X of the chi-square c distribution. c Input range: [0, 1]. c DOUBLE PRECISION P c c Q <--> 1-P. c Input range: (0, 1]. c P + Q = 1.0. c DOUBLE PRECISION Q c c X <--> Upper limit of integration of the non-central c chi-square distribution. c Input range: [0, +infinity). c Search range: [0,1E100] c DOUBLE PRECISION X c c DF <--> Degrees of freedom of the c chi-square distribution. c Input range: (0, +infinity). c Search range: [ 1E-100, 1E100] c DOUBLE PRECISION DF c c STATUS <-- 0 if calculation completed correctly c -I if input parameter number I is out of range c 1 if answer appears to be lower than lowest c search bound c 2 if answer appears to be higher than greatest c search bound c 3 if P + Q .ne. 1 c 10 indicates error returned from cumgam. See c references in cdfgam c INTEGER STATUS c c BOUND <-- Undefined if STATUS is 0 c c Bound exceeded by parameter number I if STATUS c is negative. c c Lower search bound if STATUS is 1. c c Upper search bound if STATUS is 2. c c c Method c c c Formula 26.4.19 of Abramowitz and Stegun, Handbook of c Mathematical Functions (1966) is used to reduce the chisqure c distribution to the incomplete distribution. c c Computation of other parameters involve a seach for a value that c produces the desired value of P. The search relies on the c monotinicity of P with the other parameter. c c********************************************************************** c .. Parameters .. DOUBLE PRECISION tol PARAMETER (tol=1.0D-8) DOUBLE PRECISION atol PARAMETER (atol=1.0D-50) DOUBLE PRECISION zero,inf PARAMETER (zero=1.0D-100,inf=1.0D100) c .. c .. Scalar Arguments .. DOUBLE PRECISION bound,df,p,q,x INTEGER status,which c .. c .. Local Scalars .. DOUBLE PRECISION ccum,cum,fx,porq,pq LOGICAL qhi,qleft,qporq c .. c .. External Functions .. DOUBLE PRECISION spmpar EXTERNAL spmpar c .. c .. External Subroutines .. EXTERNAL cumchi,dinvr,dstinv c .. c .. Intrinsic Functions .. INTRINSIC abs c .. porq = 0.0D0 !I assume that this is the desired initial value - LEW IF (.NOT. ((which.LT.1).OR. (which.GT.3))) GO TO 30 IF (.NOT. (which.LT.1)) GO TO 10 bound = 1.0D0 GO TO 20 10 bound = 3.0D0 20 status = -1 RETURN 30 IF (which.EQ.1) GO TO 70 IF (.NOT. ((p.LT.0.0D0).OR. (p.GT.1.0D0))) GO TO 60 IF (.NOT. (p.LT.0.0D0)) GO TO 40 bound = 0.0D0 GO TO 50 40 bound = 1.0D0 50 status = -2 RETURN 60 CONTINUE 70 IF (which.EQ.1) GO TO 110 IF (.NOT. ((q.LE.0.0D0).OR. (q.GT.1.0D0))) GO TO 100 IF (.NOT. (q.LE.0.0D0)) GO TO 80 bound = 0.0D0 GO TO 90 80 bound = 1.0D0 90 status = -3 RETURN 100 CONTINUE 110 IF (which.EQ.2) GO TO 130 IF (.NOT. (x.LT.0.0D0)) GO TO 120 bound = 0.0D0 status = -4 RETURN 120 CONTINUE 130 IF (which.EQ.3) GO TO 150 IF (.NOT. (df.LE.0.0D0)) GO TO 140 bound = 0.0D0 status = -5 RETURN 140 CONTINUE 150 IF (which.EQ.1) GO TO 190 pq = p + q IF (.NOT. (abs(((pq)-0.5D0)-0.5D0).GT. + (3.0D0*spmpar(1)))) GO TO 180 IF (.NOT. (pq.LT.0.0D0)) GO TO 160 bound = 0.0D0 GO TO 170 160 bound = 1.0D0 170 status = 3 RETURN 180 CONTINUE 190 IF (which.EQ.1) GO TO 220 qporq = p .LE. q IF (.NOT. (qporq)) GO TO 200 porq = p GO TO 210 200 porq = q 210 CONTINUE 220 IF ((1).EQ. (which)) THEN status = 0 CALL cumchi(x,df,p,q) IF (porq.GT.1.5D0) THEN status = 10 RETURN END IF ELSE IF ((2).EQ. (which)) THEN x = 5.0D0 CALL dstinv(0.0D0,inf,0.5D0,0.5D0,5.0D0,atol,tol) status = 0 CALL dinvr(status,x,fx,qleft,qhi) 230 IF (.NOT. (status.EQ.1)) GO TO 270 CALL cumchi(x,df,cum,ccum) IF (.NOT. (qporq)) GO TO 240 fx = cum - p GO TO 250 240 fx = ccum - q 250 IF (.NOT. ((fx+porq).GT.1.5D0)) GO TO 260 status = 10 RETURN 260 CALL dinvr(status,x,fx,qleft,qhi) GO TO 230 270 IF (.NOT. (status.EQ.-1)) GO TO 300 IF (.NOT. (qleft)) GO TO 280 status = 1 bound = 0.0D0 GO TO 290 280 status = 2 bound = inf 290 CONTINUE 300 CONTINUE ELSE IF ((3).EQ. (which)) THEN df = 5.0D0 CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol) status = 0 CALL dinvr(status,df,fx,qleft,qhi) 310 IF (.NOT. (status.EQ.1)) GO TO 350 CALL cumchi(x,df,cum,ccum) IF (.NOT. (qporq)) GO TO 320 fx = cum - p GO TO 330 320 fx = ccum - q 330 IF (.NOT. ((fx+porq).GT.1.5D0)) GO TO 340 status = 10 RETURN 340 CALL dinvr(status,df,fx,qleft,qhi) GO TO 310 350 IF (.NOT. (status.EQ.-1)) GO TO 380 IF (.NOT. (qleft)) GO TO 360 status = 1 bound = zero GO TO 370 360 status = 2 bound = inf 370 CONTINUE 380 END IF RETURN END c c SUBROUTINE cumchi(x,df,cum,ccum) c********************************************************************** c c SUBROUTINE FUNCTION CUMCHI(X,DF,CUM,CCUM) c CUMulative of the CHi-square distribution c c c Function c c c Calculates the cumulative chi-square distribution. c c c Arguments c c c X --> Upper limit of integration of the c chi-square distribution. c X is DOUBLE PRECISION c c DF --> Degrees of freedom of the c chi-square distribution. c DF is DOUBLE PRECISION c c CUM <-- Cumulative chi-square distribution. c CUM is DOUBLE PRECISIO c c CCUM <-- Compliment of Cumulative chi-square distribution. c CCUM is DOUBLE PRECISI c c c Method c c c Calls incomplete gamma function (CUMGAM) c c********************************************************************** c .. Scalar Arguments .. DOUBLE PRECISION df,x,cum,ccum c .. c .. Local Scalars .. DOUBLE PRECISION a,xx c .. c .. External Subroutines .. EXTERNAL cumgam c .. c .. Executable Statements .. a = df*0.5D0 xx = x*0.5D0 CALL cumgam(xx,a,cum,ccum) RETURN END c c SUBROUTINE cumgam(x,a,cum,ccum) c********************************************************************** c c SUBROUTINE CUMGAM(X,A,CUM,CCUM) c Double precision cUMulative incomplete GAMma distribution c c c Function c c c Computes the cumulative of the incomplete gamma c distribution, i.e., the integral from 0 to X of c (1/GAM(A))*EXP(-T)*T**(A-1) DT c where GAM(A) is the complete gamma function of A, i.e., c GAM(A) = integral from 0 to infinity of c EXP(-T)*T**(A-1) DT c c c Arguments c c c X --> The upper limit of integration of the incomplete gamma. c X is DOUBLE PRECISION c c A --> The shape parameter of the incomplete gamma. c A is DOUBLE PRECISION c c CUM <-- Cumulative incomplete gamma distribution. c CUM is DOUBLE PRECISION c c CCUM <-- Compliment of Cumulative incomplete gamma distribution. c CCUM is DOUBLE PRECISIO c c c Method c c c Calls the routine GRATIO. c c********************************************************************** c c .. c .. Scalar Arguments .. DOUBLE PRECISION a,x,cum,ccum c .. c .. External Routines .. EXTERNAL gratio c .. c .. Executable Statements .. IF (.NOT. (x.LE.0.0D0)) GO TO 10 cum = 0.0D0 ccum = 1.0D0 RETURN 10 CALL gratio(a,x,cum,ccum,0) c Call gratio routine RETURN END c c SUBROUTINE dinvr(status,x,fx,qleft,qhi) c********************************************************************** c c SUBROUTINE DINVR(STATUS, X, FX, QLEFT, QHI) c Double precision c bounds the zero of the function and invokes zror c Reverse Communication c c c Function c c c Bounds the function and invokes ZROR to perform the zero c finding. STINVR must have been called before this routine c in order to set its parameters. c c c Arguments c c c STATUS <--> At the beginning of a zero finding problem, STATUS c should be set to 0 and INVR invoked. (The value c of parameters other than X will be ignored on this cal c c When INVR needs the function evaluated, it will set c STATUS to 1 and return. The value of the function c should be set in FX and INVR again called without c changing any of its other parameters. c c When INVR has finished without error, it will return c with STATUS 0. In that case X is approximately a root c of F(X). c c If INVR cannot bound the function, it returns status c -1 and sets QLEFT and QHI. c INTEGER STATUS c c X <-- The value of X at which F(X) is to be evaluated. c DOUBLE PRECISION X c c FX --> The value of F(X) calculated when INVR returns with c STATUS = 1. c DOUBLE PRECISION FX c c QLEFT <-- Defined only if QMFINV returns .FALSE. In that c case it is .TRUE. If the stepping search terminated c unsucessfully at SMALL. If it is .FALSE. the search c terminated unsucessfully at BIG. c QLEFT is LOGICAL c c QHI <-- Defined only if QMFINV returns .FALSE. In that c case it is .TRUE. if F(X) .GT. Y at the termination c of the search and .FALSE. if F(X) .LT. Y at the c termination of the search. c QHI is LOGICAL c c********************************************************************** c .. Scalar Arguments .. DOUBLE PRECISION fx,x,zabsst,zabsto,zbig,zrelst,zrelto,zsmall, + zstpmu INTEGER status LOGICAL qhi,qleft c .. c .. Local Scalars .. DOUBLE PRECISION absstp,abstol,big,fbig,fsmall,relstp,reltol, + small,step,stpmul,xhi,xlb,xlo,xsave,xub,yy,zx,zy, + zz INTEGER i99999 LOGICAL qbdd,qcond,qdum1,qdum2,qincr,qlim,qok,qup c .. c .. External Subroutines .. EXTERNAL dstzr,dzror c .. c .. Intrinsic Functions .. INTRINSIC abs,max,min c .. c .. Statement Functions .. LOGICAL qxmon c .. c .. Save statement .. SAVE c .. c .. Statement Function definitions .. qxmon(zx,zy,zz) = zx .LE. zy .AND. zy .LE. zz c .. c .. Executable Statements .. IF (status.GT.0) GO TO 310 qcond = .NOT. qxmon(small,x,big) IF (qcond) STOP ' SMALL, X, BIG not monotone in INVR' xsave = x c c See that SMALL and BIG bound the zero and set QINCR c x = small c GET-FUNCTION-VALUE ASSIGN 10 TO i99999 GO TO 300 10 fsmall = fx x = big c GET-FUNCTION-VALUE ASSIGN 20 TO i99999 GO TO 300 20 fbig = fx qincr = fbig .GT. fsmall IF (.NOT. (qincr)) GO TO 50 IF (fsmall.LE.0.0D0) GO TO 30 status = -1 qleft = .TRUE. qhi = .TRUE. RETURN 30 IF (fbig.GE.0.0D0) GO TO 40 status = -1 qleft = .FALSE. qhi = .FALSE. RETURN 40 GO TO 80 50 IF (fsmall.GE.0.0D0) GO TO 60 status = -1 qleft = .TRUE. qhi = .FALSE. RETURN 60 IF (fbig.LE.0.0D0) GO TO 70 status = -1 qleft = .FALSE. qhi = .TRUE. RETURN 70 CONTINUE 80 x = xsave step = max(absstp,relstp*abs(x)) c YY = F(X) - Y c GET-FUNCTION-VALUE ASSIGN 90 TO i99999 GO TO 300 90 yy = fx IF (.NOT. (yy.EQ.0.0D0)) GO TO 100 status = 0 qok = .TRUE. RETURN 100 qup = (qincr .AND. (yy.LT.0.0D0)) .OR. + (.NOT.qincr .AND. (yy.GT.0.0D0)) c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c HANDLE CASE IN WHICH WE MUST STEP HIGHER c c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IF (.NOT. (qup)) GO TO 170 xlb = xsave xub = min(xlb+step,big) GO TO 120 110 IF (qcond) GO TO 150 c YY = F(XUB) - Y 120 x = xub c GET-FUNCTION-VALUE ASSIGN 130 TO i99999 GO TO 300 130 yy = fx qbdd = (qincr .AND. (yy.GE.0.0D0)) .OR. + (.NOT.qincr .AND. (yy.LE.0.0D0)) qlim = xub .GE. big qcond = qbdd .OR. qlim IF (qcond) GO TO 140 step = stpmul*step xlb = xub xub = min(xlb+step,big) 140 GO TO 110 150 IF (.NOT. (qlim.AND..NOT.qbdd)) GO TO 160 status = -1 qleft = .FALSE. qhi = .NOT. qincr x = big RETURN 160 GO TO 240 c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c HANDLE CASE IN WHICH WE MUST STEP LOWER c c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 170 xub = xsave xlb = max(xub-step,small) GO TO 190 180 IF (qcond) GO TO 220 c YY = F(XLB) - Y 190 x = xlb c GET-FUNCTION-VALUE ASSIGN 200 TO i99999 GO TO 300 200 yy = fx qbdd = (qincr .AND. (yy.LE.0.0D0)) .OR. + (.NOT.qincr .AND. (yy.GE.0.0D0)) qlim = xlb .LE. small qcond = qbdd .OR. qlim IF (qcond) GO TO 210 step = stpmul*step xub = xlb xlb = max(xub-step,small) 210 GO TO 180 220 IF (.NOT. (qlim.AND..NOT.qbdd)) GO TO 230 status = -1 qleft = .TRUE. qhi = qincr x = small RETURN 230 CONTINUE 240 CALL dstzr(xlb,xub,abstol,reltol) c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c IF WE REACH HERE, XLB AND XUB BOUND THE ZERO OF F. c c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ status = 0 GO TO 260 250 IF (.NOT. (status.EQ.1)) GO TO 290 260 CALL dzror(status,x,fx,xlo,xhi,qdum1,qdum2) IF (.NOT. (status.EQ.1)) GO TO 280 c GET-FUNCTION-VALUE ASSIGN 270 TO i99999 GO TO 300 270 CONTINUE 280 GO TO 250 290 x = xlo status = 0 RETURN ENTRY dstinv(zsmall,zbig,zabsst,zrelst,zstpmu,zabsto,zrelto) c********************************************************************** c c SUBROUTINE DSTINV( SMALL, BIG, ABSSTP, RELSTP, STPMUL, c + ABSTOL, RELTOL ) c Double Precision - SeT INverse finder - Reverse Communication c c c Function c c c Concise Description - Given a monotone function F finds X c such that F(X) = Y. Uses Reverse communication -- see invr. c This routine sets quantities needed by INVR. c c More Precise Description of INVR - c c F must be a monotone function, the results of QMFINV are c otherwise undefined. QINCR must be .TRUE. if F is non- c decreasing and .FALSE. if F is non-increasing. c c QMFINV will return .TRUE. if and only if F(SMALL) and c F(BIG) bracket Y, i. e., c QINCR is .TRUE. and F(SMALL).LE.Y.LE.F(BIG) or c QINCR is .FALSE. and F(BIG).LE.Y.LE.F(SMALL) c c if QMFINV returns .TRUE., then the X returned satisfies c the following condition. let c TOL(X) = MAX(ABSTOL,RELTOL*ABS(X)) c then if QINCR is .TRUE., c F(X-TOL(X)) .LE. Y .LE. F(X+TOL(X)) c and if QINCR is .FALSE. c F(X-TOL(X)) .GE. Y .GE. F(X+TOL(X)) c c c Arguments c c c SMALL --> The left endpoint of the interval to be c searched for a solution. c SMALL is DOUBLE PRECISION c c BIG --> The right endpoint of the interval to be c searched for a solution. c BIG is DOUBLE PRECISION c c ABSSTP, RELSTP --> The initial step size in the search c is MAX(ABSSTP,RELSTP*ABS(X)). See algorithm. c ABSSTP is DOUBLE PRECISION c RELSTP is DOUBLE PRECISION c c STPMUL --> When a step doesn't bound the zero, the step c size is multiplied by STPMUL and another step c taken. A popular value is 2.0 c DOUBLE PRECISION STPMUL c c ABSTOL, RELTOL --> Two numbers that determine the accuracy c of the solution. See function for a precise definition. c ABSTOL is DOUBLE PRECISION c RELTOL is DOUBLE PRECISION c c c Method c c c Compares F(X) with Y for the input value of X then uses QINCR c to determine whether to step left or right to bound the c desired x. the initial step size is c MAX(ABSSTP,RELSTP*ABS(S)) for the input value of X. c Iteratively steps right or left until it bounds X. c At each step which doesn't bound X, the step size is doubled. c The routine is careful never to step beyond SMALL or BIG. If c it hasn't bounded X at SMALL or BIG, QMFINV returns .FALSE. c after setting QLEFT and QHI. c c If X is successfully bounded then Algorithm R of the paper c 'Two Efficient Algorithms with Guaranteed Convergence for c Finding a Zero of a Function' by J. C. P. Bus and c T. J. Dekker in ACM Transactions on Mathematical c Software, Volume 1, No. 4 page 330 (DEC. '75) is employed c to find the zero of the function F(X)-Y. This is routine c QRZERO. c c********************************************************************** small = zsmall big = zbig absstp = zabsst relstp = zrelst stpmul = zstpmu abstol = zabsto reltol = zrelto RETURN cc STOP '*** EXECUTION FLOWING INTO FLECS PROCEDURES ***' c TO GET-FUNCTION-VALUE 300 status = 1 RETURN 310 CONTINUE GO TO i99999 END c c SUBROUTINE dzror(status,x,fx,xlo,xhi,qleft,qhi) c********************************************************************** c c SUBROUTINE DZROR(STATUS, X, FX, XLO, XHI, QLEFT, QHI) c Double precision ZeRo of a function -- Reverse Communication c c c Function c c c Performs the zero finding. STZROR must have been called before c this routine in order to set its parameters. c c c Arguments c c c STATUS <--> At the beginning of a zero finding problem, STATUS c should be set to 0 and ZROR invoked. (The value c of other parameters will be ignored on this call.) c c When ZROR needs the function evaluated, it will set c STATUS to 1 and return. The value of the function c should be set in FX and ZROR again called without c changing any of its other parameters. c c When ZROR has finished without error, it will return c with STATUS 0. In that case (XLO,XHI) bound the answe c c If ZROR finds an error (which implies that F(XLO)-Y an c F(XHI)-Y have the same sign, it returns STATUS -1. In c this case, XLO and XHI are undefined. c INTEGER STATUS c c X <-- The value of X at which F(X) is to be evaluated. c DOUBLE PRECISION X c c FX --> The value of F(X) calculated when ZROR returns with c STATUS = 1. c DOUBLE PRECISION FX c c XLO <-- When ZROR returns with STATUS = 0, XLO bounds the c inverval in X containing the solution below. c DOUBLE PRECISION XLO c c XHI <-- When ZROR returns with STATUS = 0, XHI bounds the c inverval in X containing the solution above. c DOUBLE PRECISION XHI c c QLEFT <-- .TRUE. if the stepping search terminated unsucessfully c at XLO. If it is .FALSE. the search terminated c unsucessfully at XHI. c QLEFT is LOGICAL c c QHI <-- .TRUE. if F(X) .GT. Y at the termination of the c search and .FALSE. if F(X) .LT. Y at the c termination of the search. c QHI is LOGICAL c c********************************************************************** c .. Scalar Arguments .. DOUBLE PRECISION fx,x,xhi,xlo,zabstl,zreltl,zxhi,zxlo INTEGER status LOGICAL qhi,qleft c .. c .. Save statement .. SAVE c .. c .. Local Scalars .. DOUBLE PRECISION a,abstol,b,c,d,fa,fb,fc,fd,fda,fdb,m,mb,p,q, + reltol,tol,w,xxhi,xxlo,zx INTEGER ext,i99999 LOGICAL first,qrzero c .. c .. Intrinsic Functions .. INTRINSIC abs,max,sign c .. c .. Statement Functions .. DOUBLE PRECISION ftol c .. c .. Statement Function definitions .. ftol(zx) = 0.5D0*max(abstol,reltol*abs(zx)) c .. c .. Executable Statements .. IF (status.GT.0) GO TO 280 xlo = xxlo xhi = xxhi b = xlo x = xlo c GET-FUNCTION-VALUE ASSIGN 10 TO i99999 GO TO 270 10 fb = fx xlo = xhi a = xlo x = xlo c GET-FUNCTION-VALUE ASSIGN 20 TO i99999 GO TO 270 c c Check that F(ZXLO) < 0 < F(ZXHI) or c F(ZXLO) > 0 > F(ZXHI) c 20 IF (.NOT. (fb.LT.0.0D0)) GO TO 40 IF (.NOT. (fx.LT.0.0D0)) GO TO 30 status = -1 qleft = fx .LT. fb qhi = .FALSE. RETURN 30 CONTINUE 40 IF (.NOT. (fb.GT.0.0D0)) GO TO 60 IF (.NOT. (fx.GT.0.0D0)) GO TO 50 status = -1 qleft = fx .GT. fb qhi = .TRUE. RETURN 50 CONTINUE 60 fa = fx c first = .TRUE. 70 c = a fc = fa ext = 0 80 IF (.NOT. (abs(fc).LT.abs(fb))) GO TO 100 IF (.NOT. (c.NE.a)) GO TO 90 d = a fd = fa 90 a = b fa = fb xlo = c b = xlo fb = fc c = a fc = fa 100 tol = ftol(xlo) m = (c+b)*.5D0 mb = m - b IF (.NOT. (abs(mb).GT.tol)) GO TO 240 IF (.NOT. (ext.GT.3)) GO TO 110 w = mb GO TO 190 110 tol = sign(tol,mb) p = (b-a)*fb IF (.NOT. (first)) GO TO 120 q = fa - fb first = .FALSE. GO TO 130 120 fdb = (fd-fb)/ (d-b) fda = (fd-fa)/ (d-a) p = fda*p q = fdb*fa - fda*fb 130 IF (.NOT. (p.LT.0.0D0)) GO TO 140 p = -p q = -q 140 IF (ext.EQ.3) p = p*2.0D0 IF (.NOT. ((p*1.0D0).EQ.0.0D0.OR.p.LE. (q*tol))) GO TO 150 w = tol GO TO 180 150 IF (.NOT. (p.LT. (mb*q))) GO TO 160 w = p/q GO TO 170 160 w = mb 170 CONTINUE 180 CONTINUE 190 d = a fd = fa a = b fa = fb b = b + w xlo = b x = xlo c GET-FUNCTION-VALUE ASSIGN 200 TO i99999 GO TO 270 200 fb = fx IF (.NOT. ((fc*fb).GE.0.0D0)) GO TO 210 GO TO 70 210 IF (.NOT. (w.EQ.mb)) GO TO 220 ext = 0 GO TO 230 220 ext = ext + 1 230 GO TO 80 240 xhi = c qrzero = (fc.GE.0.0D0 .AND. fb.LE.0.0D0) .OR. + (fc.LT.0.0D0 .AND. fb.GE.0.0D0) IF (.NOT. (qrzero)) GO TO 250 status = 0 GO TO 260 250 status = -1 260 RETURN ENTRY dstzr(zxlo,zxhi,zabstl,zreltl) c********************************************************************** c c SUBROUTINE DSTZR( XLO, XHI, ABSTOL, RELTOL ) c Double precision SeT ZeRo finder - Reverse communication version c c c Function c c c c Sets quantities needed by ZROR. The function of ZROR c and the quantities set is given here. c c Concise Description - Given a function F c find XLO such that F(XLO) = 0. c c More Precise Description - c c Input condition. F is a double precision function of a single c double precision argument and XLO and XHI are such that c F(XLO)*F(XHI) .LE. 0.0 c c If the input condition is met, QRZERO returns .TRUE. c and output values of XLO and XHI satisfy the following c F(XLO)*F(XHI) .LE. 0. c ABS(F(XLO) .LE. ABS(F(XHI) c ABS(XLO-XHI) .LE. TOL(X) c where c TOL(X) = MAX(ABSTOL,RELTOL*ABS(X)) c c If this algorithm does not find XLO and XHI satisfying c these conditions then QRZERO returns .FALSE. This c implies that the input condition was not met. c c c Arguments c c c XLO --> The left endpoint of the interval to be c searched for a solution. c XLO is DOUBLE PRECISION c c XHI --> The right endpoint of the interval to be c for a solution. c XHI is DOUBLE PRECISION c c ABSTOL, RELTOL --> Two numbers that determine the accuracy c of the solution. See function for a c precise definition. c ABSTOL is DOUBLE PRECISION c RELTOL is DOUBLE PRECISION c c c Method c c c Algorithm R of the paper 'Two Efficient Algorithms with c Guaranteed Convergence for Finding a Zero of a Function' c by J. C. P. Bus and T. J. Dekker in ACM Transactions on c Mathematical Software, Volume 1, no. 4 page 330 c (Dec. '75) is employed to find the zero of F(X)-Y. c c********************************************************************** xxlo = zxlo xxhi = zxhi abstol = zabstl reltol = zreltl RETURN cc STOP '*** EXECUTION FLOWING INTO FLECS PROCEDURES ***' c TO GET-FUNCTION-VALUE 270 status = 1 RETURN 280 CONTINUE GO TO i99999 END c c DOUBLE PRECISION FUNCTION erf(x) c----------------------------------------------------------------------- c EVALUATION OF THE REAL ERROR FUNCTION c----------------------------------------------------------------------- c .. Scalar Arguments .. DOUBLE PRECISION x c .. c .. Local Scalars .. DOUBLE PRECISION ax,bot,c,t,top,x2 c .. c .. Local Arrays .. DOUBLE PRECISION a(5),b(3),p(8),q(8),r(5),s(4) c .. c .. Intrinsic Functions .. INTRINSIC abs,exp,sign c .. c .. Data statements .. c------------------------- c------------------------- c------------------------- c------------------------- DATA c/.564189583547756D0/ DATA a(1)/.771058495001320D-04/,a(2)/-.133733772997339D-02/, + a(3)/.323076579225834D-01/,a(4)/.479137145607681D-01/, + a(5)/.128379167095513D+00/ DATA b(1)/.301048631703895D-02/,b(2)/.538971687740286D-01/, + b(3)/.375795757275549D+00/ DATA p(1)/-1.36864857382717D-07/,p(2)/5.64195517478974D-01/, + p(3)/7.21175825088309D+00/,p(4)/4.31622272220567D+01/, + p(5)/1.52989285046940D+02/,p(6)/3.39320816734344D+02/, + p(7)/4.51918953711873D+02/,p(8)/3.00459261020162D+02/ DATA q(1)/1.00000000000000D+00/,q(2)/1.27827273196294D+01/, + q(3)/7.70001529352295D+01/,q(4)/2.77585444743988D+02/, + q(5)/6.38980264465631D+02/,q(6)/9.31354094850610D+02/, + q(7)/7.90950925327898D+02/,q(8)/3.00459260956983D+02/ DATA r(1)/2.10144126479064D+00/,r(2)/2.62370141675169D+01/, + r(3)/2.13688200555087D+01/,r(4)/4.65807828718470D+00/, + r(5)/2.82094791773523D-01/ DATA s(1)/9.41537750555460D+01/,s(2)/1.87114811799590D+02/, + s(3)/9.90191814623914D+01/,s(4)/1.80124575948747D+01/ c .. c .. Executable Statements .. c------------------------- ax = abs(x) IF (ax.GT.0.5D0) GO TO 10 t = x*x top = ((((a(1)*t+a(2))*t+a(3))*t+a(4))*t+a(5)) + 1.0D0 bot = ((b(1)*t+b(2))*t+b(3))*t + 1.0D0 erf = x* (top/bot) RETURN c 10 IF (ax.GT.4.0D0) GO TO 20 top = ((((((p(1)*ax+p(2))*ax+p(3))*ax+p(4))*ax+p(5))*ax+p(6))*ax+ + p(7))*ax + p(8) bot = ((((((q(1)*ax+q(2))*ax+q(3))*ax+q(4))*ax+q(5))*ax+q(6))*ax+ + q(7))*ax + q(8) erf = 0.5D0 + (0.5D0-exp(-x*x)*top/bot) IF (x.LT.0.0D0) erf = -erf RETURN c 20 IF (ax.GE.5.8D0) GO TO 30 x2 = x*x t = 1.0D0/x2 top = (((r(1)*t+r(2))*t+r(3))*t+r(4))*t + r(5) bot = (((s(1)*t+s(2))*t+s(3))*t+s(4))*t + 1.0D0 erf = (c-top/ (x2*bot))/ax erf = 0.5D0 + (0.5D0-exp(-x2)*erf) IF (x.LT.0.0D0) erf = -erf RETURN c 30 erf = sign(1.0D0,x) RETURN END c c DOUBLE PRECISION FUNCTION erfc1(ind,x) c----------------------------------------------------------------------- c EVALUATION OF THE COMPLEMENTARY ERROR FUNCTION c c ERFC1(IND,X) = ERFC(X) IF IND = 0 c ERFC1(IND,X) = EXP(X*X)*ERFC(X) OTHERWISE c----------------------------------------------------------------------- c .. Scalar Arguments .. DOUBLE PRECISION x INTEGER ind c .. c .. Local Scalars .. DOUBLE PRECISION ax,bot,c,e,t,top,w c .. c .. Local Arrays .. DOUBLE PRECISION a(5),b(3),p(8),q(8),r(5),s(4) c .. c .. External Functions .. DOUBLE PRECISION exparg EXTERNAL exparg c .. c .. Intrinsic Functions .. INTRINSIC abs,dble,exp c .. c .. Data statements .. c------------------------- c------------------------- c------------------------- c------------------------- DATA c/.564189583547756D0/ DATA a(1)/.771058495001320D-04/,a(2)/-.133733772997339D-02/, + a(3)/.323076579225834D-01/,a(4)/.479137145607681D-01/, + a(5)/.128379167095513D+00/ DATA b(1)/.301048631703895D-02/,b(2)/.538971687740286D-01/, + b(3)/.375795757275549D+00/ DATA p(1)/-1.36864857382717D-07/,p(2)/5.64195517478974D-01/, + p(3)/7.21175825088309D+00/,p(4)/4.31622272220567D+01/, + p(5)/1.52989285046940D+02/,p(6)/3.39320816734344D+02/, + p(7)/4.51918953711873D+02/,p(8)/3.00459261020162D+02/ DATA q(1)/1.00000000000000D+00/,q(2)/1.27827273196294D+01/, + q(3)/7.70001529352295D+01/,q(4)/2.77585444743988D+02/, + q(5)/6.38980264465631D+02/,q(6)/9.31354094850610D+02/, + q(7)/7.90950925327898D+02/,q(8)/3.00459260956983D+02/ DATA r(1)/2.10144126479064D+00/,r(2)/2.62370141675169D+01/, + r(3)/2.13688200555087D+01/,r(4)/4.65807828718470D+00/, + r(5)/2.82094791773523D-01/ DATA s(1)/9.41537750555460D+01/,s(2)/1.87114811799590D+02/, + s(3)/9.90191814623914D+01/,s(4)/1.80124575948747D+01/ c .. c .. Executable Statements .. c------------------------- c c ABS(X) .LE. 0.5 c ax = abs(x) IF (ax.GT.0.5D0) GO TO 10 t = x*x top = ((((a(1)*t+a(2))*t+a(3))*t+a(4))*t+a(5)) + 1.0D0 bot = ((b(1)*t+b(2))*t+b(3))*t + 1.0D0 erfc1 = 0.5D0 + (0.5D0-x* (top/bot)) IF (ind.NE.0) erfc1 = exp(t)*erfc1 RETURN c c 0.5 .LT. ABS(X) .LE. 4 c 10 IF (ax.GT.4.0D0) GO TO 20 top = ((((((p(1)*ax+p(2))*ax+p(3))*ax+p(4))*ax+p(5))*ax+p(6))*ax+ + p(7))*ax + p(8) bot = ((((((q(1)*ax+q(2))*ax+q(3))*ax+q(4))*ax+q(5))*ax+q(6))*ax+ + q(7))*ax + q(8) erfc1 = top/bot GO TO 40 c c ABS(X) .GT. 4 c 20 IF (x.LE.-5.6D0) GO TO 60 IF (ind.NE.0) GO TO 30 IF (x.GT.100.0D0) GO TO 70 IF (x*x.GT.-exparg(1)) GO TO 70 c 30 t = (1.0D0/x)**2 top = (((r(1)*t+r(2))*t+r(3))*t+r(4))*t + r(5) bot = (((s(1)*t+s(2))*t+s(3))*t+s(4))*t + 1.0D0 erfc1 = (c-t*top/bot)/ax c c FINAL ASSEMBLY c 40 IF (ind.EQ.0) GO TO 50 IF (x.LT.0.0D0) erfc1 = 2.0D0*exp(x*x) - erfc1 RETURN 50 w = dble(x)*dble(x) t = w e = w - dble(t) erfc1 = ((0.5D0+ (0.5D0-e))*exp(-t))*erfc1 IF (x.LT.0.0D0) erfc1 = 2.0D0 - erfc1 RETURN c c LIMIT VALUE FOR LARGE NEGATIVE X c 60 erfc1 = 2.0D0 IF (ind.NE.0) erfc1 = 2.0D0*exp(x*x) RETURN c c LIMIT VALUE FOR LARGE POSITIVE X c WHEN IND = 0 c 70 erfc1 = 0.0D0 RETURN END c c DOUBLE PRECISION FUNCTION exparg(l) c-------------------------------------------------------------------- c IF L = 0 THEN EXPARG(L) = THE LARGEST POSITIVE W FOR WHICH c EXP(W) CAN BE COMPUTED. c c IF L IS NONZERO THEN EXPARG(L) = THE LARGEST NEGATIVE W FOR c WHICH THE COMPUTED VALUE OF EXP(W) IS NONZERO. c c NOTE... ONLY AN APPROXIMATE VALUE FOR EXPARG(L) IS NEEDED. c-------------------------------------------------------------------- c .. Scalar Arguments .. INTEGER l c .. c .. Local Scalars .. DOUBLE PRECISION lnb INTEGER b,m c .. c .. External Functions .. INTEGER ipmpar EXTERNAL ipmpar c .. c .. Intrinsic Functions .. INTRINSIC dble,dlog c .. c .. Executable Statements .. c b = ipmpar(4) IF (b.NE.2) GO TO 10 lnb = .69314718055995D0 GO TO 40 10 IF (b.NE.8) GO TO 20 lnb = 2.0794415416798D0 GO TO 40 20 IF (b.NE.16) GO TO 30 lnb = 2.7725887222398D0 GO TO 40 30 lnb = dlog(dble(b)) c 40 IF (l.EQ.0) GO TO 50 m = ipmpar(9) - 1 exparg = 0.99999D0* (m*lnb) RETURN 50 m = ipmpar(10) exparg = 0.99999D0* (m*lnb) RETURN END c c DOUBLE PRECISION FUNCTION gam1(a) c ------------------------------------------------------------------ c COMPUTATION OF 1/GAMMA(A+1) - 1 FOR -0.5 .LE. A .LE. 1.5 c ------------------------------------------------------------------ c .. Scalar Arguments .. DOUBLE PRECISION a c .. c .. Local Scalars .. DOUBLE PRECISION bot,d,s1,s2,t,top,w c .. c .. Local Arrays .. DOUBLE PRECISION p(7),q(5),r(9) c .. c .. Data statements .. c ------------------- c ------------------- c ------------------- c ------------------- DATA p(1)/.577215664901533D+00/,p(2)/-.409078193005776D+00/, + p(3)/-.230975380857675D+00/,p(4)/.597275330452234D-01/, + p(5)/.766968181649490D-02/,p(6)/-.514889771323592D-02/, + p(7)/.589597428611429D-03/ DATA q(1)/.100000000000000D+01/,q(2)/.427569613095214D+00/, + q(3)/.158451672430138D+00/,q(4)/.261132021441447D-01/, + q(5)/.423244297896961D-02/ DATA r(1)/-.422784335098468D+00/,r(2)/-.771330383816272D+00/, + r(3)/-.244757765222226D+00/,r(4)/.118378989872749D+00/, + r(5)/.930357293360349D-03/,r(6)/-.118290993445146D-01/, + r(7)/.223047661158249D-02/,r(8)/.266505979058923D-03/, + r(9)/-.132674909766242D-03/ DATA s1/.273076135303957D+00/,s2/.559398236957378D-01/ c .. c .. Executable Statements .. c ------------------- t = a d = a - 0.5D0 IF (d.GT.0.0D0) t = d - 0.5D0 IF (t) 40,10,20 c 10 gam1 = 0.0D0 RETURN c 20 top = (((((p(7)*t+p(6))*t+p(5))*t+p(4))*t+p(3))*t+p(2))*t + p(1) bot = (((q(5)*t+q(4))*t+q(3))*t+q(2))*t + 1.0D0 w = top/bot IF (d.GT.0.0D0) GO TO 30 gam1 = a*w RETURN 30 gam1 = (t/a)* ((w-0.5D0)-0.5D0) RETURN c 40 top = (((((((r(9)*t+r(8))*t+r(7))*t+r(6))*t+r(5))*t+r(4))*t+r(3))* + t+r(2))*t + r(1) bot = (s2*t+s1)*t + 1.0D0 w = top/bot IF (d.GT.0.0D0) GO TO 50 gam1 = a* ((w+0.5D0)+0.5D0) RETURN 50 gam1 = t*w/a RETURN END c c DOUBLE PRECISION FUNCTION gamma(a) c----------------------------------------------------------------------- c c EVALUATION OF THE GAMMA FUNCTION FOR REAL ARGUMENTS c c ----------- c c GAMMA(A) IS ASSIGNED THE VALUE 0 WHEN THE GAMMA FUNCTION CANNOT c BE COMPUTED. c c----------------------------------------------------------------------- c WRITTEN BY ALFRED H. MORRIS, JR. c NAVAL SURFACE WEAPONS CENTER c DAHLGREN, VIRGINIA c----------------------------------------------------------------------- c .. Scalar Arguments .. DOUBLE PRECISION a c .. c .. Local Scalars .. DOUBLE PRECISION bot,d,g,lnx,pi,r1,r2,r3,r4,r5,s,t,top,w,x,z INTEGER i,j,m,n c .. c .. Local Arrays .. DOUBLE PRECISION p(7),q(7) c .. c .. External Functions .. DOUBLE PRECISION exparg,spmpar EXTERNAL exparg,spmpar c .. c .. Intrinsic Functions .. INTRINSIC abs,dble,dlog,exp,int,mod,sin c .. c .. Data statements .. c-------------------------- c D = 0.5*(LN(2*PI) - 1) c-------------------------- c-------------------------- c-------------------------- DATA pi/3.1415926535898D0/ DATA d/.41893853320467274178D0/ DATA p(1)/.539637273585445D-03/,p(2)/.261939260042690D-02/, + p(3)/.204493667594920D-01/,p(4)/.730981088720487D-01/, + p(5)/.279648642639792D+00/,p(6)/.553413866010467D+00/, + p(7)/1.0D0/ DATA q(1)/-.832979206704073D-03/,q(2)/.470059485860584D-02/, + q(3)/.225211131035340D-01/,q(4)/-.170458969313360D+00/, + q(5)/-.567902761974940D-01/,q(6)/.113062953091122D+01/, + q(7)/1.0D0/ DATA r1/.820756370353826D-03/,r2/-.595156336428591D-03/, + r3/.793650663183693D-03/,r4/-.277777777770481D-02/, + r5/.833333333333333D-01/ c .. c .. Executable Statements .. c-------------------------- gamma = 0.0D0 x = a IF (abs(a).GE.15.0D0) GO TO 110 c----------------------------------------------------------------------- c EVALUATION OF GAMMA(A) FOR ABS(A) .LT. 15 c----------------------------------------------------------------------- t = 1.0D0 m = int(a) - 1 c c LET T BE THE PRODUCT OF A-J WHEN A .GE. 2 c IF (m) 40,30,10 10 DO 20 j = 1,m x = x - 1.0D0 t = x*t 20 CONTINUE 30 x = x - 1.0D0 GO TO 80 c c LET T BE THE PRODUCT OF A+J WHEN A .LT. 1 c 40 t = a IF (a.GT.0.0D0) GO TO 70 m = -m - 1 IF (m.EQ.0) GO TO 60 DO 50 j = 1,m x = x + 1.0D0 t = x*t 50 CONTINUE 60 x = (x+0.5D0) + 0.5D0 t = x*t IF (t.EQ.0.0D0) RETURN c 70 CONTINUE c c THE FOLLOWING CODE CHECKS IF 1/T CAN OVERFLOW. THIS c CODE MAY BE OMITTED IF DESIRED. c IF (abs(t).GE.1.D-30) GO TO 80 IF (abs(t)*spmpar(3).LE.1.0001D0) RETURN gamma = 1.0D0/t RETURN c c COMPUTE GAMMA(1 + X) FOR 0 .LE. X .LT. 1 c 80 top = p(1) bot = q(1) DO 90 i = 2,7 top = p(i) + x*top bot = q(i) + x*bot 90 CONTINUE gamma = top/bot c c TERMINATION c IF (a.LT.1.0D0) GO TO 100 gamma = gamma*t RETURN 100 gamma = gamma/t RETURN c----------------------------------------------------------------------- c EVALUATION OF GAMMA(A) FOR ABS(A) .GE. 15 c----------------------------------------------------------------------- 110 IF (abs(a).GE.1.D3) RETURN IF (a.GT.0.0D0) GO TO 120 x = -a n = x t = x - n IF (t.GT.0.9D0) t = 1.0D0 - t s = sin(pi*t)/pi IF (mod(n,2).EQ.0) s = -s IF (s.EQ.0.0D0) RETURN c c COMPUTE THE MODIFIED ASYMPTOTIC SUM c 120 t = 1.0D0/ (x*x) g = ((((r1*t+r2)*t+r3)*t+r4)*t+r5)/x c c ONE MAY REPLACE THE NEXT STATEMENT WITH LNX = ALOG(X) c BUT LESS ACCURACY WILL NORMALLY BE OBTAINED. c lnx = dlog(x) c c FINAL ASSEMBLY c z = x g = (d+g) + (z-0.5D0)* (lnx-1.D0) w = g t = g - dble(w) IF (w.GT.0.99999D0*exparg(0)) RETURN gamma = exp(w)* (1.0D0+t) IF (a.LT.0.0D0) gamma = (1.0D0/ (gamma*s))/x RETURN END c c SUBROUTINE gratio(a,x,ans,qans,ind) c ---------------------------------------------------------------------- c EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS c P(A,X) AND Q(A,X) c c ---------- c c IT IS ASSUMED THAT A AND X ARE NONNEGATIVE, WHERE A AND X c ARE NOT BOTH 0. c c ANS AND QANS ARE VARIABLES. GRATIO ASSIGNS ANS THE VALUE c P(A,X) AND QANS THE VALUE Q(A,X). IND MAY BE ANY INTEGER. c IF IND = 0 THEN THE USER IS REQUESTING AS MUCH ACCURACY AS c POSSIBLE (UP TO 14 SIGNIFICANT DIGITS). OTHERWISE, IF c IND = 1 THEN ACCURACY IS REQUESTED TO WITHIN 1 UNIT OF THE c 6-TH SIGNIFICANT DIGIT, AND IF IND .NE. 0,1 THEN ACCURACY c IS REQUESTED TO WITHIN 1 UNIT OF THE 3RD SIGNIFICANT DIGIT. c c ERROR RETURN ... c ANS IS ASSIGNED THE VALUE 2 WHEN A OR X IS NEGATIVE, c WHEN A*X = 0, OR WHEN P(A,X) AND Q(A,X) ARE INDETERMINANT. c P(A,X) AND Q(A,X) ARE COMPUTATIONALLY INDETERMINANT WHEN c X IS EXCEEDINGLY CLOSE TO A AND A IS EXTREMELY LARGE. c ---------------------------------------------------------------------- c WRITTEN BY ALFRED H. MORRIS, JR. c NAVAL SURFACE WEAPONS CENTER c DAHLGREN, VIRGINIA c -------------------- c .. Scalar Arguments .. DOUBLE PRECISION a,ans,qans,x INTEGER ind c .. c .. Local Scalars .. DOUBLE PRECISION a2n,a2nm1,acc,alog10,am0,amn,an,an0,apn,b2n, + b2nm1,c,c0,c1,c2,c3,c4,c5,c6,cma,d10,d20,d30,d40, + d50,d60,d70,e,e0,g,h,j,l,r,rt2pin,rta,rtpi,rtx,s, + sum,t,t1,third,tol,twoa,u,w,x0,y,z INTEGER i,iop,m,max,n c .. c .. Local Arrays .. DOUBLE PRECISION acc0(3),big(3),d0(13),d1(12),d2(10),d3(8),d4(6), + d5(4),d6(2),e00(3),wk(20),x00(3) c .. c .. External Functions .. DOUBLE PRECISION erf,erfc1,gam1,gamma,rexp,rlog,spmpar EXTERNAL erf,erfc1,gam1,gamma,rexp,rlog,spmpar c .. c .. Intrinsic Functions .. INTRINSIC abs,dble,dlog,dmax1,exp,int,sqrt c .. c .. Data statements .. c -------------------- c -------------------- c ALOG10 = LN(10) c RT2PIN = 1/SQRT(2*PI) c RTPI = SQRT(PI) c -------------------- c -------------------- c -------------------- c -------------------- c -------------------- c -------------------- c -------------------- c -------------------- c -------------------- DATA acc0(1)/5.D-15/,acc0(2)/5.D-7/,acc0(3)/5.D-4/ DATA big(1)/20.0D0/,big(2)/14.0D0/,big(3)/10.0D0/ DATA e00(1)/.25D-3/,e00(2)/.25D-1/,e00(3)/.14D0/ DATA x00(1)/31.0D0/,x00(2)/17.0D0/,x00(3)/9.7D0/ DATA alog10/2.30258509299405D0/ DATA rt2pin/.398942280401433D0/ DATA rtpi/1.77245385090552D0/ DATA third/.333333333333333D0/ DATA d0(1)/.833333333333333D-01/,d0(2)/-.148148148148148D-01/, + d0(3)/.115740740740741D-02/,d0(4)/.352733686067019D-03/, + d0(5)/-.178755144032922D-03/,d0(6)/.391926317852244D-04/, + d0(7)/-.218544851067999D-05/,d0(8)/-.185406221071516D-05/, + d0(9)/.829671134095309D-06/,d0(10)/-.176659527368261D-06/, + d0(11)/.670785354340150D-08/,d0(12)/.102618097842403D-07/, + d0(13)/-.438203601845335D-08/ DATA d10/-.185185185185185D-02/,d1(1)/-.347222222222222D-02/, + d1(2)/.264550264550265D-02/,d1(3)/-.990226337448560D-03/, + d1(4)/.205761316872428D-03/,d1(5)/-.401877572016461D-06/, + d1(6)/-.180985503344900D-04/,d1(7)/.764916091608111D-05/, + d1(8)/-.161209008945634D-05/,d1(9)/.464712780280743D-08/, + d1(10)/.137863344691572D-06/,d1(11)/-.575254560351770D-07/, + d1(12)/.119516285997781D-07/ DATA d20/.413359788359788D-02/,d2(1)/-.268132716049383D-02/, + d2(2)/.771604938271605D-03/,d2(3)/.200938786008230D-05/, + d2(4)/-.107366532263652D-03/,d2(5)/.529234488291201D-04/, + d2(6)/-.127606351886187D-04/,d2(7)/.342357873409614D-07/, + d2(8)/.137219573090629D-05/,d2(9)/-.629899213838006D-06/, + d2(10)/.142806142060642D-06/ DATA d30/.649434156378601D-03/,d3(1)/.229472093621399D-03/, + d3(2)/-.469189494395256D-03/,d3(3)/.267720632062839D-03/, + d3(4)/-.756180167188398D-04/,d3(5)/-.239650511386730D-06/, + d3(6)/.110826541153473D-04/,d3(7)/-.567495282699160D-05/, + d3(8)/.142309007324359D-05/ DATA d40/-.861888290916712D-03/,d4(1)/.784039221720067D-03/, + d4(2)/-.299072480303190D-03/,d4(3)/-.146384525788434D-05/, + d4(4)/.664149821546512D-04/,d4(5)/-.396836504717943D-04/, + d4(6)/.113757269706784D-04/ DATA d50/-.336798553366358D-03/,d5(1)/-.697281375836586D-04/, + d5(2)/.277275324495939D-03/,d5(3)/-.199325705161888D-03/, + d5(4)/.679778047793721D-04/ DATA d60/.531307936463992D-03/,d6(1)/-.592166437353694D-03/, + d6(2)/.270878209671804D-03/ DATA d70/.344367606892378D-03/ c .. c .. Executable Statements .. c -------------------- c ****** E IS A MACHINE DEPENDENT CONSTANT. E IS THE SMALLEST c FLOATING POINT NUMBER FOR WHICH 1.0 + E .GT. 1.0 . c e = spmpar(1) c c -------------------- IF (a.LT.0.0D0 .OR. x.LT.0.0D0) GO TO 430 IF (a.EQ.0.0D0 .AND. x.EQ.0.0D0) GO TO 430 IF (a*x.EQ.0.0D0) GO TO 420 c iop = ind + 1 IF (iop.NE.1 .AND. iop.NE.2) iop = 3 acc = dmax1(acc0(iop),e) e0 = e00(iop) x0 = x00(iop) c c SELECT THE APPROPRIATE ALGORITHM c IF (a.GE.1.0D0) GO TO 10 IF (a.EQ.0.5D0) GO TO 390 IF (x.LT.1.1D0) GO TO 160 t1 = a*dlog(x) - x u = a*exp(t1) IF (u.EQ.0.0D0) GO TO 380 r = u* (1.0D0+gam1(a)) GO TO 250 c 10 IF (a.GE.big(iop)) GO TO 30 IF (a.GT.x .OR. x.GE.x0) GO TO 20 twoa = a + a m = int(twoa) IF (twoa.NE.dble(m)) GO TO 20 i = m/2 IF (a.EQ.dble(i)) GO TO 210 GO TO 220 20 t1 = a*dlog(x) - x r = exp(t1)/gamma(a) GO TO 40 c 30 l = x/a IF (l.EQ.0.0D0) GO TO 370 s = 0.5D0 + (0.5D0-l) z = rlog(l) IF (z.GE.700.0D0/a) GO TO 410 y = a*z rta = sqrt(a) IF (abs(s).LE.e0/rta) GO TO 330 IF (abs(s).LE.0.4D0) GO TO 270 c t = (1.0D0/a)**2 t1 = (((0.75D0*t-1.0D0)*t+3.5D0)*t-105.0D0)/ (a*1260.0D0) t1 = t1 - y r = rt2pin*rta*exp(t1) c 40 IF (r.EQ.0.0D0) GO TO 420 IF (x.LE.dmax1(a,alog10)) GO TO 50 IF (x.LT.x0) GO TO 250 GO TO 100 c c TAYLOR SERIES FOR P/R c 50 apn = a + 1.0D0 t = x/apn wk(1) = t DO 60 n = 2,20 apn = apn + 1.0D0 t = t* (x/apn) IF (t.LE.1.D-3) GO TO 70 wk(n) = t 60 CONTINUE n = 20 c 70 sum = t tol = 0.5D0*acc 80 apn = apn + 1.0D0 t = t* (x/apn) sum = sum + t IF (t.GT.tol) GO TO 80 c max = n - 1 DO 90 m = 1,max n = n - 1 sum = sum + wk(n) 90 CONTINUE ans = (r/a)* (1.0D0+sum) qans = 0.5D0 + (0.5D0-ans) RETURN c c ASYMPTOTIC EXPANSION c 100 amn = a - 1.0D0 t = amn/x wk(1) = t DO 110 n = 2,20 amn = amn - 1.0D0 t = t* (amn/x) IF (abs(t).LE.1.D-3) GO TO 120 wk(n) = t 110 CONTINUE n = 20 c 120 sum = t 130 IF (abs(t).LE.acc) GO TO 140 amn = amn - 1.0D0 t = t* (amn/x) sum = sum + t GO TO 130 c 140 max = n - 1 DO 150 m = 1,max n = n - 1 sum = sum + wk(n) 150 CONTINUE qans = (r/x)* (1.0D0+sum) ans = 0.5D0 + (0.5D0-qans) RETURN c c TAYLOR SERIES FOR P(A,X)/X**A c 160 an = 3.0D0 c = x sum = x/ (a+3.0D0) tol = 3.0D0*acc/ (a+1.0D0) 170 an = an + 1.0D0 c = -c* (x/an) t = c/ (a+an) sum = sum + t IF (abs(t).GT.tol) GO TO 170 j = a*x* ((sum/6.0D0-0.5D0/ (a+2.0D0))*x+1.0D0/ (a+1.0D0)) c z = a*dlog(x) h = gam1(a) g = 1.0D0 + h IF (x.LT.0.25D0) GO TO 180 IF (a.LT.x/2.59D0) GO TO 200 GO TO 190 180 IF (z.GT.-.13394D0) GO TO 200 c 190 w = exp(z) ans = w*g* (0.5D0+ (0.5D0-j)) qans = 0.5D0 + (0.5D0-ans) RETURN c 200 l = rexp(z) w = 0.5D0 + (0.5D0+l) qans = (w*j-l)*g - h IF (qans.LT.0.0D0) GO TO 380 ans = 0.5D0 + (0.5D0-qans) RETURN c c FINITE SUMS FOR Q WHEN A .GE. 1 c AND 2*A IS AN INTEGER c 210 sum = exp(-x) t = sum n = 1 c = 0.0D0 GO TO 230 c 220 rtx = sqrt(x) sum = erfc1(0,rtx) t = exp(-x)/ (rtpi*rtx) n = 0 c = -0.5D0 c 230 IF (n.EQ.i) GO TO 240 n = n + 1 c = c + 1.0D0 t = (x*t)/c sum = sum + t GO TO 230 240 qans = sum ans = 0.5D0 + (0.5D0-qans) RETURN c c CONTINUED FRACTION EXPANSION c 250 tol = dmax1(5.0D0*e,acc) a2nm1 = 1.0D0 a2n = 1.0D0 b2nm1 = x b2n = x + (1.0D0-a) c = 1.0D0 260 a2nm1 = x*a2n + c*a2nm1 b2nm1 = x*b2n + c*b2nm1 am0 = a2nm1/b2nm1 c = c + 1.0D0 cma = c - a a2n = a2nm1 + cma*a2n b2n = b2nm1 + cma*b2n an0 = a2n/b2n IF (abs(an0-am0).GE.tol*an0) GO TO 260 c qans = r*an0 ans = 0.5D0 + (0.5D0-qans) RETURN c c GENERAL TEMME EXPANSION c 270 IF (abs(s).LE.2.0D0*e .AND. a*e*e.GT.3.28D-3) GO TO 430 c = exp(-y) w = 0.5D0*erfc1(1,sqrt(y)) u = 1.0D0/a z = sqrt(z+z) IF (l.LT.1.0D0) z = -z IF (iop-2) 280,290,300 c 280 IF (abs(s).LE.1.D-3) GO TO 340 c0 = ((((((((((((d0(13)*z+d0(12))*z+d0(11))*z+d0(10))*z+d0(9))*z+ + d0(8))*z+d0(7))*z+d0(6))*z+d0(5))*z+d0(4))*z+d0(3))*z+d0(2))* + z+d0(1))*z - third c1 = (((((((((((d1(12)*z+d1(11))*z+d1(10))*z+d1(9))*z+d1(8))*z+ + d1(7))*z+d1(6))*z+d1(5))*z+d1(4))*z+d1(3))*z+d1(2))*z+d1(1))* + z + d10 c2 = (((((((((d2(10)*z+d2(9))*z+d2(8))*z+d2(7))*z+d2(6))*z+ + d2(5))*z+d2(4))*z+d2(3))*z+d2(2))*z+d2(1))*z + d20 c3 = (((((((d3(8)*z+d3(7))*z+d3(6))*z+d3(5))*z+d3(4))*z+d3(3))*z+ + d3(2))*z+d3(1))*z + d30 c4 = (((((d4(6)*z+d4(5))*z+d4(4))*z+d4(3))*z+d4(2))*z+d4(1))*z + + d40 c5 = (((d5(4)*z+d5(3))*z+d5(2))*z+d5(1))*z + d50 c6 = (d6(2)*z+d6(1))*z + d60 t = ((((((d70*u+c6)*u+c5)*u+c4)*u+c3)*u+c2)*u+c1)*u + c0 GO TO 310 c 290 c0 = (((((d0(6)*z+d0(5))*z+d0(4))*z+d0(3))*z+d0(2))*z+d0(1))*z - + third c1 = (((d1(4)*z+d1(3))*z+d1(2))*z+d1(1))*z + d10 c2 = d2(1)*z + d20 t = (c2*u+c1)*u + c0 GO TO 310 c 300 t = ((d0(3)*z+d0(2))*z+d0(1))*z - third c 310 IF (l.LT.1.0D0) GO TO 320 qans = c* (w+rt2pin*t/rta) ans = 0.5D0 + (0.5D0-qans) RETURN 320 ans = c* (w-rt2pin*t/rta) qans = 0.5D0 + (0.5D0-ans) RETURN c c TEMME EXPANSION FOR L = 1 c 330 IF (a*e*e.GT.3.28D-3) GO TO 430 c = 0.5D0 + (0.5D0-y) w = (0.5D0-sqrt(y)* (0.5D0+ (0.5D0-y/3.0D0))/rtpi)/c u = 1.0D0/a z = sqrt(z+z) IF (l.LT.1.0D0) z = -z IF (iop-2) 340,350,360 c 340 c0 = ((((((d0(7)*z+d0(6))*z+d0(5))*z+d0(4))*z+d0(3))*z+d0(2))*z+ + d0(1))*z - third c1 = (((((d1(6)*z+d1(5))*z+d1(4))*z+d1(3))*z+d1(2))*z+d1(1))*z + + d10 c2 = ((((d2(5)*z+d2(4))*z+d2(3))*z+d2(2))*z+d2(1))*z + d20 c3 = (((d3(4)*z+d3(3))*z+d3(2))*z+d3(1))*z + d30 c4 = (d4(2)*z+d4(1))*z + d40 c5 = (d5(2)*z+d5(1))*z + d50 c6 = d6(1)*z + d60 t = ((((((d70*u+c6)*u+c5)*u+c4)*u+c3)*u+c2)*u+c1)*u + c0 GO TO 310 c 350 c0 = (d0(2)*z+d0(1))*z - third c1 = d1(1)*z + d10 t = (d20*u+c1)*u + c0 GO TO 310 c 360 t = d0(1)*z - third GO TO 310 c c SPECIAL CASES c 370 ans = 0.0D0 qans = 1.0D0 RETURN c 380 ans = 1.0D0 qans = 0.0D0 RETURN c 390 IF (x.GE.0.25D0) GO TO 400 ans = erf(sqrt(x)) qans = 0.5D0 + (0.5D0-ans) RETURN 400 qans = erfc1(0,sqrt(x)) ans = 0.5D0 + (0.5D0-qans) RETURN c 410 IF (abs(s).LE.2.0D0*e) GO TO 430 420 IF (x.LE.a) GO TO 370 GO TO 380 c c ERROR RETURN c 430 ans = 2.0D0 RETURN END c c INTEGER FUNCTION ipmpar(i) c----------------------------------------------------------------------- c c IPMPAR PROVIDES THE INTEGER MACHINE CONSTANTS FOR THE COMPUTER c THAT IS USED. IT IS ASSUMED THAT THE ARGUMENT I IS AN INTEGER c HAVING ONE OF THE VALUES 1-10. IPMPAR(I) HAS THE VALUE ... c c INTEGERS. c c ASSUME INTEGERS ARE REPRESENTED IN THE N-DIGIT, BASE-A FORM c c SIGN ( X(N-1)*A**(N-1) + ... + X(1)*A + X(0) ) c c WHERE 0 .LE. X(I) .LT. A FOR I=0,...,N-1. c c IPMPAR(1) = A, THE BASE. c c IPMPAR(2) = N, THE NUMBER OF BASE-A DIGITS. c c IPMPAR(3) = A**N - 1, THE LARGEST MAGNITUDE. c c FLOATING-POINT NUMBERS. c c IT IS ASSUMED THAT THE SINGLE AND DOUBLE PRECISION FLOATING c POINT ARITHMETICS HAVE THE SAME BASE, SAY B, AND THAT THE c NONZERO NUMBERS ARE REPRESENTED IN THE FORM c c SIGN (B**E) * (X(1)/B + ... + X(M)/B**M) c c WHERE X(I) = 0,1,...,B-1 FOR I=1,...,M, c X(1) .GE. 1, AND EMIN .LE. E .LE. EMAX. c c IPMPAR(4) = B, THE BASE. c c SINGLE-PRECISION c c IPMPAR(5) = M, THE NUMBER OF BASE-B DIGITS. c c IPMPAR(6) = EMIN, THE SMALLEST EXPONENT E. c c IPMPAR(7) = EMAX, THE LARGEST EXPONENT E. c c DOUBLE-PRECISION c c IPMPAR(8) = M, THE NUMBER OF BASE-B DIGITS. c c IPMPAR(9) = EMIN, THE SMALLEST EXPONENT E. c c IPMPAR(10) = EMAX, THE LARGEST EXPONENT E. c c----------------------------------------------------------------------- c c TO DEFINE THIS FUNCTION FOR THE COMPUTER BEING USED, ACTIVATE c THE DATA STATMENTS FOR THE COMPUTER BY REMOVING THE C FROM c COLUMN 1. (ALL THE OTHER DATA STATEMENTS SHOULD HAVE C IN c COLUMN 1.) c c----------------------------------------------------------------------- c c IPMPAR IS AN ADAPTATION OF THE FUNCTION I1MACH, WRITTEN BY c P.A. FOX, A.D. HALL, AND N.L. SCHRYER (BELL LABORATORIES). c IPMPAR WAS FORMED BY A.H. MORRIS (NSWC). THE CONSTANTS ARE c FROM BELL LABORATORIES, NSWC, AND OTHER SOURCES. c c----------------------------------------------------------------------- c .. Scalar Arguments .. INTEGER i c .. c .. Local Arrays .. INTEGER imach(10) c .. c .. Data statements .. c c MACHINE CONSTANTS FOR AMDAHL MACHINES. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 31 / c DATA IMACH( 3) / 2147483647 / c DATA IMACH( 4) / 16 / c DATA IMACH( 5) / 6 / c DATA IMACH( 6) / -64 / c DATA IMACH( 7) / 63 / c DATA IMACH( 8) / 14 / c DATA IMACH( 9) / -64 / c DATA IMACH(10) / 63 / c c MACHINE CONSTANTS FOR THE AT&T 3B SERIES, AT&T c PC 7300, AND AT&T 6300. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 31 / c DATA IMACH( 3) / 2147483647 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 24 / c DATA IMACH( 6) / -125 / c DATA IMACH( 7) / 128 / c DATA IMACH( 8) / 53 / c DATA IMACH( 9) / -1021 / c DATA IMACH(10) / 1024 / c c MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 33 / c DATA IMACH( 3) / 8589934591 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 24 / c DATA IMACH( 6) / -256 / c DATA IMACH( 7) / 255 / c DATA IMACH( 8) / 60 / c DATA IMACH( 9) / -256 / c DATA IMACH(10) / 255 / c c MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 39 / c DATA IMACH( 3) / 549755813887 / c DATA IMACH( 4) / 8 / c DATA IMACH( 5) / 13 / c DATA IMACH( 6) / -50 / c DATA IMACH( 7) / 76 / c DATA IMACH( 8) / 26 / c DATA IMACH( 9) / -50 / c DATA IMACH(10) / 76 / c c MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 39 / c DATA IMACH( 3) / 549755813887 / c DATA IMACH( 4) / 8 / c DATA IMACH( 5) / 13 / c DATA IMACH( 6) / -50 / c DATA IMACH( 7) / 76 / c DATA IMACH( 8) / 26 / c DATA IMACH( 9) / -32754 / c DATA IMACH(10) / 32780 / c c MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES c 60 BIT ARITHMETIC, AND THE CDC CYBER 995 64 BIT c ARITHMETIC (NOS OPERATING SYSTEM). c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 48 / c DATA IMACH( 3) / 281474976710655 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 48 / c DATA IMACH( 6) / -974 / c DATA IMACH( 7) / 1070 / c DATA IMACH( 8) / 95 / c DATA IMACH( 9) / -926 / c DATA IMACH(10) / 1070 / c c MACHINE CONSTANTS FOR THE CDC CYBER 995 64 BIT c ARITHMETIC (NOS/VE OPERATING SYSTEM). c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 63 / c DATA IMACH( 3) / 9223372036854775807 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 48 / c DATA IMACH( 6) / -4096 / c DATA IMACH( 7) / 4095 / c DATA IMACH( 8) / 96 / c DATA IMACH( 9) / -4096 / c DATA IMACH(10) / 4095 / c c MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 63 / c DATA IMACH( 3) / 9223372036854775807 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 47 / c DATA IMACH( 6) / -8189 / c DATA IMACH( 7) / 8190 / c DATA IMACH( 8) / 94 / c DATA IMACH( 9) / -8099 / c DATA IMACH(10) / 8190 / c c MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 15 / c DATA IMACH( 3) / 32767 / c DATA IMACH( 4) / 16 / c DATA IMACH( 5) / 6 / c DATA IMACH( 6) / -64 / c DATA IMACH( 7) / 63 / c DATA IMACH( 8) / 14 / c DATA IMACH( 9) / -64 / c DATA IMACH(10) / 63 / c c MACHINE CONSTANTS FOR THE HARRIS 220. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 23 / c DATA IMACH( 3) / 8388607 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 23 / c DATA IMACH( 6) / -127 / c DATA IMACH( 7) / 127 / c DATA IMACH( 8) / 38 / c DATA IMACH( 9) / -127 / c DATA IMACH(10) / 127 / c c MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 c AND DPS 8/70 SERIES. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 35 / c DATA IMACH( 3) / 34359738367 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 27 / c DATA IMACH( 6) / -127 / c DATA IMACH( 7) / 127 / c DATA IMACH( 8) / 63 / c DATA IMACH( 9) / -127 / c DATA IMACH(10) / 127 / c c MACHINE CONSTANTS FOR THE HP 2100 c 3 WORD DOUBLE PRECISION OPTION WITH FTN4 c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 15 / c DATA IMACH( 3) / 32767 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 23 / c DATA IMACH( 6) / -128 / c DATA IMACH( 7) / 127 / c DATA IMACH( 8) / 39 / c DATA IMACH( 9) / -128 / c DATA IMACH(10) / 127 / c c MACHINE CONSTANTS FOR THE HP 2100 c 4 WORD DOUBLE PRECISION OPTION WITH FTN4 c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 15 / c DATA IMACH( 3) / 32767 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 23 / c DATA IMACH( 6) / -128 / c DATA IMACH( 7) / 127 / c DATA IMACH( 8) / 55 / c DATA IMACH( 9) / -128 / c DATA IMACH(10) / 127 / c c MACHINE CONSTANTS FOR THE HP 9000. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 31 / c DATA IMACH( 3) / 2147483647 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 24 / c DATA IMACH( 6) / -126 / c DATA IMACH( 7) / 128 / c DATA IMACH( 8) / 53 / c DATA IMACH( 9) / -1021 / c DATA IMACH(10) / 1024 / c c MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, c THE ICL 2900, THE ITEL AS/6, THE XEROX SIGMA c 5/7/9 AND THE SEL SYSTEMS 85/86. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 31 / c DATA IMACH( 3) / 2147483647 / c DATA IMACH( 4) / 16 / c DATA IMACH( 5) / 6 / c DATA IMACH( 6) / -64 / c DATA IMACH( 7) / 63 / c DATA IMACH( 8) / 14 / c DATA IMACH( 9) / -64 / c DATA IMACH(10) / 63 / c c MACHINE CONSTANTS FOR THE IBM PC. c c DATA imach(1)/2/ c DATA imach(2)/31/ c DATA imach(3)/2147483647/ c DATA imach(4)/2/ c DATA imach(5)/24/ c DATA imach(6)/-125/ c DATA imach(7)/128/ c DATA imach(8)/53/ c DATA imach(9)/-1021/ c DATA imach(10)/1024/ c c MACHINE CONSTANTS FOR THE MACINTOSH II - ABSOFT c MACFORTRAN II. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 31 / c DATA IMACH( 3) / 2147483647 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 24 / c DATA IMACH( 6) / -125 / c DATA IMACH( 7) / 128 / c DATA IMACH( 8) / 53 / c DATA IMACH( 9) / -1021 / c DATA IMACH(10) / 1024 / c c MACHINE CONSTANTS FOR THE MICROVAX - S FORTRAN. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 31 / c DATA IMACH( 3) / 2147483647 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 24 / c DATA IMACH( 6) / -127 / c DATA IMACH( 7) / 127 / c DATA IMACH( 8) / 56 / c DATA IMACH( 9) / -127 / c DATA IMACH(10) / 127 / c c MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 35 / c DATA IMACH( 3) / 34359738367 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 27 / c DATA IMACH( 6) / -128 / c DATA IMACH( 7) / 127 / c DATA IMACH( 8) / 54 / c DATA IMACH( 9) / -101 / c DATA IMACH(10) / 127 / c c MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 35 / c DATA IMACH( 3) / 34359738367 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 27 / c DATA IMACH( 6) / -128 / c DATA IMACH( 7) / 127 / c DATA IMACH( 8) / 62 / c DATA IMACH( 9) / -128 / c DATA IMACH(10) / 127 / c c MACHINE CONSTANTS FOR THE PDP-11 FORTRAN SUPPORTING c 32-BIT INTEGER ARITHMETIC. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 31 / c DATA IMACH( 3) / 2147483647 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 24 / c DATA IMACH( 6) / -127 / c DATA IMACH( 7) / 127 / c DATA IMACH( 8) / 56 / c DATA IMACH( 9) / -127 / c DATA IMACH(10) / 127 / c c MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 31 / c DATA IMACH( 3) / 2147483647 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 24 / c DATA IMACH( 6) / -125 / c DATA IMACH( 7) / 128 / c DATA IMACH( 8) / 53 / c DATA IMACH( 9) / -1021 / c DATA IMACH(10) / 1024 / c c MACHINE CONSTANTS FOR THE SILICON GRAPHICS IRIS-4D c SERIES (MIPS R3000 PROCESSOR). c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 31 / c DATA IMACH( 3) / 2147483647 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 24 / c DATA IMACH( 6) / -125 / c DATA IMACH( 7) / 128 / c DATA IMACH( 8) / 53 / c DATA IMACH( 9) / -1021 / c DATA IMACH(10) / 1024 / c c MACHINE CONSTANTS FOR 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 DATA IMACH( 1) / 2 / DATA IMACH( 2) / 31 / DATA IMACH( 3) / 2147483647 / DATA IMACH( 4) / 2 / DATA IMACH( 5) / 24 / DATA IMACH( 6) / -125 / DATA IMACH( 7) / 128 / DATA IMACH( 8) / 53 / DATA IMACH( 9) / -1021 / DATA IMACH(10) / 1024 / c c MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 35 / c DATA IMACH( 3) / 34359738367 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 27 / c DATA IMACH( 6) / -128 / c DATA IMACH( 7) / 127 / c DATA IMACH( 8) / 60 / c DATA IMACH( 9) /-1024 / c DATA IMACH(10) / 1023 / c c MACHINE CONSTANTS FOR THE VAX 11/780. c c DATA IMACH( 1) / 2 / c DATA IMACH( 2) / 31 / c DATA IMACH( 3) / 2147483647 / c DATA IMACH( 4) / 2 / c DATA IMACH( 5) / 24 / c DATA IMACH( 6) / -127 / c DATA IMACH( 7) / 127 / c DATA IMACH( 8) / 56 / c DATA IMACH( 9) / -127 / c DATA IMACH(10) / 127 / c ipmpar = imach(i) RETURN END c c DOUBLE PRECISION FUNCTION rexp(x) c----------------------------------------------------------------------- c EVALUATION OF THE FUNCTION EXP(X) - 1 c----------------------------------------------------------------------- c .. Scalar Arguments .. DOUBLE PRECISION x c .. c .. Local Scalars .. DOUBLE PRECISION p1,p2,q1,q2,q3,q4,w c .. c .. Intrinsic Functions .. INTRINSIC abs,exp c .. c .. Data statements .. DATA p1/.914041914819518D-09/,p2/.238082361044469D-01/, + q1/-.499999999085958D+00/,q2/.107141568980644D+00/, + q3/-.119041179760821D-01/,q4/.595130811860248D-03/ c .. c .. Executable Statements .. c----------------------- IF (abs(x).GT.0.15D0) GO TO 10 rexp = x* (((p2*x+p1)*x+1.0D0)/ ((((q4*x+q3)*x+q2)*x+q1)*x+1.0D0)) RETURN c 10 w = exp(x) IF (x.GT.0.0D0) GO TO 20 rexp = (w-0.5D0) - 0.5D0 RETURN 20 rexp = w* (0.5D0+ (0.5D0-1.0D0/w)) RETURN END c c DOUBLE PRECISION FUNCTION rlog(x) c ------------------- c COMPUTATION OF X - 1 - LN(X) c ------------------- c .. Scalar Arguments .. DOUBLE PRECISION x c .. c .. Local Scalars .. DOUBLE PRECISION a,b,p0,p1,p2,q1,q2,r,t,u,w,w1 c .. c .. Intrinsic Functions .. INTRINSIC dble,dlog c .. c .. Data statements .. c ------------------- DATA a/.566749439387324D-01/ DATA b/.456512608815524D-01/ DATA p0/.333333333333333D+00/,p1/-.224696413112536D+00/, + p2/.620886815375787D-02/ DATA q1/-.127408923933623D+01/,q2/.354508718369557D+00/ c .. c .. Executable Statements .. c ------------------- IF (x.LT.0.61D0 .OR. x.GT.1.57D0) GO TO 40 IF (x.LT.0.82D0) GO TO 10 IF (x.GT.1.18D0) GO TO 20 c c ARGUMENT REDUCTION c u = (x-0.5D0) - 0.5D0 w1 = 0.0D0 GO TO 30 c 10 u = dble(x) - 0.7D0 u = u/0.7D0 w1 = a - u*0.3D0 GO TO 30 c 20 u = 0.75D0*dble(x) - 1.D0 w1 = b + u/3.0D0 c c SERIES EXPANSION c 30 r = u/ (u+2.0D0) t = r*r w = ((p2*t+p1)*t+p0)/ ((q2*t+q1)*t+1.0D0) rlog = 2.0D0*t* (1.0D0/ (1.0D0-r)-r*w) + w1 RETURN c c 40 r = (x-0.5D0) - 0.5D0 rlog = r - dlog(x) RETURN END c c DOUBLE PRECISION FUNCTION spmpar(i) c----------------------------------------------------------------------- c c SPMPAR PROVIDES THE SINGLE PRECISION MACHINE CONSTANTS FOR c THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT c I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE c SINGLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND c ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN c c SPMPAR(1) = B**(1 - M), THE MACHINE PRECISION, c c SPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE, c c SPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE. c c----------------------------------------------------------------------- c WRITTEN BY c ALFRED H. MORRIS, JR. c NAVAL SURFACE WARFARE CENTER c DAHLGREN VIRGINIA c----------------------------------------------------------------------- c----------------------------------------------------------------------- c MODIFIED BY BARRY W. BROWN TO RETURN DOUBLE PRECISION MACHINE c CONSTANTS FOR THE COMPUTER BEING USED. THIS MODIFICATION WAS c MADE AS PART OF CONVERTING BRATIO TO DOUBLE PRECISION c----------------------------------------------------------------------- c .. Scalar Arguments .. INTEGER i c .. c .. Local Scalars .. DOUBLE PRECISION b,binv,bm1,one,w,z INTEGER emax,emin,ibeta,m c .. c .. External Functions .. INTEGER ipmpar EXTERNAL ipmpar c .. c .. Intrinsic Functions .. INTRINSIC dble c .. c .. Executable Statements .. c IF (i.GT.1) GO TO 10 b = ipmpar(4) m = ipmpar(8) spmpar = b** (1-m) RETURN c 10 IF (i.GT.2) GO TO 20 b = ipmpar(4) emin = ipmpar(9) one = dble(1) binv = one/b w = b** (emin+2) spmpar = ((w*binv)*binv)*binv RETURN c 20 ibeta = ipmpar(4) m = ipmpar(8) emax = ipmpar(10) c b = ibeta bm1 = ibeta - 1 one = dble(1) z = b** (m-1) w = ((z-one)*b+bm1)/ (b*z) c z = b** (emax-2) spmpar = ((w*z)*b)*b RETURN END