Performance / internal needed precision of regularized incomplete beta function

[JeromeRF]

Huge precision is needed to compute arb_hypgeom_beta_lower(a, b, x, 1) for medium a and b. For large a and b computation is not possible.

To compute arb_hypgeom_beta_lower(10^5, 10^5, 4999/10000, 1) with a precision of 20 digits, internal precision has to be increased up to 262144 binary digits :

64 precision bits, beta(10^5, 10^5, 4999/10000, 1) = nan
...
131072 precision bits, beta(10^5, 10^5, 4999/10000, 1) = nan
262144 precision bits, beta(10^5, 10^5, 4999/10000, 1) = [0.4643650813520205199889814761064372417362 +/- 9.49e-42]

beta(10^10, 10^10, 4999/10000, 1) can not be computed.

It comes from the fact that arb_hypgeom_beta_lower relies on acb_hypgeom_beta_lower which relies on generic acb_hypgeom_2f1 implementation.

Still for positive a and b there exists an algorithm based on continued fraction expansion which does not seem to need such high internal precision. I have done a quick and dirty translation in Pari/GP of the implementation found in GNU GSL and it seems to give instantaneous correct results with Pari default precision of 38 digits :
beta.zip

beta(10^5, 10^5, 0.4999)
0.46436508135202051998898147610643731977

Larger computation can be computed as fast :

beta(10^10, 10^10, 0.4999)
2.6979112225279322861095016259203628939 E-176

This continued fraction expansion is listed on https://functions.wolfram.com/GammaBetaErf/BetaRegularized/10/ but its validity domain is not given.

[fredrik-johansson]

This problem affects all hypergeometric functions, Bessel functions, etc.

[fredrik-johansson]

Specifically for the continued fraction expansion of the incomplete beta function with real parameters and variable, I do believe that error bounds are available, but I don't have a reference at hand.

[fredrik-johansson]

This is also interesting: https://core.ac.uk/download/pdf/187723002.pdf

[JeromeRF]

Since your answer I read many papers dedicated to the computation of the incomplete beta:

• "Continued fractions for the incomplete beta function" by Aroian,
• "Continued fractions for the incomplete beta function: additions and corrections" by Marietta, Tretter and Walster.
• "Significant digit computation of the incomplete beta function ratios" by Didonato and Morris.
• "Certification of Algorithm 708: Significant-Digit computation of the incomplete beta" by Brown and Levy.

The 2nd one seems to say that convergence and error bound can be proven for a modified version of Mueller continued fraction expansion.
Still algorithm 708 described in later papers is far more complex and is split over many different cases (2000 lines of Fortran code).

[fredrik-johansson]

More or less fixed:

#include "acb_hypgeom.h"

int main()
{
    slong prec;
    acb_t z, a, b, res;

    acb_init(z);
    acb_init(a);
    acb_init(b);
    acb_init(res);

    prec = 64;
    acb_set_d(a, 1e5);
    acb_set_d(b, 1e5);
    acb_set_ui(z, 4999);
    acb_div_ui(z, z, 10000, prec);
    acb_hypgeom_beta_lower(res, a, b, z, 1, prec);
    acb_printn(res, 30, 0); printf("\n");

    prec = 128;
    acb_set_d(a, 1e10);
    acb_set_d(b, 1e10);
    acb_set_ui(z, 4999);
    acb_div_ui(z, z, 10000, prec);
    acb_hypgeom_beta_lower(res, a, b, z, 1, prec);
    acb_printn(res, 30, 0); printf("\n");

    prec = 128;
    acb_set_d(a, 1e15);
    acb_set_d(b, 1e15);
    acb_set_ui(z, 4999);
    acb_div_ui(z, z, 10000, prec);
    acb_hypgeom_beta_lower(res, a, b, z, 1, prec);
    acb_printn(res, 30, 0); printf("\n");
}

now prints:

[0.4643650813520 +/- 5.17e-14]
[2.69791122252793228610950163e-176 +/- 2.47e-203]
[1.0612052723416047758478e-17371784 +/- 7.21e-17371807]

Still, the current algorithm breaks down for parameters larger than 1e16 or so, and is not very efficient, so there is more to do.