Putting hypergeometric function in gnuplot with GSL (patching)

My gnuplot version with hypergeometric function using GNU Scientific Library
hyperg_2F1.patch

--- ./src/eval.c.orig    2016-08-12 17:06:35.250145505 -0300
+++ ./src/eval.c    2016-05-03 12:34:08.000000000 -0300
@@ -191,6 +191,7 @@
     {"lambertw",  f_lambertw}, /* HBB, from G.Kuhnle 20001107 */
     {"airy",  f_airy},         /* janert, 20090905 */
     {"expint",  f_expint},     /* Jim Van Zandt, 20101010 */
+    {"hyperg_2F1",f_hyperg_2F1},

 #ifdef HAVE_LIBCERF
     {"cerf", f_cerf},        /* complex error function */
--- ./src/specfun.c.orig    2016-08-12 17:06:35.248145450 -0300
+++ ./src/specfun.c    2016-05-03 12:20:14.000000000 -0300
@@ -49,6 +49,7 @@
 #include "specfun.h"
 #include "stdfn.h"
 #include "util.h"
+#include "gsl/gsl_sf_hyperg.h"

 #define ITMAX   200

@@ -3031,9 +3032,11 @@
 extern double p1evl ( double, void *, int );
 extern double sin ( double );
 extern double cos ( double );
+extern double gsl_sf_hyperg_2F1( double , double, double, double);
 #else
 double fabs(), exp(), sqrt();
 double polevl(), p1evl(), sin(), cos();
+double gsl_sf_hyperg_2F1();
 #endif

 int
--- ./src/standard.c.orig    2016-08-12 17:06:35.255145644 -0300
+++ ./src/standard.c    2016-05-03 14:40:24.741731493 -0300
@@ -39,6 +39,7 @@
 #include "gadgets.h"        /* for 'ang2rad' and 'zero' */
 #include "gp_time.h"        /* needed by f_tmsec() and friendsa */
 #include "util.h"        /* for int_error() */
+#include "gsl/gsl_sf_hyperg.h"

 static double carlson_elliptic_rc(double x,double y);
 static double carlson_elliptic_rf(double x,double y,double z);
@@ -57,7 +58,7 @@
 static double yone __PROTO((double x));
 static double rj1 __PROTO((double x));
 static double ry1 __PROTO((double x));
-
+extern double gsl_sf_hyperg_2F1(double,double,double,double);
 /* The bessel function approximations here are from
  * "Computer Approximations"
  * by Hart, Cheney et al.
@@ -841,6 +842,30 @@
     }
 }

+void
+f_hyperg_2F1(union argument *arg)
+{
+    struct value a1;
+    double    aa,bb,cc,zz;
+
+    (void) arg;            /* avoid -Wunused warning */
+    zz=real(pop(&a1));
+    cc=real(pop(&a1));
+    bb=real(pop(&a1));
+    aa=real(pop(&a1));
+//    if (fabs(imag(&a1)) > zero || fabs(imag(&a2)) > zero)
+//    int_error(NO_CARET, "can only do elliptic integrals of reals");
+
+    if (zz < 1.0 && zz > 0)
+    push(Gcomplex(&a1, gsl_sf_hyperg_2F1(aa,bb,cc,zz), 0.0));
+    else if (zz < 0) {
+    push(Gcomplex(&a1, pow(1-zz,-aa)*gsl_sf_hyperg_2F1(aa,cc-bb,cc,zz/(zz-1)), 0.0));
+    } else {
+    undefined=TRUE;
+    push(&a1);
+    }
+
+}
 /* Terminate the autoconversion from string to numeric values */
 #undef pop

--- ./src/standard.h.orig    2016-08-12 17:06:35.252145560 -0300
+++ ./src/standard.h    2016-05-03 12:37:23.000000000 -0300
@@ -93,5 +93,6 @@
 void f_tmyear __PROTO((union argument *x));
 void f_tmwday __PROTO((union argument *x));
 void f_tmyday __PROTO((union argument *x));
+void f_hyperg_2F1 __PROTO((union argument *x));

 #endif /* GNUPLOT_STANDARD_H */