      COMPLEX FUNCTION CPSI(ZIN)
C***BEGIN PROLOGUE  CPSI
C***DATE WRITTEN   780501   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***REVISION HISTORY (YYMMDD)
C   000601  Changed CLOG to generic LOG
C***CATEGORY NO.  C7C
C***KEYWORDS  COMPLEX,DIGAMMA FUNCTION,PSI FUNCTION,SPECIAL FUNCTION
C***AUTHOR  FULLERTON, W., (LANL)
C***PURPOSE  Computes the Psi function of complex argument.
C***DESCRIPTION
C
C PSI(X) calculates the psi (or digamma) function of X.  PSI(X)
C is the logarithmic derivative of the gamma function of X.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CCOT,R1MACH,XERROR
C***END PROLOGUE  CPSI
      COMPLEX ZIN, Z, Z2INV, CORR,  CCOT
      DIMENSION BERN(13)
      DATA BERN( 1) /   .8333333333 3333333 E-1 /
      DATA BERN( 2) /  -.8333333333 3333333 E-2 /
      DATA BERN( 3) /   .3968253968 2539683 E-2 /
      DATA BERN( 4) /  -.4166666666 6666667 E-2 /
      DATA BERN( 5) /   .7575757575 7575758 E-2 /
      DATA BERN( 6) /  -.2109279609 2796093 E-1 /
      DATA BERN( 7) /   .8333333333 3333333 E-1 /
      DATA BERN( 8) /  -.4432598039 2156863 E0 /
      DATA BERN( 9) /   .3053954330 2701197 E1 /
      DATA BERN(10) /  -.2645621212 1212121 E2 /
      DATA BERN(11) /   .2814601449 2753623 E3 /
      DATA BERN(12) /  -.3454885393 7728938 E4 /
      DATA BERN(13) /   .5482758333 3333333 E5 /
      DATA PI / 3.141592653 589793 E0 /
      DATA NTERM, BOUND, DXREL, RMIN, RBIG / 0, 4*0.0 /
C***FIRST EXECUTABLE STATEMENT  CPSI
      IF (NTERM.NE.0) GO TO 10
      NTERM = -0.30*ALOG(R1MACH(3))
C MAYBE BOUND = N*(0.1*EPS)**(-1/(2*N-1)) / (PI*EXP(1))
      BOUND = 0.1171*FLOAT(NTERM) *
     1  (0.1*R1MACH(3))**(-1.0/(2.0*FLOAT(NTERM)-1.0))
      DXREL = SQRT(R1MACH(4))
      RMIN = EXP (AMAX1 (ALOG(R1MACH(1)), -ALOG(R1MACH(2))) + 0.011 )
      RBIG = 1.0/R1MACH(3)
C
 10   Z = ZIN
      X = REAL(Z)
      Y = AIMAG(Z)
      IF (Y.LT.0.0) Z = CONJG(Z)
C
      CORR = (0.0, 0.0)
      CABSZ = CABS(Z)
      IF (X.GE.0.0 .AND. CABSZ.GT.BOUND) GO TO 50
      IF (X.LT.0.0 .AND. ABS(Y).GT.BOUND) GO TO 50
C
      IF (CABSZ.LT.BOUND) GO TO 20
C
C USE THE REFLECTION FORMULA FOR REAL(Z) NEGATIVE, CABS(Z) LARGE, AND
C ABS(AIMAG(Y)) SMALL.
C
      CORR = -PI*CCOT(PI*Z)
      Z = 1.0 - Z
      GO TO 50
C
C USE THE RECURSION RELATION FOR CABS(Z) SMALL.
C
 20   IF (CABSZ.LT.RMIN) CALL XERROR ( 'CPSI    CPSI CALLED WITH Z SO NE
     1AR 0 THAT CPSI OVERFLOWS',      56, 2, 2)
C
      IF (X.GE.(-0.5) .OR. ABS(Y).GT.DXREL) GO TO 30
      IF (CABS((Z-AINT(X-0.5))/X).LT.DXREL) CALL XERROR ( 'CPSI    ANSWE
     1R LT HALF PRECISION BECAUSE Z TOO NEAR NEGATIVE INTEGER', 68, 1, 1
     2)
      IF (Y.EQ.0.0 .AND. X.EQ.AINT(X)) CALL XERROR ( 'CPSI    Z IS A NEG
     1ATIVE INTEGER', 31, 3, 2)
C
 30   N = SQRT(BOUND**2-Y**2) - X + 1.0
      DO 40 I=1,N
        CORR = CORR - 1.0/Z
        Z = Z + 1.0
 40   CONTINUE
C
C NOW EVALUATE THE ASYMPTOTIC SERIES FOR SUITABLY LARGE Z.
C
 50   IF (CABSZ.GT.RBIG) CPSI = LOG(Z) + CORR
      IF (CABSZ.GT.RBIG) GO TO 70
C
      CPSI = (0.0, 0.0)
      Z2INV = 1.0/Z**2
      DO 60 I=1,NTERM
        NDX = NTERM + 1 - I
        CPSI = BERN(NDX) + Z2INV*CPSI
 60   CONTINUE
      CPSI = LOG(Z) - 0.5/Z - CPSI*Z2INV + CORR
C
 70   IF (Y.LT.0.0) CPSI = CONJG(CPSI)
C
      RETURN
      END
