Mercurial > octave-antonio

annotate liboctave/cruft/slatec-fn/dpsifn.f @ 20161:65e22ba879f0

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