/* gcc 2F1.c -lgsl -lgslcblas -o 2F1 -lm */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_sf_hyperg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_version.h>
extern double  my2F1(int a, int b, double c, double z);
#define max(a,b) (a>b ? a : b)

/*  hyperg_2F1(-1, -10, 1.0, 0.5) = 6 */

int main(int ac, char** av) {
  int a,b;
  double c,z,rv;
  if(ac > 1) a = atoi(av[1]); else a = -1;
  if(ac > 2) b = atoi(av[2]); else b = -10;
  if(ac > 3) c = atof(av[3]); else c =  1.0;
  if(ac > 4) z = atof(av[4]); else z =  0.5;
  rv = my2F1(a, b, c, z);
  printf("My 2F1 returns 2F1(%d, %d, %4.2f, %4.2f)=%1f.\n", a,b,c,z,rv);

  printf("Now for GSL version %s:\n", GSL_VERSION);

  rv = gsl_sf_hyperg_2F1((double)a, (double)b, c, z);
  printf("GSL version %s  returns %f.\n", GSL_VERSION, rv);
}

/* Routine to calculate 2F1 for non-positive args a,b.  Note
   same idea works if just one a,b is a non-positive integer. */

double my2F1(int a, int b, double c, double z) {
  int i, n;
  double s, t;
  n = -max(a,b);
  if(n<0) {
    fprintf(stderr, "my2F1: args a,b must be non-positive\n");
    exit(-1);
  }
  for(s=1.0, t=1.0, i=0; i<n; ++i) {
    s += ( t *= (z*(a+i)*(b+i)/((c+i)*(1+i))) );
  }
  return(s);
}