dpsifn.f

Go to the documentation of this file.

*DECK DPSIFN
      SUBROUTINE dpsifn (X, N, KODE, M, ANS, NZ, IERR)
C***BEGIN PROLOGUE  DPSIFN
C***PURPOSE  Compute derivatives of the Psi function.
C***LIBRARY   SLATEC
C***CATEGORY  C7C
C***TYPE      DOUBLE 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 DPSIFN:
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      DPSIFN 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, DPSIFN returns
C      the scaled derivatives as described.  KODE=2 is operative only
C      when K=0 and in that case DPSIFN returns -PSI(X) + LN(X).  That
C      is, the logarithmic behavior for large X is removed when KODE=2
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 DPSIFN(X,0,1,1,ANS) results in
C                   ANS = -PSI(X)
C
C     Input      X is DOUBLE PRECISION
C           X      - Argument, X .gt. 0.0D0
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     ANS is DOUBLE PRECISION
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 (=D1MACH(4)) and 1.0D-18 since critical constants
C         are given to only 18 digits.
C
C         PSIFN is the single precision version of DPSIFN.
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  D1MACH, I1MACH
C***REVISION HISTORY  (YYMMDD)
C   820601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890911  Removed unnecessary intrinsics.  (WRB)
C   891006  Cosmetic changes to prologue.  (WRB)
C   891006  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  DPSIFN
      INTEGER I, IERR, J, K, KODE, M, MM, MX, N, NMAX, NN, NP, NX, NZ,
     *  FN
      INTEGER I1MACH
      DOUBLE PRECISION ANS, ARG, B, DEN, ELIM, EPS, FLN,
     * FX, RLN, RXSQ, R1M4, R1M5, S, SLOPE, T, TA, TK, TOL, TOLS, TRM,
     * TRMR, TSS, TST, TT, T1, T2, WDTOL, X, XDMLN, XDMY, XINC, XLN,
     * XM, XMIN, XQ, YINT
      DOUBLE PRECISION D1MACH
      dimension b(22), trm(22), trmr(100), ans(*)
      SAVE nmax, b
      DATA nmax /100/
C-----------------------------------------------------------------------
C             BERNOULLI NUMBERS
C-----------------------------------------------------------------------
      DATA b(1), b(2), b(3), b(4), b(5), b(6), b(7), b(8), b(9), b(10),
     * b(11), b(12), b(13), b(14), b(15), b(16), b(17), b(18), b(19),
     * b(20), b(21), b(22) /1.00000000000000000d+00,
     * -5.00000000000000000d-01,1.66666666666666667d-01,
     * -3.33333333333333333d-02,2.38095238095238095d-02,
     * -3.33333333333333333d-02,7.57575757575757576d-02,
     * -2.53113553113553114d-01,1.16666666666666667d+00,
     * -7.09215686274509804d+00,5.49711779448621554d+01,
     * -5.29124242424242424d+02,6.19212318840579710d+03,
     * -8.65802531135531136d+04,1.42551716666666667d+06,
     * -2.72982310678160920d+07,6.01580873900642368d+08,
     * -1.51163157670921569d+10,4.29614643061166667d+11,
     * -1.37116552050883328d+13,4.88332318973593167d+14,
     * -1.92965793419400681d+16/
C
C***FIRST EXECUTABLE STATEMENT  DPSIFN
      ierr = 0
      nz=0
      IF (x.LE.0.0d0) ierr=1
      IF (n.LT.0) ierr=1
      IF (kode.LT.1 .OR. kode.GT.2) ierr=1
      IF (m.LT.1) ierr=1
      IF (ierr.NE.0) RETURN
      mm=m
      nx = min(-i1mach(15),i1mach(16))
      r1m5 = d1mach(5)
      r1m4 = d1mach(4)*0.5d0
      wdtol = max(r1m4,0.5d-18)
C-----------------------------------------------------------------------
C     ELIM = APPROXIMATE EXPONENTIAL OVER AND UNDERFLOW LIMIT
C-----------------------------------------------------------------------
      elim = 2.302d0*(nx*r1m5-3.0d0)
      xln = log(x)
   41 CONTINUE
      nn = n + mm - 1
      fn = nn
      t = (fn+1)*xln
C-----------------------------------------------------------------------
C     OVERFLOW AND UNDERFLOW TEST FOR SMALL AND LARGE X
C-----------------------------------------------------------------------
      IF (abs(t).GT.elim) GO TO 290
      IF (x.LT.wdtol) GO TO 260
C-----------------------------------------------------------------------
C     COMPUTE XMIN AND THE NUMBER OF TERMS OF THE SERIES, FLN+1
C-----------------------------------------------------------------------
      rln = r1m5*i1mach(14)
      rln = min(rln,18.06d0)
      fln = max(rln,3.0d0) - 3.0d0
      yint = 3.50d0 + 0.40d0*fln
      slope = 0.21d0 + fln*(0.0006038d0*fln+0.008677d0)
      xm = yint + slope*fn
      mx = int(xm) + 1
      xmin = mx
      IF (n.EQ.0) GO TO 50
      xm = -2.302d0*rln - min(0.0d0,xln)
      arg = xm/n
      arg = min(0.0d0,arg)
      eps = exp(arg)
      xm = 1.0d0 - eps
      IF (abs(arg).LT.1.0d-3) xm = -arg
      fln = x*xm/eps
      xm = xmin - x
      IF (xm.GT.7.0d0 .AND. fln.LT.15.0d0) GO TO 200
   50 CONTINUE
      xdmy = x
      xdmln = xln
      xinc = 0.0d0
      IF (x.GE.xmin) GO TO 60
      nx = int(x)
      xinc = xmin - nx
      xdmy = x + xinc
      xdmln = log(xdmy)
   60 CONTINUE
C-----------------------------------------------------------------------
C     GENERATE W(N+MM-1,X) BY THE ASYMPTOTIC EXPANSION
C-----------------------------------------------------------------------
      t = fn*xdmln
      t1 = xdmln + xdmln
      t2 = t + xdmln
      tk = max(abs(t),abs(t1),abs(t2))
      IF (tk.GT.elim) GO TO 380
      tss = exp(-t)
      tt = 0.5d0/xdmy
      t1 = tt
      tst = wdtol*tt
      IF (nn.NE.0) t1 = tt + 1.0d0/fn
      rxsq = 1.0d0/(xdmy*xdmy)
      ta = 0.5d0*rxsq
      t = (fn+1)*ta
      s = t*b(3)
      IF (abs(s).LT.tst) GO TO 80
      tk = 2.0d0
      DO 70 k=4,22
        t = t*((tk+fn+1)/(tk+1.0d0))*((tk+fn)/(tk+2.0d0))*rxsq
        trm(k) = t*b(k)
        IF (abs(trm(k)).LT.tst) GO TO 80
        s = s + trm(k)
        tk = tk + 2.0d0
   70 CONTINUE
   80 CONTINUE
      s = (s+t1)*tss
      IF (xinc.EQ.0.0d0) GO TO 100
C-----------------------------------------------------------------------
C     BACKWARD RECUR FROM XDMY TO X
C-----------------------------------------------------------------------
      nx = int(xinc)
      np = nn + 1
      IF (nx.GT.nmax) GO TO 390
      IF (nn.EQ.0) GO TO 160
      xm = xinc - 1.0d0
      fx = x + xm
C-----------------------------------------------------------------------
C     THIS LOOP SHOULD NOT BE CHANGED. FX IS ACCURATE WHEN X IS SMALL
C-----------------------------------------------------------------------
      DO 90 i=1,nx
        trmr(i) = fx**(-np)
        s = s + trmr(i)
        xm = xm - 1.0d0
        fx = x + xm
   90 CONTINUE
  100 CONTINUE
      ans(mm) = s
      IF (fn.EQ.0) GO TO 180
C-----------------------------------------------------------------------
C     GENERATE LOWER DERIVATIVES, J.LT.N+MM-1
C-----------------------------------------------------------------------
      IF (mm.EQ.1) RETURN
      DO 150 j=2,mm
        fn = fn - 1
        tss = tss*xdmy
        t1 = tt
        IF (fn.NE.0) t1 = tt + 1.0d0/fn
        t = (fn+1)*ta
        s = t*b(3)
        IF (abs(s).LT.tst) GO TO 120
        tk = 4 + fn
        DO 110 k=4,22
          trm(k) = trm(k)*(fn+1)/tk
          IF (abs(trm(k)).LT.tst) GO TO 120
          s = s + trm(k)
          tk = tk + 2.0d0
  110   CONTINUE
  120   CONTINUE
        s = (s+t1)*tss
        IF (xinc.EQ.0.0d0) GO TO 140
        IF (fn.EQ.0) GO TO 160
        xm = xinc - 1.0d0
        fx = x + xm
        DO 130 i=1,nx
          trmr(i) = trmr(i)*fx
          s = s + trmr(i)
          xm = xm - 1.0d0
          fx = x + xm
  130   CONTINUE
  140   CONTINUE
        mx = mm - j + 1
        ans(mx) = s
        IF (fn.EQ.0) GO TO 180
  150 CONTINUE
      RETURN
C-----------------------------------------------------------------------
C     RECURSION FOR N = 0
C-----------------------------------------------------------------------
  160 CONTINUE
      DO 170 i=1,nx
        s = s + 1.0d0/(x+nx-i)
  170 CONTINUE
  180 CONTINUE
      IF (kode.EQ.2) GO TO 190
      ans(1) = s - xdmln
      RETURN
  190 CONTINUE
      IF (xdmy.EQ.x) RETURN
      xq = xdmy/x
      ans(1) = s - log(xq)
      RETURN
C-----------------------------------------------------------------------
C     COMPUTE BY SERIES (X+K)**(-(N+1)) , K=0,1,2,...
C-----------------------------------------------------------------------
  200 CONTINUE
      nn = int(fln) + 1
      np = n + 1
      t1 = (n+1)*xln
      t = exp(-t1)
      s = t
      den = x
      DO 210 i=1,nn
        den = den + 1.0d0
        trm(i) = den**(-np)
        s = s + trm(i)
  210 CONTINUE
      ans(1) = s
      IF (n.NE.0) GO TO 220
      IF (kode.EQ.2) ans(1) = s + xln
  220 CONTINUE
      IF (mm.EQ.1) RETURN
C-----------------------------------------------------------------------
C     GENERATE HIGHER DERIVATIVES, J.GT.N
C-----------------------------------------------------------------------
      tol = wdtol/5.0d0
      DO 250 j=2,mm
        t = t/x
        s = t
        tols = t*tol
        den = x
        DO 230 i=1,nn
          den = den + 1.0d0
          trm(i) = trm(i)/den
          s = s + trm(i)
          IF (trm(i).LT.tols) GO TO 240
  230   CONTINUE
  240   CONTINUE
        ans(j) = s
  250 CONTINUE
      RETURN
C-----------------------------------------------------------------------
C     SMALL X.LT.UNIT ROUND OFF
C-----------------------------------------------------------------------
  260 CONTINUE
      ans(1) = x**(-n-1)
      IF (mm.EQ.1) GO TO 280
      k = 1
      DO 270 i=2,mm
        ans(k+1) = ans(k)/x
        k = k + 1
  270 CONTINUE
  280 CONTINUE
      IF (n.NE.0) RETURN
      IF (kode.EQ.2) ans(1) = ans(1) + xln
      RETURN
  290 CONTINUE
      IF (t.GT.0.0d0) GO TO 380
      nz=0
      ierr=2
      RETURN
  380 CONTINUE
      nz=nz+1
      ans(mm)=0.0d0
      mm=mm-1
      IF (mm.EQ.0) RETURN
      GO TO 41
  390 CONTINUE
      nz=0
      ierr=3
      RETURN
      END