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

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?

• Have a look at DLMF 15.8.2, for example. Nov 27, 2021 at 23:35

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)

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= 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.

– user20857
Dec 1, 2021 at 19:17

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
}
}