NAG Library Routine Document

s22bff (hyperg_gauss_real_scaled)

1
Purpose

s22bff returns a value for the Gauss hypergeometric function F 1 2 a,b;c;x  for real parameters a,b and c, and real argument x. The result is returned in the scaled form F 1 2 a,b;c;x = ffr × 2 fsc .

2
Specification

Fortran Interface
Subroutine s22bff ( ani, adr, bni, bdr, cni, cdr, x, frf, scf, ifail)
Integer, Intent (Inout):: ifail
Integer, Intent (Out):: scf
Real (Kind=nag_wp), Intent (In):: ani, adr, bni, bdr, cni, cdr, x
Real (Kind=nag_wp), Intent (Out):: frf
C Header Interface
#include <nagmk26.h>
void  s22bff_ (const double *ani, const double *adr, const double *bni, const double *bdr, const double *cni, const double *cdr, const double *x, double *frf, Integer *scf, Integer *ifail)

3
Description

s22bff returns a value for the Gauss hypergeometric function F1 2 a,b;c;x  for real parameters a, b and c, and for real argument x.
The Gauss hypergeometric function is a solution to the hypergeometric differential equation,
x1-x d2 f dx2 + c- a+b+1 x d f dx - ab f = 0 . (1)
For x<1, it may be defined by the Gauss series,
F1 2 a,b;c;x = s=0 as bs cs s! xs = 1+ ab c x + aa+1 bb+1 cc+1 2! x2 + , (2)
where as = 1 a a+1 a+2 a+s-1  is the rising factorial of a . F1 2 a,b;c;x  is undefined for c=0 or c a negative integer.
For x<1, the series is absolutely convergent and F1 2 a,b;c;x  is finite.
For x<1, linear transformations of the form,
F1 2 a,b;c;x = C1 a1,b1,c1,x1 F1 2 a1, b1 ;c1;x1 + C2 a2,b2,c2,x2 F1 2 a2, b2 ;c2;x2 (3)
exist, where x1, x20,1. C1 and C2 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, F1 2 a,b;c;x  is defined by analytic continuation, and for x<1, F1 2 a,b;c;x  is real and finite.
For x=1, the following apply:
In the complex plane, the principal branch of F1 2 a,b;c;z  is taken along the real axis from x=1.0 increasing. F1 2 a,b;c;z  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<x0.5, sets of safe parameters αi,j,βi,j,ζi,j, χj 1j2 ,1i4  are determined, such that the values of F 1 2 aj, bj ;cj;xj  required for an appropriate transformation of the type (3) may be calculated either directly or using recurrence relations from the solutions of F 1 2 αi,j , βi,j ;ζi,j;χj . If c is positive, then only transformations with C2=0.0 will be used, implying only F 1 2 a1, b1 ;c1;x1  will be required, with the transformed argument x1=x. If c is negative, in some cases a transformation with C20.0 will be used, with the argument x2=1.0-x. 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 xj=1.0-x is employed, and the above strategy for 0<x0.5 is utilized on both components. Further transformations in these sub-computations are however, limited to single terms with no argument transformation.
For x<0, a linear transformation mapping the argument x to the interval 0,0.5 is first employed. The strategy for 0<x0.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 0.5,1 may also be used.
For improved precision in the final result, this routine accepts a,b and c split into an integral and a decimal fractional component. Specifically, a=ai+ar, where ar0.5 and ai=a-ar is integral. The other parameters b and c are similarly deconstructed.
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. arbnd=0.0001×Imax, where Imax 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 non-finite values of the parameters and argument on entry, and assigns non-finite 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

NIST Digital Library of Mathematical Functions
Pearson J (2009) Computation of hypergeometric functions MSc Dissertation, Mathematical Institute, University of Oxford

5
Arguments

1:     ani – Real (Kind=nag_wp)Input
On entry: ai, the nearest integer to a, satisfying ai = a-ar.
Constraints:
  • ani=ani;
  • aniarbnd.
2:     adr – Real (Kind=nag_wp)Input
On entry: ar, the signed decimal remainder satisfying ar = a-ai  and ar 0.5.
Constraint: adr0.5.
3:     bni – Real (Kind=nag_wp)Input
On entry: bi, the nearest integer to b, satisfying bi = b-br.
Constraints:
  • bni=bni;
  • bniarbnd.
4:     bdr – Real (Kind=nag_wp)Input
On entry: br, the signed decimal remainder satisfying br = b-bi  and br 0.5.
Constraint: bdr0.5.
5:     cni – Real (Kind=nag_wp)Input
On entry: ci, the nearest integer to c, satisfying ci=c-cr.
Constraints:
  • cni=cni;
  • cniarbnd;
  • if cdr<16.0ε, cni1.0.
6:     cdr – Real (Kind=nag_wp)Input
On entry: cr, the signed decimal remainder satisfying cr = c-ci and cr 0.5.
Constraint: cdr0.5.
7:     x – Real (Kind=nag_wp)Input
On entry: the argument x.
Constraint: -arbnd<x1.
8:     frf – Real (Kind=nag_wp)Output
On exit: ffr, the scaled real component of the solution satisfying ffr = F 1 2 a,b;c;x × 2 -fsc , i.e., F 1 2 a,b;c;x = ffr × 2fsc. See Section 9 for the behaviour of ffr  when a finite or non-finite answer is returned.
9:     scf – IntegerOutput
On exit: fsc, the scaling power of two, satisfying fsc = log2 F 1 2 a,b;c;x ffr , i.e., F 1 2 a,b;c;x = ffr × 2fsc. See Section 9 for the behaviour of fsc when a non-finite answer is returned.
10:   ifail – IntegerInput/Output
On entry: ifail must be set to 0, -1 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 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 -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6
Error Indicators and Warnings

If on entry 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:
ifail=1
Underflow occurred during the evaluation of F12a,b;c;x. The returned value may be inaccurate.
ifail=2
All approximations have completed, and the final residual estimate indicates some precision may have been lost.
Relative residual=value.
ifail=3
All approximations have completed, and the final residual estimate indicates no accuracy can be guaranteed.
Relative residual=value.
ifail=4
On entry, x=value, c=value, a+b=value.
F12a,b;c;1 is infinite in the case ca+b.
ifail=5
On completion, overflow occurred in the evaluation of F12a,b;c;x.
ifail=6
Overflow occurred in a subcalculation of F12a,b;c;x. The answer may be completely incorrect.
ifail=9
An internal calculation has resulted in an undefined result.
ifail=11
On entry, ani does not satisfy aniarbnd=value.
ifail=13
ani is non-integral.
On entry, ani=value.
Constraint: ani=ani.
ifail=21
On entry, adr does not satisfy adr0.5.
ifail=31
On entry, bni does not satisfy bniarbnd=value.
ifail=33
bni is non-integral.
On entry, bni=value.
Constraint: bni=bni.
ifail=41
On entry, bdr does not satisfy bdr0.5.
ifail=51
On entry, cni does not satisfy cniarbnd=value.
ifail=52
On entry, c=cni+cdr=value.
F12a,b;c;x is undefined when c is zero or a negative integer.
ifail=53
cni is non-integral.
On entry, cni=value.
Constraint: cni=cni.
ifail=61
On entry, cdr does not satisfy cdr0.5.
ifail=71
On entry, x does not satisfy xarbnd=value.
ifail=72
On entry, x=value.
In general, F12a,b;c;x is not real valued when x>1.
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.
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.
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 ifail=0, the value of F 1 2 a,b;c;x  may be assumed accurate, with the possible loss of one or two decimal places. Assuming the result does not overflow, an error estimate res is made internally using equation (1). If the magnitude of this residual res is sufficiently large, a nonzero ifail will be returned. Specifically,
ifail=0 or 1 res1000ε
ifail=2 1000ε<res0.1
ifail=3 res>0.1
where ε is the machine precision as returned by x02ajf. Note that underflow may also have occurred if ifail=2 or 3.
A further estimate of the residual can be constructed using equation (1), and the differential identity,
d F 1 2 a,b;c;x dx = ab c F 1 2 a+1, b+1 ;c+1;x d2 F 1 2 a,b;c;x dx2 = aa+1 bb+1 cc+1 F 1 2 a+2, b+2 ;c+2;x (4)
This estimate is however, dependent upon the error involved in approximating F 1 2 a+1, b+1 ;c+1;x  and F 1 2 a+2, b+2 ;c+2;x .

8
Parallelism and Performance

s22bff is not threaded in any implementation.

9
Further Comments

s22bff returns non-finite values when appropriate. See Chapter X07 for more information on the definitions of non-finite values.
Should a non-finite value be returned, this will be indicated in the value of ifail, as detailed in the following cases.
If ifail=0 or ifail=1, 2 or 3, a finite value will have been returned with approximate accuracy as detailed in Section 7.
The values of ffr and fsc are implementation dependent. In most cases, if F 1 2 a,b;c;x = 0 , f fr = 0 and f sc = 0 will be returned, and if F 1 2 a,b;c;x  is finite, the fractional component will be bound by 0.5 f fr < 1, with f sc  chosen accordingly.
The values returned in frf (ffr) and scf (fsc) may be used to explicitly evaluate F 1 2 a,b;c;x , and may also be used to evaluate products and ratios of multiple values of F 1 2  as follows,
F 1 2 a,b;c;x = ffr × 2 fsc F 1 2 a1, b1 ;c1;x1 × F 1 2 a2, b2 ;c2;x2 = ffr1 × ffr2 × 2 fsc1 + fsc2 F 1 2 a1, b1 ;c1;x1 F 1 2 a2, b2 ;c2;x2 = ffr1 ffr2 × 2 fsc1 - fsc2 ln F 1 2 a,b;c;x = lnffr + fsc × ln2 .  
If ifail=4 then F1 2 a,b;c;x  is infinite. A signed infinity will have been returned for frf, and scf=0. The sign of frf should be correct when taking the limit as x approaches 1 from below.
If ifail=5 then upon completion, F1 2 a,b;c;x > 2 Imax , where Imax is given by x02bbf, and hence is too large to be representable even in the scaled form. The scaled real component returned in frf may still be correct, whilst scf=Imax will have been returned.
If ifail=6 then overflow occurred during a subcalculation of F1 2 a,b;c;x . The same result as for ifail=5 will have been returned, however, there is no guarantee that this is representative of either the magnitude of the scaling power fsc, or the scaled component ffr of F1 2 a,b;c;x .
For all other error exits, scf=0 will be returned and frf will be returned as a signalling NaN (see x07bbf).
If ifail=9 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 ifail=52 then c is too close to a negative integer or zero on entry, and F1 2 a,b;c;x  is 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 cj not being negative integers,
(ii) for any cj which remain as negative integers, one of the corresponding parameters aj or bj is a negative integer of magnitude less than cj.
In the first case, the transformation coefficients Cj aj,bj,cj,xj  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 ifail=11, 31, 51 or 71 then no computation will have been performed due to the risk of integer overflow. The actual solution may however, be finite.
ifail=72 indicates x>1, and hence the requested solution is on the boundary of the principal branch of F1 2 a,b;c;x . Hence it is multivalued, typically with a nonzero imaginary component. It is however, strictly finite.

10
Example

This example evaluates the Gauss hypergeometric function at two points in scaled form using s22bff, and subsequently calculates their product and ratio implicitly.

10.1
Program Text

Program Text (s22bffe.f90)

10.2
Program Data

None.

10.3
Program Results

Program Results (s22bffe.r)