GNU Octave  8.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dpsifn.f
Go to the documentation of this file.
1 *DECK DPSIFN
2  SUBROUTINE dpsifn (X, N, KODE, M, ANS, NZ, IERR)
3 C***BEGIN PROLOGUE DPSIFN
4 C***PURPOSE Compute derivatives of the Psi function.
5 C***LIBRARY SLATEC
6 C***CATEGORY C7C
7 C***TYPE DOUBLE PRECISION (PSIFN-S, DPSIFN-D)
8 C***KEYWORDS DERIVATIVES OF THE GAMMA FUNCTION, POLYGAMMA FUNCTION,
9 C PSI FUNCTION
10 C***AUTHOR Amos, D. E., (SNLA)
11 C***DESCRIPTION
12 C
13 C The following definitions are used in DPSIFN:
14 C
15 C Definition 1
16 C PSI(X) = d/dx (ln(GAMMA(X)), the first derivative of
17 C the log GAMMA function.
18 C Definition 2
19 C K K
20 C PSI(K,X) = d /dx (PSI(X)), the K-th derivative of PSI(X).
21 C ___________________________________________________________________
22 C DPSIFN computes a sequence of SCALED derivatives of
23 C the PSI function; i.e. for fixed X and M it computes
24 C the M-member sequence
25 C
26 C ((-1)**(K+1)/GAMMA(K+1))*PSI(K,X)
27 C for K = N,...,N+M-1
28 C
29 C where PSI(K,X) is as defined above. For KODE=1, DPSIFN returns
30 C the scaled derivatives as described. KODE=2 is operative only
31 C when K=0 and in that case DPSIFN returns -PSI(X) + LN(X). That
32 C is, the logarithmic behavior for large X is removed when KODE=2
33 C and K=0. When sums or differences of PSI functions are computed
34 C the logarithmic terms can be combined analytically and computed
35 C separately to help retain significant digits.
36 C
37 C Note that CALL DPSIFN(X,0,1,1,ANS) results in
38 C ANS = -PSI(X)
39 C
40 C Input X is DOUBLE PRECISION
41 C X - Argument, X .gt. 0.0D0
42 C N - First member of the sequence, 0 .le. N .le. 100
43 C N=0 gives ANS(1) = -PSI(X) for KODE=1
44 C -PSI(X)+LN(X) for KODE=2
45 C KODE - Selection parameter
46 C KODE=1 returns scaled derivatives of the PSI
47 C function.
48 C KODE=2 returns scaled derivatives of the PSI
49 C function EXCEPT when N=0. In this case,
50 C ANS(1) = -PSI(X) + LN(X) is returned.
51 C M - Number of members of the sequence, M.ge.1
52 C
53 C Output ANS is DOUBLE PRECISION
54 C ANS - A vector of length at least M whose first M
55 C components contain the sequence of derivatives
56 C scaled according to KODE.
57 C NZ - Underflow flag
58 C NZ.eq.0, A normal return
59 C NZ.ne.0, Underflow, last NZ components of ANS are
60 C set to zero, ANS(M-K+1)=0.0, K=1,...,NZ
61 C IERR - Error flag
62 C IERR=0, A normal return, computation completed
63 C IERR=1, Input error, no computation
64 C IERR=2, Overflow, X too small or N+M-1 too
65 C large or both
66 C IERR=3, Error, N too large. Dimensioned
67 C array TRMR(NMAX) is not large enough for N
68 C
69 C The nominal computational accuracy is the maximum of unit
70 C roundoff (=D1MACH(4)) and 1.0D-18 since critical constants
71 C are given to only 18 digits.
72 C
73 C PSIFN is the single precision version of DPSIFN.
74 C
75 C *Long Description:
76 C
77 C The basic method of evaluation is the asymptotic expansion
78 C for large X.ge.XMIN followed by backward recursion on a two
79 C term recursion relation
80 C
81 C W(X+1) + X**(-N-1) = W(X).
82 C
83 C This is supplemented by a series
84 C
85 C SUM( (X+K)**(-N-1) , K=0,1,2,... )
86 C
87 C which converges rapidly for large N. Both XMIN and the
88 C number of terms of the series are calculated from the unit
89 C roundoff of the machine environment.
90 C
91 C***REFERENCES Handbook of Mathematical Functions, National Bureau
92 C of Standards Applied Mathematics Series 55, edited
93 C by M. Abramowitz and I. A. Stegun, equations 6.3.5,
94 C 6.3.18, 6.4.6, 6.4.9 and 6.4.10, pp.258-260, 1964.
95 C D. E. Amos, A portable Fortran subroutine for
96 C derivatives of the Psi function, Algorithm 610, ACM
97 C Transactions on Mathematical Software 9, 4 (1983),
98 C pp. 494-502.
99 C***ROUTINES CALLED D1MACH, I1MACH
100 C***REVISION HISTORY (YYMMDD)
101 C 820601 DATE WRITTEN
102 C 890531 Changed all specific intrinsics to generic. (WRB)
103 C 890911 Removed unnecessary intrinsics. (WRB)
104 C 891006 Cosmetic changes to prologue. (WRB)
105 C 891006 REVISION DATE from Version 3.2
106 C 891214 Prologue converted to Version 4.0 format. (BAB)
107 C 920501 Reformatted the REFERENCES section. (WRB)
108 C***END PROLOGUE DPSIFN
109  INTEGER I, IERR, J, K, KODE, M, MM, MX, N, NMAX, NN, NP, NX, NZ,
110  * FN
111  INTEGER I1MACH
112  DOUBLE PRECISION ANS, ARG, B, DEN, ELIM, EPS, FLN,
113  * FX, RLN, RXSQ, R1M4, R1M5, S, SLOPE, T, TA, TK, TOL, TOLS, TRM,
114  * TRMR, TSS, TST, TT, T1, T2, WDTOL, X, XDMLN, XDMY, XINC, XLN,
115  * XM, XMIN, XQ, YINT
116  DOUBLE PRECISION D1MACH
117  dimension b(22), trm(22), trmr(100), ans(*)
118  SAVE nmax, b
119  DATA nmax /100/
120 C-----------------------------------------------------------------------
121 C BERNOULLI NUMBERS
122 C-----------------------------------------------------------------------
123  DATA b(1), b(2), b(3), b(4), b(5), b(6), b(7), b(8), b(9), b(10),
124  * b(11), b(12), b(13), b(14), b(15), b(16), b(17), b(18), b(19),
125  * b(20), b(21), b(22) /1.00000000000000000d+00,
126  * -5.00000000000000000d-01,1.66666666666666667d-01,
127  * -3.33333333333333333d-02,2.38095238095238095d-02,
128  * -3.33333333333333333d-02,7.57575757575757576d-02,
129  * -2.53113553113553114d-01,1.16666666666666667d+00,
130  * -7.09215686274509804d+00,5.49711779448621554d+01,
131  * -5.29124242424242424d+02,6.19212318840579710d+03,
132  * -8.65802531135531136d+04,1.42551716666666667d+06,
133  * -2.72982310678160920d+07,6.01580873900642368d+08,
134  * -1.51163157670921569d+10,4.29614643061166667d+11,
135  * -1.37116552050883328d+13,4.88332318973593167d+14,
136  * -1.92965793419400681d+16/
137 C
138 C***FIRST EXECUTABLE STATEMENT DPSIFN
139  ierr = 0
140  nz=0
141  IF (x.LE.0.0d0) ierr=1
142  IF (n.LT.0) ierr=1
143  IF (kode.LT.1 .OR. kode.GT.2) ierr=1
144  IF (m.LT.1) ierr=1
145  IF (ierr.NE.0) RETURN
146  mm=m
147  nx = min(-i1mach(15),i1mach(16))
148  r1m5 = d1mach(5)
149  r1m4 = d1mach(4)*0.5d0
150  wdtol = max(r1m4,0.5d-18)
151 C-----------------------------------------------------------------------
152 C ELIM = APPROXIMATE EXPONENTIAL OVER AND UNDERFLOW LIMIT
153 C-----------------------------------------------------------------------
154  elim = 2.302d0*(nx*r1m5-3.0d0)
155  xln = log(x)
156  41 CONTINUE
157  nn = n + mm - 1
158  fn = nn
159  t = (fn+1)*xln
160 C-----------------------------------------------------------------------
161 C OVERFLOW AND UNDERFLOW TEST FOR SMALL AND LARGE X
162 C-----------------------------------------------------------------------
163  IF (abs(t).GT.elim) GO TO 290
164  IF (x.LT.wdtol) GO TO 260
165 C-----------------------------------------------------------------------
166 C COMPUTE XMIN AND THE NUMBER OF TERMS OF THE SERIES, FLN+1
167 C-----------------------------------------------------------------------
168  rln = r1m5*i1mach(14)
169  rln = min(rln,18.06d0)
170  fln = max(rln,3.0d0) - 3.0d0
171  yint = 3.50d0 + 0.40d0*fln
172  slope = 0.21d0 + fln*(0.0006038d0*fln+0.008677d0)
173  xm = yint + slope*fn
174  mx = int(xm) + 1
175  xmin = mx
176  IF (n.EQ.0) GO TO 50
177  xm = -2.302d0*rln - min(0.0d0,xln)
178  arg = xm/n
179  arg = min(0.0d0,arg)
180  eps = exp(arg)
181  xm = 1.0d0 - eps
182  IF (abs(arg).LT.1.0d-3) xm = -arg
183  fln = x*xm/eps
184  xm = xmin - x
185  IF (xm.GT.7.0d0 .AND. fln.LT.15.0d0) GO TO 200
186  50 CONTINUE
187  xdmy = x
188  xdmln = xln
189  xinc = 0.0d0
190  IF (x.GE.xmin) GO TO 60
191  nx = int(x)
192  xinc = xmin - nx
193  xdmy = x + xinc
194  xdmln = log(xdmy)
195  60 CONTINUE
196 C-----------------------------------------------------------------------
197 C GENERATE W(N+MM-1,X) BY THE ASYMPTOTIC EXPANSION
198 C-----------------------------------------------------------------------
199  t = fn*xdmln
200  t1 = xdmln + xdmln
201  t2 = t + xdmln
202  tk = max(abs(t),abs(t1),abs(t2))
203  IF (tk.GT.elim) GO TO 380
204  tss = exp(-t)
205  tt = 0.5d0/xdmy
206  t1 = tt
207  tst = wdtol*tt
208  IF (nn.NE.0) t1 = tt + 1.0d0/fn
209  rxsq = 1.0d0/(xdmy*xdmy)
210  ta = 0.5d0*rxsq
211  t = (fn+1)*ta
212  s = t*b(3)
213  IF (abs(s).LT.tst) GO TO 80
214  tk = 2.0d0
215  DO 70 k=4,22
216  t = t*((tk+fn+1)/(tk+1.0d0))*((tk+fn)/(tk+2.0d0))*rxsq
217  trm(k) = t*b(k)
218  IF (abs(trm(k)).LT.tst) GO TO 80
219  s = s + trm(k)
220  tk = tk + 2.0d0
221  70 CONTINUE
222  80 CONTINUE
223  s = (s+t1)*tss
224  IF (xinc.EQ.0.0d0) GO TO 100
225 C-----------------------------------------------------------------------
226 C BACKWARD RECUR FROM XDMY TO X
227 C-----------------------------------------------------------------------
228  nx = int(xinc)
229  np = nn + 1
230  IF (nx.GT.nmax) GO TO 390
231  IF (nn.EQ.0) GO TO 160
232  xm = xinc - 1.0d0
233  fx = x + xm
234 C-----------------------------------------------------------------------
235 C THIS LOOP SHOULD NOT BE CHANGED. FX IS ACCURATE WHEN X IS SMALL
236 C-----------------------------------------------------------------------
237  DO 90 i=1,nx
238  trmr(i) = fx**(-np)
239  s = s + trmr(i)
240  xm = xm - 1.0d0
241  fx = x + xm
242  90 CONTINUE
243  100 CONTINUE
244  ans(mm) = s
245  IF (fn.EQ.0) GO TO 180
246 C-----------------------------------------------------------------------
247 C GENERATE LOWER DERIVATIVES, J.LT.N+MM-1
248 C-----------------------------------------------------------------------
249  IF (mm.EQ.1) RETURN
250  DO 150 j=2,mm
251  fn = fn - 1
252  tss = tss*xdmy
253  t1 = tt
254  IF (fn.NE.0) t1 = tt + 1.0d0/fn
255  t = (fn+1)*ta
256  s = t*b(3)
257  IF (abs(s).LT.tst) GO TO 120
258  tk = 4 + fn
259  DO 110 k=4,22
260  trm(k) = trm(k)*(fn+1)/tk
261  IF (abs(trm(k)).LT.tst) GO TO 120
262  s = s + trm(k)
263  tk = tk + 2.0d0
264  110 CONTINUE
265  120 CONTINUE
266  s = (s+t1)*tss
267  IF (xinc.EQ.0.0d0) GO TO 140
268  IF (fn.EQ.0) GO TO 160
269  xm = xinc - 1.0d0
270  fx = x + xm
271  DO 130 i=1,nx
272  trmr(i) = trmr(i)*fx
273  s = s + trmr(i)
274  xm = xm - 1.0d0
275  fx = x + xm
276  130 CONTINUE
277  140 CONTINUE
278  mx = mm - j + 1
279  ans(mx) = s
280  IF (fn.EQ.0) GO TO 180
281  150 CONTINUE
282  RETURN
283 C-----------------------------------------------------------------------
284 C RECURSION FOR N = 0
285 C-----------------------------------------------------------------------
286  160 CONTINUE
287  DO 170 i=1,nx
288  s = s + 1.0d0/(x+nx-i)
289  170 CONTINUE
290  180 CONTINUE
291  IF (kode.EQ.2) GO TO 190
292  ans(1) = s - xdmln
293  RETURN
294  190 CONTINUE
295  IF (xdmy.EQ.x) RETURN
296  xq = xdmy/x
297  ans(1) = s - log(xq)
298  RETURN
299 C-----------------------------------------------------------------------
300 C COMPUTE BY SERIES (X+K)**(-(N+1)) , K=0,1,2,...
301 C-----------------------------------------------------------------------
302  200 CONTINUE
303  nn = int(fln) + 1
304  np = n + 1
305  t1 = (n+1)*xln
306  t = exp(-t1)
307  s = t
308  den = x
309  DO 210 i=1,nn
310  den = den + 1.0d0
311  trm(i) = den**(-np)
312  s = s + trm(i)
313  210 CONTINUE
314  ans(1) = s
315  IF (n.NE.0) GO TO 220
316  IF (kode.EQ.2) ans(1) = s + xln
317  220 CONTINUE
318  IF (mm.EQ.1) RETURN
319 C-----------------------------------------------------------------------
320 C GENERATE HIGHER DERIVATIVES, J.GT.N
321 C-----------------------------------------------------------------------
322  tol = wdtol/5.0d0
323  DO 250 j=2,mm
324  t = t/x
325  s = t
326  tols = t*tol
327  den = x
328  DO 230 i=1,nn
329  den = den + 1.0d0
330  trm(i) = trm(i)/den
331  s = s + trm(i)
332  IF (trm(i).LT.tols) GO TO 240
333  230 CONTINUE
334  240 CONTINUE
335  ans(j) = s
336  250 CONTINUE
337  RETURN
338 C-----------------------------------------------------------------------
339 C SMALL X.LT.UNIT ROUND OFF
340 C-----------------------------------------------------------------------
341  260 CONTINUE
342  ans(1) = x**(-n-1)
343  IF (mm.EQ.1) GO TO 280
344  k = 1
345  DO 270 i=2,mm
346  ans(k+1) = ans(k)/x
347  k = k + 1
348  270 CONTINUE
349  280 CONTINUE
350  IF (n.NE.0) RETURN
351  IF (kode.EQ.2) ans(1) = ans(1) + xln
352  RETURN
353  290 CONTINUE
354  IF (t.GT.0.0d0) GO TO 380
355  nz=0
356  ierr=2
357  RETURN
358  380 CONTINUE
359  nz=nz+1
360  ans(mm)=0.0d0
361  mm=mm-1
362  IF (mm.EQ.0) RETURN
363  GO TO 41
364  390 CONTINUE
365  nz=0
366  ierr=3
367  RETURN
368  END
max
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:230
min
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:207
dpsifn
subroutine dpsifn(X, N, KODE, M, ANS, NZ, IERR)
Definition: dpsifn.f:3
abs
static T abs(T x)
Definition: pr-output.cc:1678