NAG Library Routine Document
s22bef (hyperg_gauss_real)
1
Purpose
s22bef returns a value for the Gauss hypergeometric function ${}_{2}F_{1}\left(a,b;c;x\right)$ for real parameters $a,b$ and $c$, and real argument $x$.
2
Specification
Fortran Interface
Real (Kind=nag_wp)  ::  s22bef  Integer, Intent (Inout)  ::  ifail  Real (Kind=nag_wp), Intent (In)  ::  a, b, c, x 

C Header Interface
#include <nagmk26.h>
double 
s22bef_ (const double *a, const double *b, const double *c, const double *x, Integer *ifail) 

3
Description
s22bef returns a value for the Gauss hypergeometric function ${}_{2}F_{1}\left(a,b;c;x\right)$ for real parameters $a$, $b$ and $c$, and for real argument $x$.
The associated routine
s22bff performs the same operations, but returns
${}_{2}F_{1}\left(a,b;c;x\right)$ in the scaled form
${}_{2}F_{1}\left(a,b;c;x\right)={f}_{\mathrm{fr}}\times {2}^{{f}_{\mathrm{sc}}}$ to allow calculations to be performed when
${}_{2}F_{1}\left(a,b;c;x\right)$ is not representable as a single working precision number. It also accepts the parameters
$a$,
$b$ and
$c$ as summations of an integer and a decimal fraction, giving higher accuracy when any are close to an integer.
The Gauss hypergeometric function is a solution to the hypergeometric differential equation,
For
$\leftx\right<1$, it may be defined by the Gauss series,
where
${\left(a\right)}_{s}=1\left(a\right)\left(a+1\right)\left(a+2\right)\dots \left(a+s1\right)$ is the rising factorial of
$a$.
${}_{2}F_{1}\left(a,b;c;x\right)$ is undefined for
$c=0$ or
$c$ a negative integer.
For $\leftx\right<1$, the series is absolutely convergent and ${}_{2}F_{1}\left(a,b;c;x\right)$ is finite.
For
$x<1$, linear transformations of the form,
exist, where
${x}_{1}$,
${x}_{2}\in \left(0,1\right]$.
${C}_{1}$ and
${C}_{2}$ are real valued functions of the parameters and argument, typically involving products of gamma functions. When these are degenerate, finite limiting cases exist. Hence for
$x<0$,
${}_{2}F_{1}\left(a,b;c;x\right)$ is defined by analytic continuation, and for
$x<1$,
${}_{2}F_{1}\left(a,b;c;x\right)$ is real and finite.
For
$x=1$, the following apply:
 If $c>a+b$, ${}_{2}F_{1}\left(a,b;c;1\right)=\frac{\Gamma \left(c\right)\Gamma \left(cab\right)}{\Gamma \left(ca\right)\Gamma \left(cb\right)}$, and hence is finite. Solutions also exist for the degenerate cases where $ca$ or $cb$ are negative integers or zero.
 If $c\le a+b$, ${}_{2}F_{1}\left(a,b;c;1\right)$ is infinite, and the sign of ${}_{2}F_{1}\left(a,b;c;1\right)$ is determinable as $x$ approaches $1$ from below.
In the complex plane, the principal branch of ${}_{2}F_{1}\left(a,b;c;z\right)$ is taken along the real axis from $x=1.0$ increasing. ${}_{2}F_{1}\left(a,b;c;z\right)$ is multivalued along this branch, and for real parameters $a,b$ and $c$ is typically not real valued. As such, this routine will not compute a solution when $x>1$.
The solution strategy used by this routine is primarily dependent upon the value of the argument $x$. Once trivial cases and the case $x=1.0$ are eliminated, this proceeds as follows.
For
$0<x\le 0.5$, sets of safe parameters
$\left\{{\alpha}_{i,j};{\beta}_{i,j};{\zeta}_{i,j};{\chi}_{j}\left1\le j\le 2\right;1\le i\le 4\right\}$ are determined, such that the values of
${}_{2}F_{1}\left({a}_{j},{b}_{j};{c}_{j};{x}_{j}\right)$ required for an appropriate transformation of the type
(3) may be calculated either directly or using recurrence relations from the solutions of
${}_{2}F_{1}\left({\alpha}_{i,j},{\beta}_{i,j};{\zeta}_{i,j};{\chi}_{j}\right)$. If
$c$ is positive, then only transformations with
${C}_{2}=0.0$ will be used, implying only
${}_{2}F_{1}\left({a}_{1},{b}_{1};{c}_{1};{x}_{1}\right)$ will be required, with the transformed argument
${x}_{1}=x$. If
$c$ is negative, in some cases a transformation with
${C}_{2}\ne 0.0$ will be used, with the argument
${x}_{2}=1.0x$. The routine then cycles through these sets until acceptable solutions are generated. If no computation produces an accurate answer, the least inaccurate answer is selected to complete the computation. See
Section 7.
For $0.5<x<1.0$, an identical approach is first used with the argument $x$. Should this fail, a linear transformation resulting in both transformed arguments satisfying ${x}_{j}=1.0x$ is employed, and the above strategy for $0<x\le 0.5$ is utilized on both components. Further transformations in these subcomputations are however, limited to single terms with no argument transformation.
For $x<0$, a linear transformation mapping the argument $x$ to the interval $\left(0,0.5\right]$ is first employed. The strategy for $0<x\le 0.5$ is then used on each component, including possible further two term transforms. To avoid some degenerate cases, a transform mapping the argument $x$ to $\left[0.5,1\right)$ may also be used.
In addition to the above restrictions on
$c$ and
$x$, an artificial bound,
arbnd, is placed on the magnitudes of
$a,b,c$ and
$x$ to minimize the occurrence of overflow in internal calculations, particularly those involving real to integer conversions.
$\mathit{arbnd}=0.0001\times {I}_{\mathrm{max}}$, where
${I}_{\mathrm{max}}$ is the largest machine integer (see
x02bbf). It should however, not be assumed that this routine will produce accurate answers for all values of
$a,b,c$ and
$x$ satisfying this criterion.
This routine also tests for nonfinite values of the parameters and argument on entry, and assigns nonfinite values upon completion if appropriate. See
Section 9 and
Chapter X07.
Please consult the
NIST Digital Library of Mathematical Functions for a detailed discussion of the Gauss hypergeometric function including special cases, transformations, relations and asymptotic approximations.
4
References
Pearson J (2009) Computation of hypergeometric functions MSc Dissertation, Mathematical Institute, University of Oxford
5
Arguments
 1: $\mathbf{a}$ – Real (Kind=nag_wp)Input

On entry: the parameter $a$.
Constraint:
$\left{\mathbf{a}}\right\le \mathit{arbnd}$.
 2: $\mathbf{b}$ – Real (Kind=nag_wp)Input

On entry: the parameter $b$.
Constraint:
$\left{\mathbf{b}}\right\le \mathit{arbnd}$.
 3: $\mathbf{c}$ – Real (Kind=nag_wp)Input

On entry: the parameter $c$.
Constraints:
 $\left{\mathbf{c}}\right\le \mathit{arbnd}$;
 ${\mathbf{c}}\ne 0,1,2,\dots $.
 4: $\mathbf{x}$ – Real (Kind=nag_wp)Input

On entry: the argument $x$.
Constraint:
$\mathit{arbnd}<{\mathbf{x}}\le 1$.
 5: $\mathbf{ifail}$ – IntegerInput/Output

On entry:
ifail must be set to
$0$,
$1\text{or}1$. If you are unfamiliar with this argument you should refer to
Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, if you are not familiar with this argument, the recommended value is
$0$.
When the value $\mathbf{1}\text{or}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}=1$

Underflow occurred during the evaluation of ${}_{2}F_{1}\left(a,b;c;x\right)$. The returned value may be inaccurate.
 ${\mathbf{ifail}}=2$

All approximations have completed, and the final residual estimate indicates some precision may have been lost.
$\text{Relative residual}=\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=3$

All approximations have completed, and the final residual estimate indicates no accuracy can be guaranteed.
$\text{Relative residual}=\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=4$

On entry, ${\mathbf{x}}=\u2329\mathit{\text{value}}\u232a$, $c=\u2329\mathit{\text{value}}\u232a$, $a+b=\u2329\mathit{\text{value}}\u232a$.
${}_{2}F_{1}\left(a,b;c;1\right)$ is infinite in the case $c\le a+b$.
 ${\mathbf{ifail}}=5$

On completion, overflow occurred in the evaluation of ${}_{2}F_{1}\left(a,b;c;x\right)$.
 ${\mathbf{ifail}}=6$

Overflow occurred in a subcalculation of ${}_{2}F_{1}\left(a,b;c;x\right)$. The result may or may not be infinite.
 ${\mathbf{ifail}}=9$

An internal calculation has resulted in an undefined result.
 ${\mathbf{ifail}}=11$

On entry,
a does not satisfy
$\left{\mathbf{a}}\right\le \mathit{arbnd}=\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=21$

On entry,
b does not satisfy
$\left{\mathbf{b}}\right\le \mathit{arbnd}=\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=31$

On entry,
c does not satisfy
$\left{\mathbf{c}}\right\le \mathit{arbnd}=\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=32$

On entry, ${\mathbf{c}}=\u2329\mathit{\text{value}}\u232a$.
${}_{2}F_{1}\left(a,b;c;x\right)$ is undefined when $c$ is zero or a negative integer.
 ${\mathbf{ifail}}=41$

On entry,
x does not satisfy
$\left{\mathbf{x}}\right\le \mathit{arbnd}=\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=42$

On entry, ${\mathbf{x}}=\u2329\mathit{\text{value}}\u232a$.
In general, ${}_{2}F_{1}\left(a,b;c;x\right)$ is not real valued when $x>1$.
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.9 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 3.8 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 3.7 in How to Use the NAG Library and its Documentation for further information.
7
Accuracy
In general, if
${\mathbf{ifail}}={\mathbf{0}}$, the value of
${}_{2}F_{1}\left(a,b;c;x\right)$ may be assumed accurate, with the possible loss of one or two decimal places. Assuming the result does not under or overflow, an error estimate
$\mathit{res}$ is made internally using equation
(1). If the magnitude of
$\mathit{res}$ is sufficiently large, a
nonzero
ifail
will be returned. Specifically,
${\mathbf{ifail}}={\mathbf{0}}$ or ${\mathbf{1}}$ 
$\mathit{res}\le 1000\epsilon $ 
${\mathbf{ifail}}={\mathbf{2}}$ 
$1000\epsilon <\mathit{res}\le 0.1$ 
${\mathbf{ifail}}={\mathbf{3}}$ 
$\mathit{res}>0.1$ 
where
$\epsilon $ is the
machine precision as returned by
x02ajf.
A further estimate of the residual can be constructed using equation
(1), and the differential identity,
This estimate is however, dependent upon the error involved in approximating ${}_{2}F_{1}\left(a+1,b+1;c+1;x\right)$ and ${}_{2}F_{1}\left(a+2,b+2;c+2;x\right)$.
Furthermore, the accuracy of the solution, and the error estimate, can be dependent upon the accuracy of the decimal fraction of the input parameters
$a$ and
$b$. For example, if
$c={c}_{i}+{c}_{r}=100+\text{1.0E\u22126}$, then on a machine with
$16$ decimal digits of precision, the internal calculation of
${c}_{r}$ will only be accurate to
$8$ decimal places. This can subsequently pollute the final solution by several decimal places without affecting the residual estimate as greatly. Should you require higher accuracy in such regions, then you should use
s22bff, which requires you to supply the correct decimal fraction.
8
Parallelism and Performance
s22bef is not threaded in any implementation.
s22bef returns nonfinite values when appropriate. See
Chapter X07 for more information on the definitions of nonfinite values.
Should a nonfinite value be returned, this will be indicated in the value of
ifail, as detailed in the following cases.
If
${\mathbf{ifail}}={\mathbf{0}}$, or
${\mathbf{ifail}}={\mathbf{1}}$,
${\mathbf{2}}$ or
${\mathbf{3}}$, a finite value will have been returned with an approximate accuracy as detailed in
Section 7.
If ${\mathbf{ifail}}={\mathbf{4}}$ then ${}_{2}F_{1}\left(a,b;c;x\right)$ is infinite, and a signed infinity will have been returned. The sign of the infinity should be correct when taking the limit as $x$ approaches $1$ from below.
If
${\mathbf{ifail}}={\mathbf{5}}$ then upon completion,
$\left{}_{2}F_{1}\left(a,b;c;x\right)\right>{R}_{\mathrm{max}}$, where
${R}_{\mathrm{max}}$ is the largest machine number given by
x02alf, and hence is too large to be representable. The result will be returned as a signed infinity. The sign should be correct.
If ${\mathbf{ifail}}={\mathbf{6}}$ then overflow occurred during a subcalculation of ${}_{2}F_{1}\left(a,b;c;x\right)$. A signed infinity will have been returned, however, there is no guarantee that this is representative of either the magnitude or the sign of ${}_{2}F_{1}\left(a,b;c;x\right)$.
For all other error exits,
s22bef will return a signalling NaN (see
x07bbf).
If ${\mathbf{ifail}}={\mathbf{9}}$ then an internal computation produced an undefined result. This may occur when two terms overflow with opposite signs, and the result is dependent upon their summation for example.
If
${\mathbf{ifail}}={\mathbf{32}}$ then
$c$ is too close to a negative integer or zero on entry, and
${}_{2}F_{1}\left(a,b;c;x\right)$ is considered undefined. Note, this will also be the case when
$c$ is a negative integer, and a (possibly trivial) linear transformation of the form
(3) would result in either:
(i) 
all ${c}_{j}$ not being negative integers, 
(ii) 
for any ${c}_{j}$ which remain as negative integers, one of the corresponding parameters ${a}_{j}$ or ${b}_{j}$ is a negative integer of magnitude less than ${c}_{j}$. 
In the first case, the transformation coefficients
${C}_{j}\left({a}_{j},{b}_{j},{c}_{j},{x}_{j}\right)$ are typically either infinite or undefined, preventing a solution being constructed. In the second case, the series
(2) will terminate before the degenerate term, resulting in a polynomial of fixed degree, and hence potentially a finite solution.
If ${\mathbf{ifail}}={\mathbf{11}}$, ${\mathbf{21}}$, ${\mathbf{31}}$ or ${\mathbf{41}}$ then no computation will have been performed. The actual solution may however, be finite.
${\mathbf{ifail}}={\mathbf{42}}$ indicates $x>1$. Hence the requested solution is on the boundary of the principal branch of ${}_{2}F_{1}\left(a,b;c;x\right)$, and hence is multivalued, typically with a nonzero imaginary component. It is however, strictly finite.
10
Example
This example evaluates ${}_{2}F_{1}\left(a,b;c;x\right)$ at a fixed set of parameters $a,b$ and $c$, and for several values for the argument $x$.
10.1
Program Text
Program Text (s22befe.f90)
10.2
Program Data
None.
10.3
Program Results
Program Results (s22befe.r)