s14aaf evaluates an approximation to the gamma function
$\Gamma \left(x\right)$. The routine is based on the Chebyshev expansion:
where
$0\le u<1,t=2u1\text{,}$ and uses the property
$\Gamma \left(1+x\right)=x\Gamma \left(x\right)$. If
$x=N+1+u$ where
$N$ is integral and
$0\le u<1$ then it follows that:
for $N>0$, 
$\Gamma \left(x\right)=\left(x1\right)\left(x2\right)\cdots \left(xN\right)\Gamma \left(1+u\right)$, 
for $N=0$, 
$\Gamma \left(x\right)=\Gamma \left(1+u\right)$, 
for $N<0$, 
$\Gamma \left(x\right)=\frac{\Gamma \left(1+u\right)}{x\left(x+1\right)\left(x+2\right)\cdots \left(xN1\right)}$. 
There are four possible failures for this routine:

(i)if $x$ is too large, there is a danger of overflow since $\Gamma \left(x\right)$ could become too large to be represented in the machine;

(ii)if $x$ is too large and negative, there is a danger of underflow;

(iii)if $x$ is equal to a negative integer, $\Gamma \left(x\right)$ would overflow since it has poles at such points;

(iv)if $x$ is too near zero, there is again the danger of overflow on some machines. For small $x$, $\Gamma \left(x\right)\simeq 1/x$, and on some machines there exists a range of nonzero but small values of $x$ for which $1/x$ is larger than the greatest representable value.
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Let
$\delta $ and
$\epsilon $ be the relative errors in the argument and the result respectively. If
$\delta $ is somewhat larger than the
machine precision (i.e., is due to data errors etc.), then
$\epsilon $ and
$\delta $ are approximately related by:
(provided
$\epsilon $ is also greater than the representation error). Here
$\Psi \left(x\right)$ is the digamma function
$\frac{{\Gamma}^{\prime}\left(x\right)}{\Gamma \left(x\right)}$.
Figure 1 shows the behaviour of the error amplification factor
$\leftx\Psi \left(x\right)\right$.
There is clearly a severe, but unavoidable, loss of accuracy for arguments close to the poles of $\Gamma \left(x\right)$ at negative integers. However, relative accuracy is preserved near the pole at $x=0$ right up to the point of failure arising from the danger of overflow.
Also, accuracy will necessarily be lost as
$x$ becomes large since in this region
However, since
$\Gamma \left(x\right)$ increases rapidly with
$x$, the routine must fail due to the danger of overflow before this loss of accuracy is too great. (For example, for
$x=20$, the amplification factor
$\text{}\simeq 60$.)
None.