This routine computes the normalized incomplete Beta function
\(I_x(a,b) = \operatorname{B}_x(a,b)/\operatorname{B}(a,b)\) where
\(\operatorname{B}_x(a,b) = \int_0^x t^{a-1} (1-t)^{b-1} dt\) for
\(0 \leq x \leq 1\). For \(a > 0\), \(b > 0\) the value is computed using
a continued fraction expansion. For all other values it is computed
using the relation
\[I_x(a,b) = (1/a) x^a {}_2F_1(a,1-b,a+1,x)/\operatorname{B}(a,b).\]