[Bug-gsl]gamma_inc_P

[Martin Jansche]

Hi,

I've run into trouble caused either by gsl_sf_gamma_inc_P() or by my
lack of understanding of how it is supposed to behave.  I'm using GSL
1.3 on Debian "Sid" x86.  Here's a program that illustrates the
problem:

#include <stdio.h>
#include <gsl/gsl_sf_gamma.h>

int main(void) {
  double x = 2.77556e-16;
  double b = 2.17;
  double y = x / b;
  printf("%.16g\n", y);
  double d = gsl_sf_gamma_inc_P(10.0, y);
  return 0;
}

It generates the following output:

$ ./a.out
1.279059907834101e-16
gsl: log.c:211: ERROR: domain error
Default GSL error handler invoked.
Aborted
$

The code in specfunc/gamma_inc.c and the documentation agree on what
the domain is (a > 0, x >= 0), so the propagated domain error comes as
a surprise to me.  What am I missing?

Thanks,

- martin