--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/cruft/slatec-fn/psifn.f	Sun May 03 22:52:07 2015 +0100
@@ -0,0 +1,368 @@
+*DECK PSIFN
+      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
+      INTEGER I, IERR, J, K, KODE, M, MM, MX, N, NMAX, NN, NP, NX, NZ
+      INTEGER I1MACH
+      REAL ANS, ARG, B, DEN, ELIM, EPS, FLN, FN, FNP, FNS, 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
+      REAL R1MACH
+      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.00000000000000000E+00,
+     * -5.00000000000000000E-01,1.66666666666666667E-01,
+     * -3.33333333333333333E-02,2.38095238095238095E-02,
+     * -3.33333333333333333E-02,7.57575757575757576E-02,
+     * -2.53113553113553114E-01,1.16666666666666667E+00,
+     * -7.09215686274509804E+00,5.49711779448621554E+01,
+     * -5.29124242424242424E+02,6.19212318840579710E+03,
+     * -8.65802531135531136E+04,1.42551716666666667E+06,
+     * -2.72982310678160920E+07,6.01580873900642368E+08,
+     * -1.51163157670921569E+10,4.29614643061166667E+11,
+     * -1.37116552050883328E+13,4.88332318973593167E+14,
+     * -1.92965793419400681E+16/
+C
+C***FIRST EXECUTABLE STATEMENT  PSIFN
+      IERR = 0
+      NZ=0
+      IF (X.LE.0.0E0) 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(12),I1MACH(13))
+      R1M5 = R1MACH(5)
+      R1M4 = R1MACH(4)*0.5E0
+      WDTOL = MAX(R1M4,0.5E-18)
+C-----------------------------------------------------------------------
+C     ELIM = APPROXIMATE EXPONENTIAL OVER AND UNDERFLOW LIMIT
+C-----------------------------------------------------------------------
+      ELIM = 2.302E0*(NX*R1M5-3.0E0)
+      XLN = LOG(X)
+   41 CONTINUE
+      NN = N + MM - 1
+      FN = NN
+      FNP = FN + 1.0E0
+      T = FNP*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(11)
+      RLN = MIN(RLN,18.06E0)
+      FLN = MAX(RLN,3.0E0) - 3.0E0
+      YINT = 3.50E0 + 0.40E0*FLN
+      SLOPE = 0.21E0 + FLN*(0.0006038E0*FLN+0.008677E0)
+      XM = YINT + SLOPE*FN
+      MX = INT(XM) + 1
+      XMIN = MX
+      IF (N.EQ.0) GO TO 50
+      XM = -2.302E0*RLN - MIN(0.0E0,XLN)
+      FNS = N
+      ARG = XM/FNS
+      ARG = MIN(0.0E0,ARG)
+      EPS = EXP(ARG)
+      XM = 1.0E0 - EPS
+      IF (ABS(ARG).LT.1.0E-3) XM = -ARG
+      FLN = X*XM/EPS
+      XM = XMIN - X
+      IF (XM.GT.7.0E0 .AND. FLN.LT.15.0E0) GO TO 200
+   50 CONTINUE
+      XDMY = X
+      XDMLN = XLN
+      XINC = 0.0E0
+      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.5E0/XDMY
+      T1 = TT
+      TST = WDTOL*TT
+      IF (NN.NE.0) T1 = TT + 1.0E0/FN
+      RXSQ = 1.0E0/(XDMY*XDMY)
+      TA = 0.5E0*RXSQ
+      T = FNP*TA
+      S = T*B(3)
+      IF (ABS(S).LT.TST) GO TO 80
+      TK = 2.0E0
+      DO 70 K=4,22
+        T = T*((TK+FN+1.0E0)/(TK+1.0E0))*((TK+FN)/(TK+2.0E0))*RXSQ
+        TRM(K) = T*B(K)
+        IF (ABS(TRM(K)).LT.TST) GO TO 80
+        S = S + TRM(K)
+        TK = TK + 2.0E0
+   70 CONTINUE
+   80 CONTINUE
+      S = (S+T1)*TSS
+      IF (XINC.EQ.0.0E0) 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.0E0
+      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.0E0
+        FX = X + XM
+   90 CONTINUE
+  100 CONTINUE
+      ANS(MM) = S
+      IF (FN.EQ.0.0E0) GO TO 180
+C-----------------------------------------------------------------------
+C     GENERATE LOWER DERIVATIVES, J.LT.N+MM-1
+C-----------------------------------------------------------------------
+      IF (MM.EQ.1) RETURN
+      DO 150 J=2,MM
+        FNP = FN
+        FN = FN - 1.0E0
+        TSS = TSS*XDMY
+        T1 = TT
+        IF (FN.NE.0.0E0) T1 = TT + 1.0E0/FN
+        T = FNP*TA
+        S = T*B(3)
+        IF (ABS(S).LT.TST) GO TO 120
+        TK = 3.0E0 + FNP
+        DO 110 K=4,22
+          TRM(K) = TRM(K)*FNP/TK
+          IF (ABS(TRM(K)).LT.TST) GO TO 120
+          S = S + TRM(K)
+          TK = TK + 2.0E0
+  110   CONTINUE
+  120   CONTINUE
+        S = (S+T1)*TSS
+        IF (XINC.EQ.0.0E0) GO TO 140
+        IF (FN.EQ.0.0E0) GO TO 160
+        XM = XINC - 1.0E0
+        FX = X + XM
+        DO 130 I=1,NX
+          TRMR(I) = TRMR(I)*FX
+          S = S + TRMR(I)
+          XM = XM - 1.0E0
+          FX = X + XM
+  130   CONTINUE
+  140   CONTINUE
+        MX = MM - J + 1
+        ANS(MX) = S
+        IF (FN.EQ.0.0E0) GO TO 180
+  150 CONTINUE
+      RETURN
+C-----------------------------------------------------------------------
+C     RECURSION FOR N = 0
+C-----------------------------------------------------------------------
+  160 CONTINUE
+      DO 170 I=1,NX
+        S = S + 1.0E0/(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 = (FNS+1.0E0)*XLN
+      T = EXP(-T1)
+      S = T
+      DEN = X
+      DO 210 I=1,NN
+        DEN = DEN + 1.0E0
+        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.0E0
+      DO 250 J=2,MM
+        T = T/X
+        S = T
+        TOLS = T*TOL
+        DEN = X
+        DO 230 I=1,NN
+          DEN = DEN + 1.0E0
+          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.0E0) GO TO 380
+      NZ=0
+      IERR=2
+      RETURN
+  380 CONTINUE
+      NZ=NZ+1
+      ANS(MM)=0.0E0
+      MM=MM-1
+      IF(MM.EQ.0) RETURN
+      GO TO 41
+  390 CONTINUE
+      IERR=3
+      NZ=0
+      RETURN
+      END