s14ad returns a sequence of values of scaled derivatives of the psi function (also known as the digamma function).
- Type: System..::..DoubleOn entry: the argument of the function.Constraint: .
- Type: System..::..Int32On entry: the index of the first member of the sequence of functions.Constraint: .
- Type: System..::..Int32On entry: the number of members required in the sequence , for .Constraint: .
- Type: array<System..::..Double>()An array of size [m]On exit: the first elements of ans contain the required values , for .
s14ad computes values of the function
for , , , where is the psi function
and denotes the th derivative of .
The method is derived from the method PSIFN in Amos (1983). The basic method of evaluation of is the asymptotic series
for large greater than a machine-dependent value , followed by backward recurrence using
for smaller values of , where when , when , and , , are the Bernoulli numbers.
When is large, the above procedure may be inefficient, and the expansion
which converges rapidly for large , is used instead.
Abramowitz M and Stegun I A (1972) Handbook of Mathematical Functions (3rd Edition) Dover Publications
Amos D E (1983) Algorithm 610: A portable FORTRAN subroutine for derivatives of the psi function ACM Trans. Math. Software 9 494–502
Errors or warnings detected by the method:
On entry, . On entry, . On entry, .
All constants in s14ad are given to approximately digits of precision. Calling the number of digits of precision in the floating-point arithmetic being used , then clearly the maximum number of correct digits in the results obtained is limited by . Empirical tests of s14ad, taking values of in the range , and in the range , have shown that the maximum relative error is a loss of approximately two decimal places of precision. Tests with , i.e., testing the function , have shown somewhat better accuracy, except at points close to the zero of , , where only absolute accuracy can be obtained.
This example reads values of the argument from a file, evaluates the function at each value of and prints the results.