Mercurial > octave-antonio
annotate liboctave/cruft/slatec-fn/psifn.f @ 20161:65e22ba879f0
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 |
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 PSIFN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
2 SUBROUTINE PSIFN (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 PSIFN |
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 SINGLE 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 PSIFN: |
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 PSIFN 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, PSIFN 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 PSIFN 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=1 |
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 PSIFN(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 |
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.0E0 |
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 |
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 (=R1MACH(4)) and 1.0E-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 DPSIFN is the Double Precision version of PSIFN. |
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 I1MACH, R1MACH |
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 890531 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
|
104 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
|
105 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
|
106 C***END PROLOGUE PSIFN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
107 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
|
108 INTEGER I1MACH |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
109 REAL ANS, ARG, B, DEN, ELIM, EPS, FLN, FN, FNP, FNS, FX, RLN, |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
110 * RXSQ, R1M4, R1M5, S, SLOPE, T, TA, TK, TOL, TOLS, TRM, TRMR, |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
111 * TSS, TST, TT, T1, T2, WDTOL, X, XDMLN, XDMY, XINC, XLN, XM, |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
112 * XMIN, XQ, YINT |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
113 REAL R1MACH |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
114 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
|
115 SAVE NMAX, B |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
116 DATA NMAX /100/ |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
117 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
118 C BERNOULLI NUMBERS |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
119 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
120 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
|
121 * 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
|
122 * B(20), B(21), B(22) /1.00000000000000000E+00, |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
123 * -5.00000000000000000E-01,1.66666666666666667E-01, |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
124 * -3.33333333333333333E-02,2.38095238095238095E-02, |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
125 * -3.33333333333333333E-02,7.57575757575757576E-02, |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
126 * -2.53113553113553114E-01,1.16666666666666667E+00, |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
127 * -7.09215686274509804E+00,5.49711779448621554E+01, |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
128 * -5.29124242424242424E+02,6.19212318840579710E+03, |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
129 * -8.65802531135531136E+04,1.42551716666666667E+06, |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
130 * -2.72982310678160920E+07,6.01580873900642368E+08, |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
131 * -1.51163157670921569E+10,4.29614643061166667E+11, |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
132 * -1.37116552050883328E+13,4.88332318973593167E+14, |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
133 * -1.92965793419400681E+16/ |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
134 C |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
135 C***FIRST EXECUTABLE STATEMENT PSIFN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
136 IERR = 0 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
137 NZ=0 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
138 IF (X.LE.0.0E0) IERR=1 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
139 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
|
140 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
|
141 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
|
142 IF (IERR.NE.0) RETURN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
143 MM=M |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
144 NX = MIN(-I1MACH(12),I1MACH(13)) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
145 R1M5 = R1MACH(5) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
146 R1M4 = R1MACH(4)*0.5E0 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
147 WDTOL = MAX(R1M4,0.5E-18) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
148 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
149 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
|
150 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
151 ELIM = 2.302E0*(NX*R1M5-3.0E0) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
152 XLN = LOG(X) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
153 41 CONTINUE |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
154 NN = N + MM - 1 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
155 FN = NN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
156 FNP = FN + 1.0E0 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
157 T = FNP*XLN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
158 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
159 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
|
160 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
161 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
|
162 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
|
163 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
164 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
|
165 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
166 RLN = R1M5*I1MACH(11) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
167 RLN = MIN(RLN,18.06E0) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
168 FLN = MAX(RLN,3.0E0) - 3.0E0 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
169 YINT = 3.50E0 + 0.40E0*FLN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
170 SLOPE = 0.21E0 + FLN*(0.0006038E0*FLN+0.008677E0) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
171 XM = YINT + SLOPE*FN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
172 MX = INT(XM) + 1 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
173 XMIN = MX |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
174 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
|
175 XM = -2.302E0*RLN - MIN(0.0E0,XLN) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
176 FNS = N |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
177 ARG = XM/FNS |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
178 ARG = MIN(0.0E0,ARG) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
179 EPS = EXP(ARG) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
180 XM = 1.0E0 - EPS |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
181 IF (ABS(ARG).LT.1.0E-3) XM = -ARG |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
182 FLN = X*XM/EPS |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
183 XM = XMIN - X |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
184 IF (XM.GT.7.0E0 .AND. FLN.LT.15.0E0) GO TO 200 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
185 50 CONTINUE |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
186 XDMY = X |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
187 XDMLN = XLN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
188 XINC = 0.0E0 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
189 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
|
190 NX = INT(X) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
191 XINC = XMIN - NX |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
192 XDMY = X + XINC |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
193 XDMLN = LOG(XDMY) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
194 60 CONTINUE |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
195 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
196 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
|
197 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
198 T = FN*XDMLN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
199 T1 = XDMLN + XDMLN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
200 T2 = T + XDMLN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
201 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
|
202 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
|
203 TSS = EXP(-T) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
204 TT = 0.5E0/XDMY |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
205 T1 = TT |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
206 TST = WDTOL*TT |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
207 IF (NN.NE.0) T1 = TT + 1.0E0/FN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
208 RXSQ = 1.0E0/(XDMY*XDMY) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
209 TA = 0.5E0*RXSQ |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
210 T = FNP*TA |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
211 S = T*B(3) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
212 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
|
213 TK = 2.0E0 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
214 DO 70 K=4,22 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
215 T = T*((TK+FN+1.0E0)/(TK+1.0E0))*((TK+FN)/(TK+2.0E0))*RXSQ |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
216 TRM(K) = T*B(K) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
217 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
|
218 S = S + TRM(K) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
219 TK = TK + 2.0E0 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
220 70 CONTINUE |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
221 80 CONTINUE |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
222 S = (S+T1)*TSS |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
223 IF (XINC.EQ.0.0E0) GO TO 100 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
224 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
225 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
|
226 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
227 NX = INT(XINC) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
228 NP = NN + 1 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
229 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
|
230 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
|
231 XM = XINC - 1.0E0 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
232 FX = X + XM |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
233 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
234 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
|
235 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
236 DO 90 I=1,NX |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
237 TRMR(I) = FX**(-NP) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
238 S = S + TRMR(I) |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
239 XM = XM - 1.0E0 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
240 FX = X + XM |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
241 90 CONTINUE |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
242 100 CONTINUE |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
243 ANS(MM) = S |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
244 IF (FN.EQ.0.0E0) GO TO 180 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
245 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
246 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
|
247 C----------------------------------------------------------------------- |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
248 IF (MM.EQ.1) RETURN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
249 DO 150 J=2,MM |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
250 FNP = FN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
251 FN = FN - 1.0E0 |
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.0E0) T1 = TT + 1.0E0/FN |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
255 T = FNP*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 = 3.0E0 + FNP |
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)*FNP/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.0E0 |
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.0E0) 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.0E0) 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.0E0 |
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.0E0 |
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.0E0) 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.0E0/(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 = (FNS+1.0E0)*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.0E0 |
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.0E0 |
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.0E0 |
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.0E0) 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.0E0 |
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 IERR=3 |
65e22ba879f0
psi: add support to compute the polygamma function (kth-derivative).
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
366 NZ=0 |
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 |