[Bug-gsl] error in hypergeometric function for large z, non-int a.

[Ivan Liu]

Hi, 
I found a strange behaviour, possibily a bug, in the routine for
hypergeometric function: gsl_sf_hyperg_1F1(a,b,z). After some
investigation I found which part in the source code is responsible for
the error.

The test code is below:
//----------- Begin code: main.cpp --------------
#include <iostream>
#include <fstream>
#include <gsl/gsl_sf_hyperg.h>
using namespace std;

int main(void)
{
  double a= -26.1; // if write instead, int a=-26, the result is fine
  int b= 2;
  ofstream file("data.dat");

  for( double z= 0.0; z<=100.0; z+= 0.1){
    double test = gsl_sf_hyperg_1F1(a,b,z);
    file << z << '\t' << test << endl;
  }
  return 0;
}
//------------- End code ---------------

Execution command:
> g++ -o run main.cpp -lgsl -lgslcblas
> run
Then plot the graph with any data plotting software, eg,
> xmgrace data.dat

Then you'll see if a is a integer the curve is smooth, but, if a is
not an integer, there is instability for large z.

In the routine, int a and non-int a are handled separately. So I guess
it's the non-int routine that gives the error. This part is from about
line 1011 in hyperg_1F1.c, beginning with
  else {
    /* x < 0
     * b < a
     */

    if(a <= 0.5*(b-x) || a >= -x) {
      /* Recurse down in b, from near the a=b line, b=a+eps,a+eps-1.
       */
      double N   = floor(a - b);
.....
until line 1577, the closure for the if. The equivalent for the case
of integer a is starts from line 1040.

I don't know how this can be fixed. Maybe someone else knows.

Best Regards,
Ivan Liu