s14ad returns a sequence of values of scaled derivatives of the psi function

ψ (x)

(also known as the digamma function).

Syntax

C#
public static void s14ad( double x, int n, int m, double[] ans, out int ifail )

Visual Basic
Public Shared Sub s14ad ( _ x As Double, _ n As Integer, _ m As Integer, _ ans As Double(), _ <OutAttribute> ByRef ifail As Integer _ )

Visual C++
public: static void s14ad( double x, int n, int m, array<double>^ ans, [OutAttribute] int% ifail )

F#
static member s14ad : x : float * n : int * m : int * ans : float[] * ifail : int byref -> unit

Parameters

x: Type: System..::..Double
On entry: the argument $x$ of the function.

Constraint: $x > 0.0$ .

n: Type: System..::..Int32
On entry: the index of the first member $n$ of the sequence of functions.

Constraint: $n \geq 0$ .

m: Type: System..::..Int32
On entry: the number of members $m$ required in the sequence $w (k, x)$ , for $k = n, \dots, n + m - 1$ .

Constraint: $m \geq 1$ .

ans: Type: array<System..::..Double>[]()[][]
An array of size [m]
On exit: the first $m$ elements of ans contain the required values $w (k, x)$ , for $k = n, \dots, n + m - 1$ .

ifail: Type: System..::..Int32%
On exit: $ifail = 0$ unless the method detects an error or a warning has been flagged (see [Error Indicators and Warnings]).

Description

s14ad computes

m

values of the function

w (k, x) = \frac{{(- 1)}^{k + 1} ψ^{(k)} (x)}{k!},

for

x > 0

k = n

n + 1, \dots, n + m - 1

, where

ψ

is the psi function

ψ (x) = \frac{d}{d x} \ln Γ (x) = \frac{Γ^{'} (x)}{Γ (x)},

and

ψ^{(k)}

denotes the

k

th derivative of

ψ

The method is derived from the method PSIFN in Amos (1983). The basic method of evaluation of

w (k, x)

is the asymptotic series

w (k, x) \sim ε (k, x) + \frac{1}{2 x^{k + 1}} + \frac{1}{x^{k}} \sum_{j = 1}^{\infty} B_{2 j} \frac{(2 j + k - 1)!}{(2 j)! k! x^{2 j}}

for large

x

greater than a machine-dependent value

x_{\min}

, followed by backward recurrence using

w (k, x) = w (k, x + 1) + x^{- k - 1}

for smaller values of

x

, where

ε (k, x) = - \ln x

when

k = 0

ε (k, x) = \frac{1}{k x^{k}}

when

k > 0

, and

B_{2 j}

j = 1, 2, \dots

, are the Bernoulli numbers.

When

k

is large, the above procedure may be inefficient, and the expansion

w (k, x) = \sum_{j = 1}^{\infty} \frac{1}{{(x + j)}^{k + 1}},

which converges rapidly for large

k

, is used instead.

References

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

Error Indicators and Warnings

Errors or warnings detected by the method:

$ifail = 1$

On entry,

x \leq 0.0

$ifail = 2$

On entry,

n < 0

$ifail = 3$

On entry,

m < 1

$ifail = 4$: No results are returned because underflow is likely. Either x or $n + m - 1$ is too large. If possible, reduce the value of m and call s14ad again.

$ifail = 5$: No results are returned because overflow is likely. Either x is too small, or $n + m - 1$ is too large. If possible, reduce the value of m and call s14ad again.

$ifail = 6$: No results are returned because there is not enough internal workspace to continue computation. $n + m - 1$ may be too large. If possible, reduce the value of m and call s14ad again.

$ifail = -9000$: An error occured, see message report.
$ifail = -8000$: Negative dimension for array $〈value〉$
$ifail = -6000$: Invalid Parameters $〈value〉$

Accuracy

All constants in s14ad are given to approximately

18

digits of precision. Calling the number of digits of precision in the floating-point arithmetic being used

t

, then clearly the maximum number of correct digits in the results obtained is limited by

p = \min (t, 18)

. Empirical tests of s14ad, taking values of

x

in the range

0.0 < x < 50.0

, and

n

in the range

1 \leq n \leq 50

, have shown that the maximum relative error is a loss of approximately two decimal places of precision. Tests with

n = 0

, i.e., testing the function

- ψ (x)

, have shown somewhat better accuracy, except at points close to the zero of

ψ (x)

x ≃ 1.461632

, where only absolute accuracy can be obtained.

Parallelism and Performance

None.

Further Comments

The time taken for a call of s14ad is approximately proportional to

m

, plus a constant. In general, it is much cheaper to call s14ad with

m

greater than

1

to evaluate the function

w (k, x)

, for

k = n, \dots, n + m - 1

, rather than to make

m

separate calls of s14ad.

Example

This example reads values of the argument

x

from a file, evaluates the function at each value of

x

and prints the results.

Example program (C#): s14ade.cs

Example program data: s14ade.d

Example program results: s14ade.r

Syntax

Parameters

Description

References

Error Indicators and Warnings

Accuracy

Parallelism and Performance

Further Comments

Example

See Also