[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Bug-gsl] overflow in gsl_ran_beta_pdf, gsl 1.14
From: |
Brian Gough |
Subject: |
Re: [Bug-gsl] overflow in gsl_ran_beta_pdf, gsl 1.14 |
Date: |
Wed, 21 Jul 2010 21:13:02 +0100 |
User-agent: |
Wanderlust/2.15.6 (Almost Unreal) Emacs/23.2 Mule/6.0 (HANACHIRUSATO) |
At Thu, 15 Jul 2010 15:17:00 -0700,
Michael Radford wrote:
> 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.
>
Thanks for the bug report and patch! I've applied it. Sorry for the
delay in replying.
--
Brian Gough