[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

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

reply via email to

[Prev in Thread] Current Thread [Next in Thread]