[Bug-gsl] overflow in gsl_ran_beta_pdf, gsl 1.14

[Michael Radford]

Hi,

The function gsl_ran_beta_pdf contains the following section,
which is causing incorrect overflows for x == 0 or 1 and when
a and b are large enough:

(at randist/beta.c:53)
      ...
      double gab = gsl_sf_lngamma (a + b);
      double ga = gsl_sf_lngamma (a);
      double gb = gsl_sf_lngamma (b);

      if (x == 0.0 || x == 1.0)
        {
          p = exp (gab - ga - gb) * pow (x, a - 1) * pow (1 - x, b - 1);
        }
      else
      ...

For example:
(gdb) call gsl_ran_beta_pdf(1,100,100)
$26 = 0
(gdb) call gsl_ran_beta_pdf(0,100,100)
$27 = 0
(gdb) call gsl_ran_beta_pdf(1,1000,1000)
$28 = -nan(0x8000000000000)
(gdb) call gsl_ran_beta_pdf(0,1000,1000)
$29 = -nan(0x8000000000000)

but the correct result is always 0 when x == 0 or 1 and a, b > 1.

The reason is that the call to exp() overflows when a and b are
large enough.  However, it's unnecessary to call exp() when either
of the two calls to pow() will return 0.

I'm attaching a minimal patch which simply checks for a, b > 1
and explicitly returns 0 in that case.

Mike

beta.patch
Description: Text Data