invIncompleteBeta fails for some (exotic) inputs

[Shimuuar]

*Numeric.SpecFunctions.Internal> invIncompleteBeta 9 2 1e-300
*** Exception: Numeric.SpecFunctions.incompletBeta_: x out of [0,1] range. p=9.0 q=2.0 x=NaN

It fails in Halley step when derivative of incomplete beta underflows and becomes zero so no next approximation could be computed, It happens for large a and small x.

[Shimuuar]

Problem is underflow in derivative of incomplete beta which becomes zero and next Halley correction becomes NaN. Naive solution is to revert to bisection. But convergence is slow so function gives answers which are several orders of magnitude off.

Attempted fix:

      -- When derivative is equal to zero (actually underflows since
      -- it's always positive) we cannot make Newton iteration so we
      -- have to switch to bisection
      | f' == 0                      = case () of
        _| f > 0                     -> loop (i+1) (x / 2)
         | f == 0                    -> x
         | otherwise                 -> loop (i+1) ((1 + x) / 2)
      -- When derivative becomes infinite we cannot continue