Examining bug #21838

[Brian Gough]

At Fri, 15 Feb 2008 00:00:55 +0100,
Frank Reininghaus wrote:
> I examined the first testcase: It seems that the origin of the NaN is in 
> beta_cont_frac () in the source file cdf/beta_inc.c. The loop aborts 
> after the maximum number of iterations because fabs (delta_frac - 1.0) 
> does not get small enough, and NaN is returned.
> 
> The funny thing is that I could influence this behaviour in a way I 
> didn't expect:

That's interesting -- it looks like there is a problem with extended
precision making the result worse.
> 
> (P - eP) [0] = nan, tolerance = 2.3283064365386963e-10
> (P - eP) [1] = nan, tolerance = 2.3283064365386963e-10
> (P - eP) [2] = nan, tolerance = 2.3283064365386963e-10
> (P - eP) [3] = 4.1520314247533996e-05, tolerance = 2.3283064365386963e-10

Here is a corrrection for the final value, it should be:

   eP [3] = 0.9167386972447996480399

The results should make more sense with that.  I'll commit that in the
test suite.

-- 
Brian Gough