Mathlib : A C Library of Special Functions
Copyright (C) 1998 Ross Ihaka
Copyright (C) 2000-2007 the R Development Core Team
Copyright (C) 2004 The R Foundation
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, a copy is available at
http://www.r-project.org/Licenses/
SYNOPSIS
#include
double[] dpsifn(double x, int n, int kode, int m)
double digamma(double x);
double trigamma(double x)
double tetragamma(double x)
double pentagamma(double x)
double psigamma(double x, double n)
DESCRIPTION
Compute the derivatives of the psi function
and polygamma functions.
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 dpsifn(x, 0, 1, 1, ans) results in ans = -psi(x).
INPUT
x - argument, x > 0.
n - first member of the sequence, 0 <= n <= 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 >= 1
OUTPUT
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 == 0, a normal return
nz != 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 1e-18 since critical constants
are given to only 18 digits.
The basic method of evaluation is the asymptotic expansion
for large x >= 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.
AUTHOR
Amos, D. E. (Fortran)
Ross Ihaka (C Translation)
Martin Maechler (x < 0, and psigamma())
Roby Joehanes (Java translation)
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, (1983). "A Portable Fortran Subroutine for
Derivatives of the Psi Function", Algorithm 610,
TOMS 9(4), pp. 494-502.