[Rd] qbeta hang (PR#2894)

[Morten Welinder]

Ok, I can confirm that it does not, in fact, loop forever.  Just a close
approximation.

> Morten: the gcc version is often crucial in these cases.

Sorry.  gcc-2.7.2 (both -O2 and -O0) on Solaris 2.9 / sparc.

The problems seem to be...

1. The initial guess underflows to zero.

2. That guess is replaced by the truly awful guess 0.5.

3. The root finding algorithm ignores all the nice properties of pbeta
   -- such as it being monotonically increasing and bounded.

> and you're trying to solve pbeta(x, 1e-8, 0.5,, TRUE, FALSE) == 0.1

Not quite.  What I was trying to do was the see what would happen with an
insanely small alpha.  Actually doing that would involve getting the
arguments right, of course, :-)  The not-quite-hang was an unplesant side
effect of swapping the args.

I notice things like

    r = sqrt(-log(a * a));

in the code.  I fail to see why that isn't just

    r = sqrt(-2 * log (a))

which ought to be faster and more accurate.  Then later

    ...log((1. - a) * qq)...

which might be better as

    ...(log1p(-a) + log (qq))...

if a is close to zero.

There are lots of other places that worry me with respect to cancellation
errors, for example

    r = 1 - pp;
    t = 1 - qq;

Morten