Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_specfun_polygamma (s14ac)

## Purpose

nag_specfun_polygamma (s14ac) returns a value of the function $\psi \left(x\right)-\mathrm{ln}x$, where $\psi$ is the psi function $\psi \left(x\right)=\frac{d}{dx}\mathrm{ln}\Gamma \left(x\right)=\frac{{\Gamma }^{\prime }\left(x\right)}{\Gamma \left(x\right)}$.

## Syntax

[result, ifail] = s14ac(x)
[result, ifail] = nag_specfun_polygamma(x)

## Description

nag_specfun_polygamma (s14ac) returns a value of the function $\psi \left(x\right)-\mathrm{ln}x$. The psi function is computed without the logarithmic term so that when $x$ is large, sums or differences of psi functions may be computed without unnecessary loss of precision, by analytically combining the logarithmic terms. For example, the difference $d=\psi \left(x+\frac{1}{2}\right)-\psi \left(x\right)$ has an asymptotic behaviour for large $x$ given by $d\sim \mathrm{ln}\left(x+\frac{1}{2}\right)-\mathrm{ln}x+\mathit{O}\left(\frac{1}{{x}^{2}}\right)\sim \mathrm{ln}\left(1+\frac{1}{2x}\right)\sim \frac{1}{2x}$.
Computing $d$ directly would amount to subtracting two large numbers which are close to $\mathrm{ln}\left(x+\frac{1}{2}\right)$ and $\mathrm{ln}x$ to produce a small number close to $\frac{1}{2x}$, resulting in a loss of significant digits. However, using this function to compute $f\left(x\right)=\psi \left(x\right)-\mathrm{ln}x$, we can compute $d=f\left(x+\frac{1}{2}\right)-f\left(x\right)+\mathrm{ln}\left(1+\frac{1}{2x}\right)$, and the dominant logarithmic term may be computed accurately from its power series when $x$ is large. Thus we avoid the unnecessary loss of precision.
The function is derived from the function PSIFN in Amos (1983).

## 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

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{x}$ – double scalar
The argument $x$ of the function.
Constraint: ${\mathbf{x}}>0.0$.

None.

### Output Parameters

1:     $\mathrm{result}$ – double scalar
The result of the function.
2:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

Errors or warnings detected by the function:
${\mathbf{ifail}}=1$
 On entry, ${\mathbf{x}}\le 0.0$. nag_specfun_polygamma (s14ac) returns the value zero.
${\mathbf{ifail}}=2$
No result is computed because underflow is likely. The value of x is too large. nag_specfun_polygamma (s14ac) returns the value zero.
${\mathbf{ifail}}=3$
No result is computed because overflow is likely. The value of x is too small. nag_specfun_polygamma (s14ac) returns the value zero.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

All constants in nag_specfun_polygamma (s14ac) 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=\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(t,18\right)$.
With the above proviso, results returned by this function should be accurate almost to full precision, except at points close to the zero of $\psi \left(x\right)$, $x\simeq 1.461632$, where only absolute rather than relative accuracy can be obtained.

None.

## Example

The example program reads values of the argument $x$ from a file, evaluates the function at each value of $x$ and prints the results.
```function s14ac_example

fprintf('s14ac example results\n\n');

x = [0.1  0.5  3.6  8];
n = size(x,2);
result = x;

for j=1:n
[result(j), ifail] = s14ac(x(j));
end

disp('        x      Psi(x)-ln(x)');
fprintf('%12.4f%12.4f\n',[x; result]);

s14ac_plot;

function s14ac_plot
x = [0.1:0.1:8];

for j=1:numel(x)
[pml(j), ifail] = s14ac(x(j));
end

fig1 = figure;
plot(x,pml,'-r');
xlabel('x');
ylabel('\Psi(x) - ln(x)');
title('\Psi(x) - ln(x)');
axis([0 8 -7 0]);

```
```s14ac example results

x      Psi(x)-ln(x)
0.1000     -8.1212
0.5000     -1.2704
3.6000     -0.1453
8.0000     -0.0638
```