public class PolyGammaextends Object
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.