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