The Gauss hypergeometric function

The flags ACB_HYPGEOM_2F1_AB (a - b \in \mathbb{Z}) and ACB_HYPGEOM_2F1_ABC (a + b - c \in \mathbb{Z}) are used in this way; AC and BC are accepted but currently unused. For example, to evaluate {}_2F_1(\sqrt{2}, 1/2, \sqrt{2}+3/2, 9/10), set ACB_HYPGEOM_2F1_ABC; otherwise the transformation z \to 1 - z divides by \sin(\pi(c-a-b)) computed on a ball containing zero and the result is indeterminate. If the parameters are exact (including exact integers and half-integers), the limits are detected automatically and the flags are unnecessary.

The function handles special parameters first (terminating series, z = 0 and z = 1, a = c or b = c for \mathbf{F}, reduction to a polynomial by the Pfaff transformation when c - a or c - b is a nonpositive integer), then chooses with acb_hypgeom_2f1_choose between direct summation, one of five transformations and the corner case, and retries real arguments with arb_hypgeom_2f1_integration if the result is inaccurate. See The Gauss hypergeometric function and The Gauss hypergeometric function.

acb_hypgeom_2f1 computes the Gauss hypergeometric function {}_2F_1(a,b,c,z), or the regularized function \mathbf{F}(a,b,c,z) = {}_2F_1(a,b,c,z)/\Gamma(c) if the flag ACB_HYPGEOM_2F1_REGULARIZED is set, for all complex a, b, c, z, with the branch cut [1, +\infty). arb_hypgeom_2f1 is the real version and accepts the same flags. The series converges only for |z| < 1, and slowly near |z| = 1; the function covers the plane with five linear fractional transformations, evaluates the limits of these formulas at integer parameter differences, continues analytically along a path near the two points e^{\pm\pi i/3} that no transformation reaches, and falls back to numerical integration for real arguments. The algorithms are described in The Gauss hypergeometric function.

Algorithm selection#

1. Non-finite input gives an indeterminate result; z = 0 gives 1 (or 1/\Gamma(c)).
2. For \mathbf{F} with a nonpositive integer c, the result is 0 if a or b is a nonpositive integer at least c (the series terminates before the poles). If a = c or b = c, the regularized function is (1-z)^{-b}/\Gamma(c) or (1-z)^{-a}/\Gamma(c).
3. If a or b is a nonpositive integer -n with n < \mathit{prec}, the function is a polynomial of degree n, summed directly for any z.
4. If c is exact and c - b (or c - a) is a nonpositive integer, the Pfaff transformation F(a,b,c,z) = (1-z)^{-a} F(a, c-b, c, z/(z-1)) turns the function into a polynomial.
5. At z = 1, Gauss's formula \Gamma(c)\Gamma(c-a-b)/(\Gamma(c-a) \Gamma(c-b)) is used if \operatorname{Re}(c-a-b) > 0; otherwise the result is indeterminate.
6. Otherwise acb_hypgeom_2f1_choose selects direct summation, one of the transformations z/(z-1), 1/z, 1/(1-z), 1-z, 1-1/z of acb_hypgeom_2f1_transform, or the corner case acb_hypgeom_2f1_corner, by minimizing the modulus of the new argument (direct summation is preferred for |z| \le 1/4).
7. If the result has less than half the requested accuracy and less than half the accuracy of the inputs, and all arguments are real and smaller than 2^{60}, arb_hypgeom_2f1_integration is tried and the better enclosure is returned.

Integer parameter differences#

The transformations to 1/z and 1/(1-z) divide by \sin(\pi(b-a)), and those to 1-z and 1-1/z by \sin(\pi(c-a-b)), which vanish when the parameter difference is an integer; the function then computes the limit (by differentiating with respect to a perturbation of a). It recognizes an integer difference only if it can prove it from the balls, which works for exact parameters but not, for instance, for a = \sqrt 2 and c = \sqrt 2 + 3/2, whose computed difference is an inexact ball containing 3/2. The flags ACB_HYPGEOM_2F1_AB and ACB_HYPGEOM_2F1_ABC declare that a - b or a + b - c is an exact integer; the caller is responsible for this being true. The flags ACB_HYPGEOM_2F1_AC and ACB_HYPGEOM_2F1_BC are accepted but currently not used.

Example#

With c - a - b = 1 known only as a ball, the transformation to 1-z needed at z = 0.9 fails unless the flag is given; near e^{\pi i/3} the corner case (method 6) is used; on the branch cut, the value is the limit from the lower half-plane:

#include <flint/acb.h>
#include <flint/acb_hypgeom.h>

int main(void)
{
    acb_t a, b, c, z, r;

    acb_init(a);
    acb_init(b);
    acb_init(c);
    acb_init(z);
    acb_init(r);

    /* a = sqrt(2) + i/2, b = 1/2, c = a + b + 1: c - a - b = 1 exactly,
       but not as balls */
    acb_set_ui(a, 2);
    acb_sqrt(a, a, 64);
    arb_set_d(acb_imagref(a), 0.5);
    acb_set_d(b, 0.5);
    acb_add_ui(c, a, 1, 64);
    acb_add(c, c, b, 64);
    acb_set_ui(z, 9);
    acb_div_ui(z, z, 10, 64);

    acb_hypgeom_2f1(r, a, b, c, z, 0, 64);
    flint_printf("no flags : "); acb_printn(r, 15, 0); flint_printf("\n");
    acb_hypgeom_2f1(r, a, b, c, z, ACB_HYPGEOM_2F1_ABC, 64);
    flint_printf("ABC flag : "); acb_printn(r, 15, 0); flint_printf("\n");

    /* near exp(pi i/3) none of the transformations converges well */
    acb_set_d(a, 0.25);
    acb_set_d(b, 1.5);
    acb_set_ui(c, 3);
    acb_set_d_d(z, 0.5, 0.875);
    flint_printf("method   : %d\n", acb_hypgeom_2f1_choose(z));
    acb_hypgeom_2f1(r, a, b, c, z, 0, 64);
    flint_printf("corner   : "); acb_printn(r, 15, 0); flint_printf("\n");

    /* on the branch cut (1, +inf): the limit from the lower half-plane */
    acb_set_ui(z, 3);
    acb_hypgeom_2f1(r, a, b, c, z, 0, 64);
    flint_printf("z = 3    : "); acb_printn(r, 15, 0); flint_printf("\n");

    acb_clear(a);
    acb_clear(b);
    acb_clear(c);
    acb_clear(z);
    acb_clear(r);
    flint_cleanup();
    return 0;
}

no flags : nan + nan*I
ABC flag : [1.45790187393845 +/- 4.64e-15] + [0.110163228098910 +/- 3.70e-16]*I
method   : 6
corner   : [1.01464883732715 +/- 1.70e-15] + [0.134531649059307 +/- 4.70e-16]*I
z = 3    : [0.995480256475680 +/- 1.90e-16] + [-0.581843941733150 +/- 3.74e-16]*I

For real parameters, the first call would not fail: acb_hypgeom_2f1 would fall back to arb_hypgeom_2f1_integration, which needs no transformation. With a complex parameter, as here, only the flag helps.