Hypergeometric function $_2F_1(z)$ with $|z| > 1$ in GSL

Question

I need to evaluate the hypergeometric function $_2F_1$ with $|z| > 1$ as in Wolfram Language with GSL but the GSL documentation says the $_2F_1$ needs $|z| < 1$.

Is there any way I can use GSL and $_2F_1$ with $|z| > 1$? Any reference about the subject?

Answer 1

The reason GSL doesn't do this is that you need an analytic continuation. That's what Mathematica does. For example, see here and here. A very similar question was asked here, but they propose a roll-your-own solution.

Edit: Useful paper reference here per my link above (thanks for pulling it out @JoseM)

Answer 2

Please, consider this as a toy or educational code... If you want to do it yourself consider the use of a variable precision library like GMP.

In my case z > 2, so reading the paper on the forum suggested by @spencer-bryngelson, Computation of Hypergeometric Functions you get Eq.4.20 (use table 13) or you can use DLMF 15.8.2 as suggested by @a-rural-reader (for me this's more clear because you can see a discontinuity if $$a - b \in \mathbb{Z}$$ see Hypergeometric2F1 too)

#include <stdio.h>
#include <stdlib.h>

#include <gsl_sf.h>
#include <gsl_math.h>

double _2F1(double a,
            double b,
            double c,
            double z)
{
  /* a - b !€ Z */
  double bd = b - 1e-15;

  double A = (gsl_sf_gamma(bd - a)*gsl_sf_gamma(c)*pow(-z, -a))/(gsl_sf_gamma(b)*gsl_sf_gamma(c - a));
  double Fa = gsl_sf_hyperg_2F1(a, a - c + 1.0, a - bd + 1.0, 1.0/z);
  double B = (gsl_sf_gamma(a - bd)*gsl_sf_gamma(c)*pow(-z, -b))/(gsl_sf_gamma(a)*gsl_sf_gamma(c - bd));
  double Fb = gsl_sf_hyperg_2F1(bd, bd - c + 1.0, bd - a + 1.0, 1.0/z);

  return A*Fa + B*Fb;
}

int main()
{
  /* Mathematica
  Hypergeometric2F1[2., 3., 4., 5.0]
  = 0.156542 + 0.150796i
  */

  double a = 2.0;
  double b = 3.0;
  double c = 4.0;
  double z = 5.0;

  double result = _2F1(a, b, c, z);
  printf("2F1(%.1f, %.1f, %.1f, %.1f): %f\n", a, b, c, z, result);
  /* result = 0.15625*/

  /*
  Hypergeometric2F1[2., 3., 4., 10.0]
  Out[30]= 0.03985 + 0.0188496 I
  */
  z = 10.0;

  result = _2F1(a, b, c, z);
  printf("2F1(%.1f, %.1f, %.1f, %.1f): %f\n", a, b, c, z, result);
  /* result = 0.0341796875*/

  return 0;
}

As you can see the results are far from accurate.

If you don't want to do it yourself, I think Arb seems a good choice.

Answer 3

You can use a transformation formula. Here is how I do in R:

library(gsl)
Gauss2F1 <- function(a,b,c,x){
    if(x>=0 & x<1){
        hyperg_2F1(a,b,c,x)
    }else{
        hyperg_2F1(c-a,b,c,1-1/(1-x))/(1-x)^b
    }
}