bug-gsl
[Top][All Lists]

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

 From: Ivan Liu Subject: [Bug-gsl] error in hypergeometric function for large z, non-int a. Date: Sat, 4 Jun 2005 16:22:58 +0200

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

```