[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
## [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

**[Bug-gsl] error in hypergeometric function for large z, non-int a.**,
*Ivan Liu* **<=**