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_gamma (s14aa)

## Purpose

nag_specfun_gamma (s14aa) returns the value of the gamma function $\Gamma \left(x\right)$, via the function name.

## Syntax

[result, ifail] = s14aa(x)
[result, ifail] = nag_specfun_gamma(x)

## Description

nag_specfun_gamma (s14aa) evaluates an approximation to the gamma function $\Gamma \left(x\right)$. The function is based on the Chebyshev expansion:
 $Γ1+u=∑r=0′arTrt, where ​ 0≤u<1,t=2u-1,$
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$, $\text{ }\Gamma \left(x\right)=\left(x-1\right)\left(x-2\right)\cdots \left(x-N\right)\Gamma \left(1+u\right)$,
• for $N=0$, $\text{ }\Gamma \left(x\right)=\Gamma \left(1+u\right)$,
• for $N<0$, $\text{ }\Gamma \left(x\right)=\frac{\Gamma \left(1+u\right)}{x\left(x+1\right)\left(x+2\right)\cdots \left(x-N-1\right)}$.
There are four possible failures for this function:
 (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.

## References

Abramowitz M and Stegun I A (1972) Handbook of Mathematical Functions (3rd Edition) Dover Publications

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{x}$ – double scalar
The argument $x$ of the function.
Constraint: ${\mathbf{x}}$ must not be zero or a negative integer.

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:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

${\mathbf{ifail}}=1$
The argument is too large. On soft failure the function returns the approximate value of $\Gamma \left(x\right)$ at the nearest valid argument.
${\mathbf{ifail}}=2$
The argument is too large and negative. On soft failure the function returns zero.
W  ${\mathbf{ifail}}=3$
The argument is too close to zero. On soft failure the function returns the approximate value of $\Gamma \left(x\right)$ at the nearest valid argument.
W  ${\mathbf{ifail}}=4$
The argument is a negative integer, at which value $\Gamma \left(x\right)$ is infinite. On soft failure the function returns a large positive value.
${\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

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:
 $ε≃xΨxδ$
(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 $\left|x\Psi \left(x\right)\right|$.
If $\delta$ is of the same order as machine precision, then rounding errors could make $\epsilon$ slightly larger than the above relation predicts.
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
 $ε≃δxln⁡x.$
However since $\Gamma \left(x\right)$ increases rapidly with $x$, the function 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$.) Figure 1

None.

## Example

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

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

x = [ 1    1.25     1.5     1.75     2       5       10     -1.5];
n = size(x,2);
result = x;

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

disp('      x        Gamma(x)');
fprintf('%12.3e%12.3e\n',[x; result]);

s14aa_plot;

function s14aa_plot
x = {[-3.99:0.01:-3.01]; [-2.99:0.05:-2.04]; [-1.9:0.1:-1.1];
[-0.9:0.1:-0.1]; [0.1:0.2:3.9]};

for k = 1:5
for j=1:numel(x{k})
[g{k}(j), ifail] = s14aa(x{k}(j));
end
end

fig1 = figure;
hold on;
for k = 1:5
plot(x{k},g{k},'-r');
end
xlabel('x');
ylabel('\Gamma(x)');
title('Gamma Function \Gamma(x)');
axis([-4 4 -5 5]);
hold off;
```
```s14aa example results

x        Gamma(x)
1.000e+00   1.000e+00
1.250e+00   9.064e-01
1.500e+00   8.862e-01
1.750e+00   9.191e-01
2.000e+00   1.000e+00
5.000e+00   2.400e+01
1.000e+01   3.629e+05
-1.500e+00   2.363e+00
``` 