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

# Syntax

C# |
---|

public static double s14aa( double x, out int ifail ) |

Visual Basic |
---|

Public Shared Function s14aa ( _ x As Double, _ <OutAttribute> ByRef ifail As Integer _ ) As Double |

Visual C++ |
---|

public: static double s14aa( double x, [OutAttribute] int% ifail ) |

F# |
---|

static member s14aa : x : float * ifail : int byref -> float |

#### Parameters

- x
- Type: System..::..Double
*On entry*: the argument $x$ of the function.*Constraint*: ${\mathbf{x}}$ must not be zero or a negative integer.

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

#### Return Value

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

# Description

s14aa evaluates an approximation to the gamma function $\Gamma \left(x\right)$. The method is based on the Chebyshev expansion:

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:

$$\Gamma \left(1+u\right)=\sum _{r=0}^{\prime}{a}_{r}{T}_{r}\left(t\right)\text{, \hspace{1em} where}0\le u<1,t=2u-1\text{,}$$ |

- for $N>0$, $\text{\hspace{1em}}\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{\hspace{1em}}\Gamma \left(x\right)=\Gamma \left(1+u\right)$,
- for $N<0$, $\text{\hspace{1em}}\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 method:

(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# Error Indicators and Warnings

Errors or warnings detected by the method:

- ${\mathbf{ifail}}=1$
- The argument is too large. On failure the method 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 failure the method returns zero.

- ${\mathbf{ifail}}=3$
- The argument is too close to zero. On failure the method returns the approximate value of $\Gamma \left(x\right)$ at the nearest valid argument.

- ${\mathbf{ifail}}=4$
- The argument is a negative integer, at which value $\Gamma \left(x\right)$ is infinite. On failure the method returns a large positive value.

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

(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|$.

$$\epsilon \simeq \left|x\Psi \left(x\right)\right|\delta $$ |

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

However since $\Gamma \left(x\right)$ increases rapidly with $x$, the method 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$.)

$$\epsilon \simeq \delta x\mathrm{ln}\u200ax\text{.}$$ |

# Parallelism and Performance

None.

# Further Comments

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.

Example program (C#): s14aae.cs