      REAL FUNCTION PSI(XX)
      REAL AUG,DEN,P1(7),P2(4),PIOV4,Q1(6),Q2(4),SGN,XSMALL,UPPER,
     1     W,X,XMX0,XMAX1,XX,Z
      INTEGER I,M,N,NQ
      DOUBLE PRECISION DX0
C
C     THIS SUBROUTINE COMPUTES SINGLE PRECISION VALUES OF THE
C     LOGARITHMIC DERIVATIVE OF THE GAMMA FUNCTION,
C
C      PSI(X) = (D GAMMA(X) / DX ) / GAMMA(X) = D( LN GAMMA(X) ) / DX,
C
C     IT IS ASSUMED THAT X.NE.0, X IS NOT A NEGATIVE INTEGER, AND
C     X.GT.-XMAX1. (XMAX1 IS A MACHINE DEPENDENT CONSTANT.)
C
C     THE FORTRAN CALLING STATEMENT FOR THIS SUBROUTINE IS
C
C                  Y = PSI(X)
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     OTHER SUBPROGRAMS REQUIRED
C
C            ABS, ALOG, COS, DBLE, FLOAT, INT, SIN, SPMPAR
C
C       ----------
C
C     ERROR RETURN
C
C     IF X = 0, X IS A NEGATIVE INTEGER, OR X.LE.-XMAX1 THEN PSI(X)
C     IS ASSIGNED THE VALUE 0.
C
C        ----------
C
C     PSI WAS WRITTEN AT ARGONNE NATIONAL LABORATORY. THE FORMAT OF
C     THE FUNCTION WAS MODIFIED BY A.H. MORRIS (NAVAL SURFACE WEAPONS
C     CENTER, DAHLGREN, VIRGINIA).
C
C---------------------------------------------------------------------
C
C     PIOV4 = PI/4
C     DX0 = ZERO OF PSI TO EXTENDED PRECISION
C
C---------------------------------------------------------------------
      DATA PIOV4/.785398163397448E0/
      DATA DX0/1.461632144968362341262659542325721325D0/
C---------------------------------------------------------------------
C
C     COEFFICIENTS FOR RATIONAL APPROXIMATION OF
C     PSI(X) / (X - X0),  0.5 .LE. X .LE. 3.0
C
C---------------------------------------------------------------------
      DATA P1(1)/.895385022981970E-02/,  P1(2)/.477762828042627E+01/,
     1     P1(3)/.142441585084029E+03/,  P1(4)/.118645200713425E+04/,
     2     P1(5)/.363351846806499E+04/,  P1(6)/.413810161269013E+04/,
     3     P1(7)/.130560269827897E+04/
      DATA Q1(1)/.448452573429826E+02/,  Q1(2)/.520752771467162E+03/,
     1     Q1(3)/.221000799247830E+04/,  Q1(4)/.364127349079381E+04/,
     2     Q1(5)/.190831076596300E+04/,  Q1(6)/.691091682714533E-05/
C---------------------------------------------------------------------
C
C     COEFFICIENTS FOR RATIONAL APPROXIMATION OF
C     PSI(X) - LN(X) + 1 / (2*X),  X .GT. 3.0
C
C---------------------------------------------------------------------
      DATA P2(1)/-.212940445131011E+01/, P2(2)/-.701677227766759E+01/,
     1     P2(3)/-.448616543918019E+01/, P2(4)/-.648157123766197E+00/
      DATA Q2(1)/ .322703493791143E+02/, Q2(2)/ .892920700481861E+02/,
     1     Q2(3)/ .546117738103215E+02/, Q2(4)/ .777788548522962E+01/
C---------------------------------------------------------------------
C
C     MACHINE DEPENDENT CONSTANTS
C
C        XMAX1 = THE SMALLEST POSITIVE FLOATING POINT CONSTANT
C                WITH ENTIRELY INTEGER REPRESENTATION.  ALSO USED
C                AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE
C                ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH
C                PSI MAY BE REPRESENTED AS ALOG(X)
C        XSMALL = 2**(-26), ABSOLUTE ARGUMENT BELOW WHICH
C                 PI*COTAN(PI*X) MAY BE REPRESENTED BY  1/X.
C
C---------------------------------------------------------------------
      XMAX1 = 1.0/SPMPAR(1)
      XSMALL = 1.490116E-08
C
C---------------------------------------------------------------------
C
      X = XX
      AUG = 0.0E0
      IF (X .GE. 0.5E0) GO TO 200
C---------------------------------------------------------------------
C     X .LT. 0.5,  USE REFLECTION FORMULA
C     PSI(1-X) = PSI(X) + PI * COTAN(PI*X)
C---------------------------------------------------------------------
      IF (ABS(X) .GT. XSMALL) GO TO 100
      IF (X .EQ. 0.0E0) GO TO 400
C---------------------------------------------------------------------
C     0 .LT. ABS(X) .LE. XSMALL.  USE 1/X AS A SUBSTITUTE
C     FOR  PI*COTAN(PI*X)
C---------------------------------------------------------------------
      AUG = -1.0E0 / X
      GO TO 150
C---------------------------------------------------------------------
C     REDUCTION OF ARGUMENT FOR COTAN
C---------------------------------------------------------------------
  100 W = - X
      SGN = PIOV4
      IF (W .GT. 0.0E0) GO TO 120
      W = - W
      SGN = -SGN
C---------------------------------------------------------------------
C     MAKE AN ERROR EXIT IF X .LE. -XMAX1
C---------------------------------------------------------------------
  120 IF (W .GE. XMAX1) GO TO 400
      NQ = INT(W)
      W = W - FLOAT(NQ)
      NQ = INT(W*4.0E0)
      W = 4.0E0 * (W - FLOAT(NQ) * .25E0)
C---------------------------------------------------------------------
C     W IS NOW RELATED TO THE FRACTIONAL PART OF  4.0 * X.
C     ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST
C     QUADRANT AND DETERMINE SIGN
C---------------------------------------------------------------------
      N = NQ / 2
      IF ((N+N) .NE. NQ) W = 1.0E0 - W
      Z = PIOV4 * W
      M = N / 2
      IF ((M+M) .NE. N) SGN = - SGN
C---------------------------------------------------------------------
C     DETERMINE FINAL VALUE FOR  -PI*COTAN(PI*X)
C---------------------------------------------------------------------
      N = (NQ + 1) / 2
      M = N / 2
      M = M + M
      IF (M .NE. N) GO TO 140
C---------------------------------------------------------------------
C     CHECK FOR SINGULARITY
C---------------------------------------------------------------------
      IF (Z .EQ. 0.0E0) GO TO 400
C---------------------------------------------------------------------
C     USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND
C     SIN/COS AS A SUBSTITUTE FOR TAN
C---------------------------------------------------------------------
      AUG = SGN * ((COS(Z) / SIN(Z)) * 4.0E0)
      GO TO 150
  140 AUG = SGN * ((SIN(Z) / COS(Z)) * 4.0E0)
  150 X = 1.0E0 - X
  200 IF (X .GT. 3.0E0) GO TO 300
C---------------------------------------------------------------------
C     0.5 .LE. X .LE. 3.0
C---------------------------------------------------------------------
      DEN = X
      UPPER = P1(1) * X
C
      DO 210 I = 1, 5
         DEN = (DEN + Q1(I)) * X
         UPPER = (UPPER + P1(I+1)) * X
  210 CONTINUE
C
      DEN = (UPPER + P1(7)) / (DEN + Q1(6))
      XMX0 = DBLE(X) - DX0
      PSI = DEN * XMX0 + AUG
      RETURN
C---------------------------------------------------------------------
C     IF X .GE. XMAX1, PSI = LN(X)
C---------------------------------------------------------------------
  300 IF (X .GE. XMAX1) GO TO 350
C---------------------------------------------------------------------
C     3.0 .LT. X .LT. XMAX1
C---------------------------------------------------------------------
      W = 1.0E0 / (X * X)
      DEN = W
      UPPER = P2(1) * W
C
      DO 310 I = 1, 3
         DEN = (DEN + Q2(I)) * W
         UPPER = (UPPER + P2(I+1)) * W
  310 CONTINUE
C
      AUG = UPPER / (DEN + Q2(4)) - 0.5E0 / X + AUG
  350 PSI = AUG + ALOG(X)
      RETURN
C---------------------------------------------------------------------
C     ERROR RETURN
C---------------------------------------------------------------------
  400 PSI = 0.0E0
      RETURN
C     ---------- LAST CARD OF PSI ---------
      END
