relibc_openlibm/slatec/cpsi.f

*DECK CPSI
      COMPLEX FUNCTION CPSI (ZIN)
C***BEGIN PROLOGUE  CPSI
C***PURPOSE  Compute the Psi (or Digamma) function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7C
C***TYPE      COMPLEX (PSI-S, DPSI-D, CPSI-C)
C***KEYWORDS  DIGAMMA FUNCTION, FNLIB, PSI FUNCTION, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
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
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CCOT, R1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   780501  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900727  Added EXTERNAL statement.  (WRB)
C***END PROLOGUE  CPSI
      COMPLEX ZIN, Z, Z2INV, CORR, CCOT
      DIMENSION BERN(13)
      LOGICAL FIRST
      EXTERNAL CCOT
      SAVE BERN, PI, NTERM, BOUND, DXREL, RMIN, RBIG, FIRST
      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 FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  CPSI
      IF (FIRST) THEN
         NTERM = -0.30*LOG(R1MACH(3))
C MAYBE BOUND = N*(0.1*EPS)**(-1/(2*N-1)) / (PI*EXP(1))
         BOUND = 0.1171*NTERM*(0.1*R1MACH(3))**(-1.0/(2*NTERM-1))
         DXREL = SQRT(R1MACH(4))
         RMIN = EXP (MAX (LOG(R1MACH(1)), -LOG(R1MACH(2))) + 0.011 )
         RBIG = 1.0/R1MACH(3)
      ENDIF
      FIRST = .FALSE.
C
      Z = ZIN
      X = REAL(Z)
      Y = AIMAG(Z)
      IF (Y.LT.0.0) Z = CONJG(Z)
C
      CORR = (0.0, 0.0)
      CABSZ = ABS(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, ABS(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 ABS(Z) SMALL.
C
 20   IF (CABSZ .LT. RMIN) CALL XERMSG ('SLATEC', 'CPSI',
     +   'CPSI CALLED WITH Z SO NEAR 0 THAT CPSI OVERFLOWS', 2, 2)
C
      IF (X.GE.(-0.5) .OR. ABS(Y).GT.DXREL) GO TO 30
      IF (ABS((Z-AINT(X-0.5))/X) .LT. DXREL) CALL XERMSG ('SLATEC',
     +   'CPSI',
     +   'ANSWER LT HALF PRECISION BECAUSE Z TOO NEAR NEGATIVE INTEGER',
     +   1, 1)
      IF (Y .EQ. 0.0 .AND. X .EQ. AINT(X)) CALL XERMSG ('SLATEC',
     +   'CPSI', 'Z IS A NEGATIVE INTEGER', 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