# psifn.f

          SUBROUTINE PSIFN (X, N, KODE, M, ANS, NZ, IERR)
    C***BEGIN PROLOGUE  PSIFN
    C***PURPOSE  Compute derivatives of the Psi function.
    C***LIBRARY   SLATEC
    C***CATEGORY  C7C
    C***TYPE      SINGLE PRECISION (PSIFN-S, DPSIFN-D)
    C***KEYWORDS  DERIVATIVES OF THE GAMMA FUNCTION, POLYGAMMA FUNCTION,
    C             PSI FUNCTION
    C***AUTHOR  Amos, D. E., (SNLA)
    C***DESCRIPTION
    C
    C         The following definitions are used in PSIFN:
    C
    C      Definition 1
    C         PSI(X) = d/dx (ln(GAMMA(X)), the first derivative of
    C                  the LOG GAMMA function.
    C      Definition 2
    C                     K   K
    C         PSI(K,X) = d /dx (PSI(X)), the K-th derivative of PSI(X).
    C   ___________________________________________________________________
    C       PSIFN computes a sequence of SCALED derivatives of
    C       the PSI function; i.e. for fixed X and M it computes
    C       the M-member sequence
    C
    C                  ((-1)**(K+1)/GAMMA(K+1))*PSI(K,X)
    C                    for K = N,...,N+M-1
    C
    C       where PSI(K,X) is as defined above.   For KODE=1, PSIFN returns
    C       the scaled derivatives as described.  KODE=2 is operative only
    C       when K=0 and in that case PSIFN returns -PSI(X) + LN(X).  That
    C       is, the logarithmic behavior for large X is removed when KODE=1
    C       and K=0.  When sums or differences of PSI functions are computed
    C       the logarithmic terms can be combined analytically and computed
    C       separately to help retain significant digits.
    C
    C         Note that CALL PSIFN(X,0,1,1,ANS) results in
    C                   ANS = -PSI(X)
    C
    C     Input
    C           X      - Argument, X .gt. 0.0E0
    C           N      - First member of the sequence, 0 .le. N .le. 100
    C                    N=0 gives ANS(1) = -PSI(X)       for KODE=1
    C                                       -PSI(X)+LN(X) for KODE=2
    C           KODE   - Selection parameter
    C                    KODE=1 returns scaled derivatives of the PSI
    C                    function.
    C                    KODE=2 returns scaled derivatives of the PSI
    C                    function EXCEPT when N=0. In this case,
    C                    ANS(1) = -PSI(X) + LN(X) is returned.
    C           M      - Number of members of the sequence, M .ge. 1
    C
    C    Output
    C           ANS    - A vector of length at least M whose first M
    C                    components contain the sequence of derivatives
    C                    scaled according to KODE.
    C           NZ     - Underflow flag
    C                    NZ.eq.0, A normal return
    C                    NZ.ne.0, Underflow, last NZ components of ANS are
    C                             set to zero, ANS(M-K+1)=0.0, K=1,...,NZ
    C           IERR   - Error flag
    C                    IERR=0, A normal return, computation completed
    C                    IERR=1, Input error,     no computation
    C                    IERR=2, Overflow,        X too small or N+M-1 too
    C                            large or both
    C                    IERR=3, Error,           N too large. Dimensioned
    C                            array TRMR(NMAX) is not large enough for N
    C
    C         The nominal computational accuracy is the maximum of unit
    C         roundoff (=R1MACH(4)) and 1.0E-18 since critical constants
    C         are given to only 18 digits.
    C
    C         DPSIFN is the Double Precision version of PSIFN.
    C
    C *Long Description:
    C
    C         The basic method of evaluation is the asymptotic expansion
    C         for large X.ge.XMIN followed by backward recursion on a two
    C         term recursion relation
    C
    C                  W(X+1) + X**(-N-1) = W(X).
    C
    C         This is supplemented by a series
    C
    C                  SUM( (X+K)**(-N-1) , K=0,1,2,... )
    C
    C         which converges rapidly for large N. Both XMIN and the
    C         number of terms of the series are calculated from the unit
    C         roundoff of the machine environment.
    C
    C***REFERENCES  Handbook of Mathematical Functions, National Bureau
    C                 of Standards Applied Mathematics Series 55, edited
    C                 by M. Abramowitz and I. A. Stegun, equations 6.3.5,
    C                 6.3.18, 6.4.6, 6.4.9 and 6.4.10, pp.258-260, 1964.
    C               D. E. Amos, A portable Fortran subroutine for
    C                 derivatives of the Psi function, Algorithm 610, ACM
    C                 Transactions on Mathematical Software 9, 4 (1983),
    C                 pp. 494-502.
    C***ROUTINES CALLED  I1MACH, R1MACH
    C***REVISION HISTORY  (YYMMDD)
    C   820601  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   920501  Reformatted the REFERENCES section.  (WRB)
    C***END PROLOGUE  PSIFN