Numerical evaluation of the Gauss hypergeometric function with the hypergeo package

This paper introduces the hypergeo package of R routines, for numerical calculation of hypergeometric functions. The package is focussed on efficient and accurate evaluation of the hypergeometric function over the whole of the complex plane within the constraints of fixed-precision arithmetic. The hypergeometric series is convergent only within the unit circle, so analytic continuation must be used to define the function outside the unit circle. This short document outlines the numerical and conceptual methods used in the package; and justifies the package philosophy, which is to maintain transparent and verifiable links between the software and AMS-55. The package is demonstrated in the context of game theory.


Introduction
The geometric series ∞ k=0 t k with t k = z k may be characterized by its first term and the constant ratio of successive terms t k+1 /t k = z, giving the familiar identity ∞ k=0 z k = (1 − z) −1 . Observe that while the series has unit radius of convergence, the right hand side is defined over the whole complex plane except for z = 1 where it has a pole. Series of this type may be generalized to a hypergeometric series in which the ratio of successive terms is a rational function of k: where P (k) and Q(k) are polynomials. If both numerator and denominator have been completely factored we would write t k+1 t k = (k + a 1 )(k + a 2 ) · · · (k + a p ) (k + b 1 )(k + b 2 ) · · · (k + b q )(k + 1) z (the final term in the denominator is due to historical reasons), and if we require t 0 = 1 then we write ∞ k=0 t k z k = a F b a 1 , a 2 , . . . , a p b 1 , b 2 , . . . , b q ; z (1) when defined. An absent factor is indicated with a dash; thus 0 F 0 − − ; z = e z . In most cases of interest one finds that p = 2, q = 1 suffices. Writing a, b, c for the two upper and one lower argument respectively, the resulting function 2 F 1 (a, b; c; z) is known as the hypergeometric function. Many functions of elementary analysis are of this form; examples would include logarithmic and trigonometric functions, Bessel functions, etc. For example, 2 F 1 1 2 , 1; 3 2 ; −z 2 = z −1 arctan z. Michel and Stoitsov (2008) state that physical applications are "plethora". In addition, naturally-occuring combinatorial series frequently have a sum expressible in terms of hypergeometric functions and an example from the author's work in the field game theory is given below.

Equivalent forms
The hypergeometric function's series representation, namely has unit radius of convergence by the ratio test but the integral form 3.1) due to Gauss, furnishes analytic continuation; it is usual to follow Riemann and define a cut along the positive real axis from 1 to ∞ and specify continuity from below [NB: equations with three-part numbers, as 15.1.1 and 15.3.1 above, are named for their reference in Abramowitz and Stegun (1965)]. This is implemented as f15.3.1() in the package and exhibits surprisingly accurate evaluation.
Gauss also provided a continued fraction form for the hypergeometric function [implemented as hypergeo_contfrac() in the package] which has superior convergence rates for parts of the complex plane at the expense of more complicated convergence properties (Cuyt et al. 2008).

The hypergeo package
The hypergeo package provides some functionality for the hypergeometric function; the emphasis is on fast vectorized R-centric code, complex z and moderate real values for the auxiliary parameters a, b, c. The package is released under GPL-2.
Observing the slow convergence of the series representation 15.1.1, the complex behaviour of the continued fraction representation, and the heavy computational expense of the integral representation 15.3.1, it is clear that non-trivial numerical techniques are required for a production package.
The package implements a generalization of the method of Forrey (1997) to the complex case. It utilizes the observation that the ratio of successive terms approaches z, and thus the strategy adopted is to seek a transformation which reduces the modulus of z to a minimum. Abramowitz and Stegun give the following transformations: Observing that the set z, z z−1 , 1 − z, 1 z , 1 1−z , 1 − 1 z forms a group under functional composition 1 we may apply each of the transformations to the primary argument z and choose the one of smallest absolute value to evaluate.
Given the appropriate transformation, the right hand side is evaluated using direct summation. If |z| < 1, the series is convergent by the ratio test, but may require a large number of terms to achieve acceptable numerical precision. Summation is dispatched to genhypergeo_series() which evaluates the generalized hypergeometric function 1; the R implementation uses multiplication by repeatedly incremented upper and lower indices a i , b i . (slightly edited in the interests of visual clarity). This is a typical internal function of the package and like all similar functions is named for its equation number in Abramowitz and Stegun (1965). Note the helper function i15.3.9(), which calculates the Gamma coefficients of the two hypergeometric terms in the identity. This structure allows transparent checking of the code.

Special cases
The methods detailed above are not applicable for all values of the parameters a, b, c. If, for example, c = a + b ± m, m ∈ N (a not uncommon case), then equation 15.3.6 is not useful because each term has a pole; and it is numerically difficult to approach the limit. In this case the package dispatches to hypergeo_cover1() which uses 15.3.4 through 15.3.9 but with 15.3.6 replaced with suitable limiting forms such as (15.3.11) (Abramowitz and Stegun give a similar expression for negative m). This equation is comparable to 15.3.6 in terms of computational complexity but requires evaluation of the digamma function ψ. Equation 15.3.11 is evaluated in the package using an algorithm similar to that for genhypergeo_series() but includes a runtime option which specifies whether to evaluate ψ (·) ab initio each time it is needed, or to use the recurrence relation ψ (z + 1) = ψ (z) + 1/z at each iteration after the first. These two options appear to be comparable in terms of both numerical accuracy and speed of execution, but further work would be needed to specify which is preferable in this context.
A similar methodology is used for the case b = a±m, m = 0, 1, 2, . . . in which case the package dispatches to hypergeo_cover2().
However, the case c − a = 0, 1, 2, . . . is not covered by Abramowitz and Stegun (1965) and the package dispatches to hypergeo_cover3() which uses formulae taken from the Wolfram functions site (Wolfram 2014). For example w07.23.06.0026.01() gives a straightforwardly implementable numerical expression for 2 F 1 as a sum of two finite series and a generalized hypergeometric function 3 F 2 with primary argument z −1 .
In all these cases, the limiting behaviour is problematic. For example, if a+b−c is close to, but not exactly equal to, an integer then equation 15.3.11 is not applicable. The analytic value of the hypergeometric function in these circumstances is typically of moderate modulus, but both terms of equation 15.3.6 have large amplitude and numerics are susceptible to cancellation errors.

Critical points
All the above methods fail when z = 1 2 ± i √ 3 2 , because none of the transformations 15.3.6-15.3.9 change the modulus of z from 1. The function is convergent at these points but numerical evaluation is difficult. This issue does not arise in the real case considered by Forrey (1997).
These points were considered by Buhring (1987) who presented a computational method for these values; however, his method is not suitable for finite-precision arithmetic (a brief discussion is presented at ?buhring) and the package employs either an iterative scheme due to Gosper (Johansson et al. 2010), or the residue theorem if z is close to either of these points.

Package testing suite
The package comes with an extensive test suite in the tests/ directory. The tests fall into two main categories, firstly comparison with either Maple or Mathematica output (although Becken and Schmelcher (2000) caution that Mathematica routines cannot be used as reference values); and secondly, verification of identities which appear in AMS-55 as named equations.

The package in use
The hypergeo package offers direct numerical functionality to the R user on the command line. One example from the author's current work is in game theory. Consider a game in which a player is given n counters each of which she must allocate into one of two boxes, A or B. At times t = 1, 2, 3 . . . a box is identified at random and, if it is not empty, a counter removed from it; box A is chosen with probability p and box B with probability 1 − p. The object of the game is to remove all counters as quickly as possible. If the player places a counters in box A and b in B, then the probability mass function of removing the final counter at time t = a + b + r is The two terms correspond to the final counter being removed from box A or B respectively. This PMF has expectation with R idiom: > expected <-function(a,b,p){ + Re( + choose(a+b,b) * p^a * (1-p)^b * ( + p *b/(1+a) * hypergeo(a+b+1,2,a+2, p) + + (1-p)*a/(1+b) * hypergeo(a+b+1,2,b+2,1-p) )) + } Thus if p = 0.8 and given n = 10 counters we might wonder whether it is preferable to allocate them (8, 2) or (9, 1): > c(expected(8,2,0.8),expected(9,1,0.8)) [1] 3.019899 1.921089 showing that the latter allocation is preferable in expectation.
The package is designed for use with R and Figure 1 shows the package being used to visualize 2 F 1 2, 1 2 ; 2 3 ; z over a region of the complex plane. Figure 1: View of the function 2 F 1 2, 1 2 ; 2 3 ; z evaluated over a part of the complex plane using the hypergeo package. Function visualization following Thaller (1998) and the elliptic package (Hankin 2006); hue corresponds to argument and saturation to modulus. Solid contour lines correspond to real function values and dotted to imaginary function values. Note the cut line along the real axis starting at (1, 0), made visible by an abrupt change in hue

Conclusions and further work
Evaluation of the hypergeometric function is hard, as evidenced by the extensive literature concerning its numerical evaluation (Becken and Schmelcher 2000;Michel and Stoitsov 2008;Forrey 1997;Buhring 1987). The hypergeo package is presented as a partial implementation, providing reasonably accurate evaluation over a large portion of the complex plane and covering moderate real values of the auxiliary parameters a, b, c. Difficulties arise when b − a or c − b − a become close to, but not exactly, integers because the terms in equations 15.3.6 through 15.3.9 become large and cancellation errors become important.
Further work might include extension to complex auxiliary parameters but Michel and Stoitsov caution that this is not straightforward.