1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
|
=pod
# void dpsifn(double x, int n, int kode, int m, double *ans, int *nz, int *ierr)
/*
* Purpose: Compute Derivatives of the Psi Function.
*
* Author: Amos, D. E., (SNLA)
* C Translation by Ross Ihaka
*
* The following definitions are used in dpsifn:
*
* Definition 1
* psi(x) = d/dx (ln(gamma(x)), the first derivative of
* the log gamma function.
* Definition 2
* k k
* psi(k,x) = d /dx (psi(x)), the k-th derivative of psi(x).
*
*
* Dpsifn computes a sequence of scaled derivatives of
* the psi function; i.e. for fixed x and m it computes
* the m-member sequence
*
* ((-1)**(k+1)/gamma(k+1))*psi(k,x)
* for k = n,...,n+m-1
*
* where psi(k,x) is as defined above. For kode=1, dpsifn returns
* the scaled derivatives as described. kode=2 is operative only
* when k=0 and in that case dpsifn returns -psi(x) + ln(x). that
* is, the logarithmic behavior for large x is removed when kode=2
* and k=0. when sums or differences of psi functions are computed
* the logarithmic terms can be combined analytically and computed
* separately to help retain significant digits.
*
* Note that call dpsifn(x,0,1,1,ans) results in ans = -psi(x)
*
* Input x is double precision
* x - argument, x .gt. 0.0d0
* n - first member of the sequence, 0 .le. n .le. 100
* n=0 gives ans(1) = -psi(x) for kode=1
* -psi(x)+ln(x) for kode=2
* kode - selection parameter
* kode=1 returns scaled derivatives of the psi
* function.
* kode=2 returns scaled derivatives of the psi
* function except when n=0. in this case,
* ans(1) = -psi(x) + ln(x) is returned.
* m - number of members of the sequence, m.ge.1
*
* Output ans is double precision
* ans - a vector of length at least m whose first m
* components contain the sequence of derivatives
* scaled according to kode.
* nz - underflow flag
* nz.eq.0, a normal return
* nz.ne.0, underflow, last nz components of ans are
* set to zero, ans(m-k+1)=0.0, k=1,...,nz
* ierr - error flag
* ierr=0, a normal return, computation completed
* ierr=1, input error, no computation
* ierr=2, overflow, x too small or n+m-1 too
* large or both
* ierr=3, error, n too large. dimensioned
* array trmr(nmax) is not large enough for n
*
* The nominal computational accuracy is the maximum of unit
* roundoff (=d1mach(4)) and 1.0d-18 since critical constants
* are given to only 18 digits.
*
* Long Description:
*
* The basic method of evaluation is the asymptotic expansion
* for large x.ge.xmin followed by backward recursion on a two
* term recursion relation
*
* w(x+1) + x**(-n-1) = w(x).
*
* this is supplemented by a series
*
* sum( (x+k)**(-n-1) , k=0,1,2,... )
*
* which converges rapidly for large n. both xmin and the
* number of terms of the series are calculated from the unit
* roundoff of the machine environment.
*
* References:
*
* Handbook of Mathematical Functions,
* National Bureau of Standards Applied Mathematics Series 55,
* Edited by M. Abramowitz and I. A. Stegun, equations 6.3.5,
* 6.3.18, 6.4.6, 6.4.9 and 6.4.10, pp.258-260, 1964.
*
* D. E. Amos, A Portable Fortran Subroutine for Derivatives
* of the Psi Function, Algorithm 610,
* ACM Transactions on Mathematical Software 9, 4 (1983), pp. 494-502.
=cut
# void dpsifn(double x, int n, int kode, int m, double *ans, int *nz, int *ierr)
pp_addpm(<<'EOD');
# This interface is OK , but it should be improved
sub dpsifn {
my($x,$n,$kode,$m) = @_;
my($ans) = zeroes(float(),$m);
my ($b,$c) = (0,0);
dpsifn_in($x,$n,$kode,$m,$ans,$b,$c);
return ($ans,$b,$c);
}
EOD
pp_def("dpsifn_in",
Pars => 'x(); n(); kode(); m(); [o]ans(n); int [o]nz(); int [o]ierr()',
GenericTypes => [D],
Code => '
extern void dpsifn(double x, int n, int kode, int m, double *ans, int *nz, int *ierr);
int ierro, nzo;
dpsifn((double)$x(),(int)$n(),(int)$kode(),(int)$m(),$P(ans),&nzo,&ierro);
$nz() = nzo;
$ierr() = ierro;
');
|
