DOUBLE PRECISION FUNCTION psi(xx) C--------------------------------------------------------------------- C C EVALUATION OF THE DIGAMMA FUNCTION C C ----------- C C PSI(XX) IS ASSIGNED THE VALUE 0 WHEN THE DIGAMMA FUNCTION CANNOT C BE COMPUTED. C C THE MAIN COMPUTATION INVOLVES EVALUATION OF RATIONAL CHEBYSHEV C APPROXIMATIONS PUBLISHED IN MATH. COMP. 27, 123-127(1973) BY C CODY, STRECOK AND THACHER. C C--------------------------------------------------------------------- C PSI WAS WRITTEN AT ARGONNE NATIONAL LABORATORY FOR THE FUNPACK C PACKAGE OF SPECIAL FUNCTION SUBROUTINES. PSI WAS MODIFIED BY C A.H. MORRIS (NSWC). C--------------------------------------------------------------------- C .. Scalar Arguments .. DOUBLE PRECISION xx C .. C .. Local Scalars .. DOUBLE PRECISION aug,den,dx0,piov4,sgn,upper,w,x,xmax1,xmx0, + xsmall,z INTEGER i,m,n,nq C .. C .. Local Arrays .. DOUBLE PRECISION p1(7),p2(4),q1(6),q2(4) C .. C .. External Functions .. DOUBLE PRECISION spmpar INTEGER ipmpar EXTERNAL spmpar,ipmpar C .. C .. Intrinsic Functions .. INTRINSIC abs,cos,dble,dlog,dmin1,int,sin C .. C .. Data statements .. C--------------------------------------------------------------------- C C PIOV4 = PI/4 C DX0 = ZERO OF PSI TO EXTENDED PRECISION C C--------------------------------------------------------------------- C--------------------------------------------------------------------- C C COEFFICIENTS FOR RATIONAL APPROXIMATION OF C PSI(X) / (X - X0), 0.5 .LE. X .LE. 3.0 C C--------------------------------------------------------------------- C--------------------------------------------------------------------- C C COEFFICIENTS FOR RATIONAL APPROXIMATION OF C PSI(X) - LN(X) + 1 / (2*X), X .GT. 3.0 C C--------------------------------------------------------------------- DATA piov4/.785398163397448D0/ DATA dx0/1.461632144968362341262659542325721325D0/ DATA p1(1)/.895385022981970D-02/,p1(2)/.477762828042627D+01/, + p1(3)/.142441585084029D+03/,p1(4)/.118645200713425D+04/, + p1(5)/.363351846806499D+04/,p1(6)/.413810161269013D+04/, + p1(7)/.130560269827897D+04/ DATA q1(1)/.448452573429826D+02/,q1(2)/.520752771467162D+03/, + q1(3)/.221000799247830D+04/,q1(4)/.364127349079381D+04/, + q1(5)/.190831076596300D+04/,q1(6)/.691091682714533D-05/ DATA p2(1)/-.212940445131011D+01/,p2(2)/-.701677227766759D+01/, + p2(3)/-.448616543918019D+01/,p2(4)/-.648157123766197D+00/ DATA q2(1)/.322703493791143D+02/,q2(2)/.892920700481861D+02/, + q2(3)/.546117738103215D+02/,q2(4)/.777788548522962D+01/ C .. C .. Executable Statements .. C--------------------------------------------------------------------- C C MACHINE DEPENDENT CONSTANTS ... C C XMAX1 = THE SMALLEST POSITIVE FLOATING POINT CONSTANT C WITH ENTIRELY INTEGER REPRESENTATION. ALSO USED C AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE C ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH C PSI MAY BE REPRESENTED AS ALOG(X). C C XSMALL = ABSOLUTE ARGUMENT BELOW WHICH PI*COTAN(PI*X) C MAY BE REPRESENTED BY 1/X. C C--------------------------------------------------------------------- xmax1 = ipmpar(3) xmax1 = dmin1(xmax1,1.0D0/spmpar(1)) xsmall = 1.D-9 C--------------------------------------------------------------------- x = xx aug = 0.0D0 IF (x.GE.0.5D0) GO TO 50 C--------------------------------------------------------------------- C X .LT. 0.5, USE REFLECTION FORMULA C PSI(1-X) = PSI(X) + PI * COTAN(PI*X) C--------------------------------------------------------------------- IF (abs(x).GT.xsmall) GO TO 10 IF (x.EQ.0.0D0) GO TO 100 C--------------------------------------------------------------------- C 0 .LT. ABS(X) .LE. XSMALL. USE 1/X AS A SUBSTITUTE C FOR PI*COTAN(PI*X) C--------------------------------------------------------------------- aug = -1.0D0/x GO TO 40 C--------------------------------------------------------------------- C REDUCTION OF ARGUMENT FOR COTAN C--------------------------------------------------------------------- 10 w = -x sgn = piov4 IF (w.GT.0.0D0) GO TO 20 w = -w sgn = -sgn C--------------------------------------------------------------------- C MAKE AN ERROR EXIT IF X .LE. -XMAX1 C--------------------------------------------------------------------- 20 IF (w.GE.xmax1) GO TO 100 nq = int(w) w = w - dble(nq) nq = int(w*4.0D0) w = 4.0D0* (w-dble(nq)*.25D0) C--------------------------------------------------------------------- C W IS NOW RELATED TO THE FRACTIONAL PART OF 4.0 * X. C ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST C QUADRANT AND DETERMINE SIGN C--------------------------------------------------------------------- n = nq/2 IF ((n+n).NE.nq) w = 1.0D0 - w z = piov4*w m = n/2 IF ((m+m).NE.n) sgn = -sgn C--------------------------------------------------------------------- C DETERMINE FINAL VALUE FOR -PI*COTAN(PI*X) C--------------------------------------------------------------------- n = (nq+1)/2 m = n/2 m = m + m IF (m.NE.n) GO TO 30 C--------------------------------------------------------------------- C CHECK FOR SINGULARITY C--------------------------------------------------------------------- IF (z.EQ.0.0D0) GO TO 100 C--------------------------------------------------------------------- C USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND C SIN/COS AS A SUBSTITUTE FOR TAN C--------------------------------------------------------------------- aug = sgn* ((cos(z)/sin(z))*4.0D0) GO TO 40 30 aug = sgn* ((sin(z)/cos(z))*4.0D0) 40 x = 1.0D0 - x 50 IF (x.GT.3.0D0) GO TO 70 C--------------------------------------------------------------------- C 0.5 .LE. X .LE. 3.0 C--------------------------------------------------------------------- den = x upper = p1(1)*x C DO 60 i = 1,5 den = (den+q1(i))*x upper = (upper+p1(i+1))*x 60 CONTINUE C den = (upper+p1(7))/ (den+q1(6)) xmx0 = dble(x) - dx0 psi = den*xmx0 + aug RETURN C--------------------------------------------------------------------- C IF X .GE. XMAX1, PSI = LN(X) C--------------------------------------------------------------------- 70 IF (x.GE.xmax1) GO TO 90 C--------------------------------------------------------------------- C 3.0 .LT. X .LT. XMAX1 C--------------------------------------------------------------------- w = 1.0D0/ (x*x) den = w upper = p2(1)*w C DO 80 i = 1,3 den = (den+q2(i))*w upper = (upper+p2(i+1))*w 80 CONTINUE C aug = upper/ (den+q2(4)) - 0.5D0/x + aug 90 psi = aug + dlog(x) RETURN C--------------------------------------------------------------------- C ERROR RETURN C--------------------------------------------------------------------- 100 psi = 0.0D0 RETURN END