GSL 不这样做的原因是您需要分析延续。这就是 Mathematica 所做的。例如，请参见此处和此处。这里提出了一个非常相似的问题，但他们提出了一个自己动手的解决方案。

编辑：根据我上面的链接在此处提供有用的论文参考（感谢您将其拉出@JoseM）

请将其视为玩具或教育代码...如果您想自己做，请考虑使用诸如 GMP 之类的可变精度库。

在我的情况下 z > 2，所以阅读@spencer-bryngelson 建议的论坛上的论文，超几何函数的计算，你会得到 Eq.4.20（使用表 13），或者你可以按照@a-rural 的建议使用DLMF 15.8.2 -reader（对我来说这更清楚，因为如果也看到Hypergeometric2F1 ，你可以看到不连续性）

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

如您所见，结果远非准确。

如果你不想自己做，我认为Arb似乎是一个不错的选择。

您可以使用转换公式。这是我在 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
    }
}