1*DECK DHKSEQ
2      SUBROUTINE DHKSEQ (X, M, H, IERR)
3C***BEGIN PROLOGUE  DHKSEQ
4C***SUBSIDIARY
5C***PURPOSE  Subsidiary to DBSKIN
6C***LIBRARY   SLATEC
7C***TYPE      DOUBLE PRECISION (HKSEQ-S, DHKSEQ-D)
8C***AUTHOR  Amos, D. E., (SNLA)
9C***DESCRIPTION
10C
11C   DHKSEQ is an adaptation of subroutine DPSIFN described in the
12C   reference below.  DHKSEQ generates the sequence
13C   H(K,X) = (-X)**(K+1)*(PSI(K,X) PSI(K,X+0.5))/GAMMA(K+1), for
14C            K=0,...,M.
15C
16C***SEE ALSO  DBSKIN
17C***REFERENCES  D. E. Amos, A portable Fortran subroutine for
18C                 derivatives of the Psi function, Algorithm 610, ACM
19C                 Transactions on Mathematical Software 9, 4 (1983),
20C                 pp. 494-502.
21C***ROUTINES CALLED  D1MACH, I1MACH
22C***REVISION HISTORY  (YYMMDD)
23C   820601  DATE WRITTEN
24C   890531  Changed all specific intrinsics to generic.  (WRB)
25C   890911  Removed unnecessary intrinsics.  (WRB)
26C   891214  Prologue converted to Version 4.0 format.  (BAB)
27C   900328  Added TYPE section.  (WRB)
28C   910722  Updated AUTHOR section.  (ALS)
29C   920528  DESCRIPTION and REFERENCES sections revised.  (WRB)
30C***END PROLOGUE  DHKSEQ
31      INTEGER I, IERR, J, K, M, MX, NX
32      INTEGER I1MACH
33      DOUBLE PRECISION B, FK, FLN, FN, FNP, H, HRX, RLN, RXSQ, R1M5, S,
34     * SLOPE, T, TK, TRM, TRMH, TRMR, TST, U, V, WDTOL, X, XDMY, XH,
35     * XINC, XM, XMIN, YINT
36      DOUBLE PRECISION D1MACH
37      DIMENSION B(22), TRM(22), TRMR(25), TRMH(25), U(25), V(25), H(*)
38      SAVE B
39C-----------------------------------------------------------------------
40C             SCALED BERNOULLI NUMBERS 2.0*B(2K)*(1-2**(-2K))
41C-----------------------------------------------------------------------
42      DATA B(1), B(2), B(3), B(4), B(5), B(6), B(7), B(8), B(9), B(10),
43     * B(11), B(12), B(13), B(14), B(15), B(16), B(17), B(18), B(19),
44     * B(20), B(21), B(22) /1.00000000000000000D+00,
45     * -5.00000000000000000D-01,2.50000000000000000D-01,
46     * -6.25000000000000000D-02,4.68750000000000000D-02,
47     * -6.64062500000000000D-02,1.51367187500000000D-01,
48     * -5.06103515625000000D-01,2.33319091796875000D+00,
49     * -1.41840972900390625D+01,1.09941936492919922D+02,
50     * -1.05824747562408447D+03,1.23842434241771698D+04,
51     * -1.73160495905935764D+05,2.85103429084961116D+06,
52     * -5.45964619322445132D+07,1.20316174668075304D+09,
53     * -3.02326315271452307D+10,8.59229286072319606D+11,
54     * -2.74233104097776039D+13,9.76664637943633248D+14,
55     * -3.85931586838450360D+16/
56C
57C***FIRST EXECUTABLE STATEMENT  DHKSEQ
58      IERR=0
59      WDTOL = MAX(D1MACH(4),1.0D-18)
60      FN = M - 1
61      FNP = FN + 1.0D0
62C-----------------------------------------------------------------------
63C     COMPUTE XMIN
64C-----------------------------------------------------------------------
65      R1M5 = D1MACH(5)
66      RLN = R1M5*I1MACH(14)
67      RLN = MIN(RLN,18.06D0)
68      FLN = MAX(RLN,3.0D0) - 3.0D0
69      YINT = 3.50D0 + 0.40D0*FLN
70      SLOPE = 0.21D0 + FLN*(0.0006038D0*FLN+0.008677D0)
71      XM = YINT + SLOPE*FN
72      MX = INT(XM) + 1
73      XMIN = MX
74C-----------------------------------------------------------------------
75C     GENERATE H(M-1,XDMY)*XDMY**(M) BY THE ASYMPTOTIC EXPANSION
76C-----------------------------------------------------------------------
77      XDMY = X
78      XINC = 0.0D0
79      IF (X.GE.XMIN) GO TO 10
80      NX = INT(X)
81      XINC = XMIN - NX
82      XDMY = X + XINC
83   10 CONTINUE
84      RXSQ = 1.0D0/(XDMY*XDMY)
85      HRX = 0.5D0/XDMY
86      TST = 0.5D0*WDTOL
87      T = FNP*HRX
88C-----------------------------------------------------------------------
89C     INITIALIZE COEFFICIENT ARRAY
90C-----------------------------------------------------------------------
91      S = T*B(3)
92      IF (ABS(S).LT.TST) GO TO 30
93      TK = 2.0D0
94      DO 20 K=4,22
95        T = T*((TK+FN+1.0D0)/(TK+1.0D0))*((TK+FN)/(TK+2.0D0))*RXSQ
96        TRM(K) = T*B(K)
97        IF (ABS(TRM(K)).LT.TST) GO TO 30
98        S = S + TRM(K)
99        TK = TK + 2.0D0
100   20 CONTINUE
101      GO TO 110
102   30 CONTINUE
103      H(M) = S + 0.5D0
104      IF (M.EQ.1) GO TO 70
105C-----------------------------------------------------------------------
106C     GENERATE LOWER DERIVATIVES, I.LT.M-1
107C-----------------------------------------------------------------------
108      DO 60 I=2,M
109        FNP = FN
110        FN = FN - 1.0D0
111        S = FNP*HRX*B(3)
112        IF (ABS(S).LT.TST) GO TO 50
113        FK = FNP + 3.0D0
114        DO 40 K=4,22
115          TRM(K) = TRM(K)*FNP/FK
116          IF (ABS(TRM(K)).LT.TST) GO TO 50
117          S = S + TRM(K)
118          FK = FK + 2.0D0
119   40   CONTINUE
120        GO TO 110
121   50   CONTINUE
122        MX = M - I + 1
123        H(MX) = S + 0.5D0
124   60 CONTINUE
125   70 CONTINUE
126      IF (XINC.EQ.0.0D0) RETURN
127C-----------------------------------------------------------------------
128C     RECUR BACKWARD FROM XDMY TO X
129C-----------------------------------------------------------------------
130      XH = X + 0.5D0
131      S = 0.0D0
132      NX = INT(XINC)
133      DO 80 I=1,NX
134        TRMR(I) = X/(X+NX-I)
135        U(I) = TRMR(I)
136        TRMH(I) = X/(XH+NX-I)
137        V(I) = TRMH(I)
138        S = S + U(I) - V(I)
139   80 CONTINUE
140      MX = NX + 1
141      TRMR(MX) = X/XDMY
142      U(MX) = TRMR(MX)
143      H(1) = H(1)*TRMR(MX) + S
144      IF (M.EQ.1) RETURN
145      DO 100 J=2,M
146        S = 0.0D0
147        DO 90 I=1,NX
148          TRMR(I) = TRMR(I)*U(I)
149          TRMH(I) = TRMH(I)*V(I)
150          S = S + TRMR(I) - TRMH(I)
151   90   CONTINUE
152        TRMR(MX) = TRMR(MX)*U(MX)
153        H(J) = H(J)*TRMR(MX) + S
154  100 CONTINUE
155      RETURN
156  110 CONTINUE
157      IERR=2
158      RETURN
159      END
160