GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
psifn.f
Go to the documentation of this file.
1 *DECK PSIFN
2  SUBROUTINE psifn (X, N, KODE, M, ANS, NZ, IERR)
3 C***BEGIN PROLOGUE PSIFN
4 C***PURPOSE Compute derivatives of the Psi function.
5 C***LIBRARY SLATEC
6 C***CATEGORY C7C
7 C***TYPE SINGLE 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 PSIFN:
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 PSIFN 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, PSIFN returns
30 C the scaled derivatives as described. KODE=2 is operative only
31 C when K=0 and in that case PSIFN returns -PSI(X) + LN(X). That
32 C is, the logarithmic behavior for large X is removed when KODE=1
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 PSIFN(X,0,1,1,ANS) results in
38 C ANS = -PSI(X)
39 C
40 C Input
41 C X - Argument, X .gt. 0.0E0
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
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 (=R1MACH(4)) and 1.0E-18 since critical constants
71 C are given to only 18 digits.
72 C
73 C DPSIFN is the Double Precision version of PSIFN.
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 I1MACH, R1MACH
100 C***REVISION HISTORY (YYMMDD)
101 C 820601 DATE WRITTEN
102 C 890531 Changed all specific intrinsics to generic. (WRB)
103 C 890531 REVISION DATE from Version 3.2
104 C 891214 Prologue converted to Version 4.0 format. (BAB)
105 C 920501 Reformatted the REFERENCES section. (WRB)
106 C***END PROLOGUE PSIFN
107  INTEGER I, IERR, J, K, KODE, M, MM, MX, N, NMAX, NN, NP, NX, NZ
108  INTEGER I1MACH
109  REAL ANS, ARG, B, DEN, ELIM, EPS, FLN, FN, FNP, FNS, FX, RLN,
110  * RXSQ, R1M4, R1M5, S, SLOPE, T, TA, TK, TOL, TOLS, TRM, TRMR,
111  * TSS, TST, TT, T1, T2, WDTOL, X, XDMLN, XDMY, XINC, XLN, XM,
112  * XMIN, XQ, YINT
113  REAL R1MACH
114  dimension b(22), trm(22), trmr(100), ans(*)
115  SAVE nmax, b
116  DATA nmax /100/
117 C-----------------------------------------------------------------------
118 C BERNOULLI NUMBERS
119 C-----------------------------------------------------------------------
120  DATA b(1), b(2), b(3), b(4), b(5), b(6), b(7), b(8), b(9), b(10),
121  * b(11), b(12), b(13), b(14), b(15), b(16), b(17), b(18), b(19),
122  * b(20), b(21), b(22) /1.00000000000000000e+00,
123  * -5.00000000000000000e-01,1.66666666666666667e-01,
124  * -3.33333333333333333e-02,2.38095238095238095e-02,
125  * -3.33333333333333333e-02,7.57575757575757576e-02,
126  * -2.53113553113553114e-01,1.16666666666666667e+00,
127  * -7.09215686274509804e+00,5.49711779448621554e+01,
128  * -5.29124242424242424e+02,6.19212318840579710e+03,
129  * -8.65802531135531136e+04,1.42551716666666667e+06,
130  * -2.72982310678160920e+07,6.01580873900642368e+08,
131  * -1.51163157670921569e+10,4.29614643061166667e+11,
132  * -1.37116552050883328e+13,4.88332318973593167e+14,
133  * -1.92965793419400681e+16/
134 C
135 C***FIRST EXECUTABLE STATEMENT PSIFN
136  ierr = 0
137  nz=0
138  IF (x.LE.0.0e0) ierr=1
139  IF (n.LT.0) ierr=1
140  IF (kode.LT.1 .OR. kode.GT.2) ierr=1
141  IF (m.LT.1) ierr=1
142  IF (ierr.NE.0) RETURN
143  mm=m
144  nx = min(-i1mach(12),i1mach(13))
145  r1m5 = r1mach(5)
146  r1m4 = r1mach(4)*0.5e0
147  wdtol = max(r1m4,0.5e-18)
148 C-----------------------------------------------------------------------
149 C ELIM = APPROXIMATE EXPONENTIAL OVER AND UNDERFLOW LIMIT
150 C-----------------------------------------------------------------------
151  elim = 2.302e0*(nx*r1m5-3.0e0)
152  xln = log(x)
153  41 CONTINUE
154  nn = n + mm - 1
155  fn = nn
156  fnp = fn + 1.0e0
157  t = fnp*xln
158 C-----------------------------------------------------------------------
159 C OVERFLOW AND UNDERFLOW TEST FOR SMALL AND LARGE X
160 C-----------------------------------------------------------------------
161  IF (abs(t).GT.elim) GO TO 290
162  IF (x.LT.wdtol) GO TO 260
163 C-----------------------------------------------------------------------
164 C COMPUTE XMIN AND THE NUMBER OF TERMS OF THE SERIES, FLN+1
165 C-----------------------------------------------------------------------
166  rln = r1m5*i1mach(11)
167  rln = min(rln,18.06e0)
168  fln = max(rln,3.0e0) - 3.0e0
169  yint = 3.50e0 + 0.40e0*fln
170  slope = 0.21e0 + fln*(0.0006038e0*fln+0.008677e0)
171  xm = yint + slope*fn
172  mx = int(xm) + 1
173  xmin = mx
174  IF (n.EQ.0) GO TO 50
175  xm = -2.302e0*rln - min(0.0e0,xln)
176  fns = n
177  arg = xm/fns
178  arg = min(0.0e0,arg)
179  eps = exp(arg)
180  xm = 1.0e0 - eps
181  IF (abs(arg).LT.1.0e-3) xm = -arg
182  fln = x*xm/eps
183  xm = xmin - x
184  IF (xm.GT.7.0e0 .AND. fln.LT.15.0e0) GO TO 200
185  50 CONTINUE
186  xdmy = x
187  xdmln = xln
188  xinc = 0.0e0
189  IF (x.GE.xmin) GO TO 60
190  nx = int(x)
191  xinc = xmin - nx
192  xdmy = x + xinc
193  xdmln = log(xdmy)
194  60 CONTINUE
195 C-----------------------------------------------------------------------
196 C GENERATE W(N+MM-1,X) BY THE ASYMPTOTIC EXPANSION
197 C-----------------------------------------------------------------------
198  t = fn*xdmln
199  t1 = xdmln + xdmln
200  t2 = t + xdmln
201  tk = max(abs(t),abs(t1),abs(t2))
202  IF (tk.GT.elim) GO TO 380
203  tss = exp(-t)
204  tt = 0.5e0/xdmy
205  t1 = tt
206  tst = wdtol*tt
207  IF (nn.NE.0) t1 = tt + 1.0e0/fn
208  rxsq = 1.0e0/(xdmy*xdmy)
209  ta = 0.5e0*rxsq
210  t = fnp*ta
211  s = t*b(3)
212  IF (abs(s).LT.tst) GO TO 80
213  tk = 2.0e0
214  DO 70 k=4,22
215  t = t*((tk+fn+1.0e0)/(tk+1.0e0))*((tk+fn)/(tk+2.0e0))*rxsq
216  trm(k) = t*b(k)
217  IF (abs(trm(k)).LT.tst) GO TO 80
218  s = s + trm(k)
219  tk = tk + 2.0e0
220  70 CONTINUE
221  80 CONTINUE
222  s = (s+t1)*tss
223  IF (xinc.EQ.0.0e0) GO TO 100
224 C-----------------------------------------------------------------------
225 C BACKWARD RECUR FROM XDMY TO X
226 C-----------------------------------------------------------------------
227  nx = int(xinc)
228  np = nn + 1
229  IF (nx.GT.nmax) GO TO 390
230  IF (nn.EQ.0) GO TO 160
231  xm = xinc - 1.0e0
232  fx = x + xm
233 C-----------------------------------------------------------------------
234 C THIS LOOP SHOULD NOT BE CHANGED. FX IS ACCURATE WHEN X IS SMALL
235 C-----------------------------------------------------------------------
236  DO 90 i=1,nx
237  trmr(i) = fx**(-np)
238  s = s + trmr(i)
239  xm = xm - 1.0e0
240  fx = x + xm
241  90 CONTINUE
242  100 CONTINUE
243  ans(mm) = s
244  IF (fn.EQ.0.0e0) GO TO 180
245 C-----------------------------------------------------------------------
246 C GENERATE LOWER DERIVATIVES, J.LT.N+MM-1
247 C-----------------------------------------------------------------------
248  IF (mm.EQ.1) RETURN
249  DO 150 j=2,mm
250  fnp = fn
251  fn = fn - 1.0e0
252  tss = tss*xdmy
253  t1 = tt
254  IF (fn.NE.0.0e0) t1 = tt + 1.0e0/fn
255  t = fnp*ta
256  s = t*b(3)
257  IF (abs(s).LT.tst) GO TO 120
258  tk = 3.0e0 + fnp
259  DO 110 k=4,22
260  trm(k) = trm(k)*fnp/tk
261  IF (abs(trm(k)).LT.tst) GO TO 120
262  s = s + trm(k)
263  tk = tk + 2.0e0
264  110 CONTINUE
265  120 CONTINUE
266  s = (s+t1)*tss
267  IF (xinc.EQ.0.0e0) GO TO 140
268  IF (fn.EQ.0.0e0) GO TO 160
269  xm = xinc - 1.0e0
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.0e0
275  fx = x + xm
276  130 CONTINUE
277  140 CONTINUE
278  mx = mm - j + 1
279  ans(mx) = s
280  IF (fn.EQ.0.0e0) 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.0e0/(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 = (fns+1.0e0)*xln
306  t = exp(-t1)
307  s = t
308  den = x
309  DO 210 i=1,nn
310  den = den + 1.0e0
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.0e0
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.0e0
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.0e0) GO TO 380
355  nz=0
356  ierr=2
357  RETURN
358  380 CONTINUE
359  nz=nz+1
360  ans(mm)=0.0e0
361  mm=mm-1
362  IF(mm.EQ.0) RETURN
363  GO TO 41
364  390 CONTINUE
365  ierr=3
366  nz=0
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
abs
static T abs(T x)
Definition: pr-output.cc:1678
psifn
subroutine psifn(X, N, KODE, M, ANS, NZ, IERR)
Definition: psifn.f:3