GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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
dpsifn
subroutine dpsifn(X, N, KODE, M, ANS, NZ, IERR)
Definition: dpsifn.f:2
log
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the base of natural logarithms The constant ex $e satisfies the equation log(e)
base_det::exp
int exp(void) const
Definition: DET.h:66
to
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
Definition: Quad-opts.cc:233
dimension
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
max
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:228
abs
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
min
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:205
Generated on Wed May 10 2017 15:42:57 for GNU Octave by   doxygen 1.8.8