Open JeromeRF opened 3 years ago
This problem affects all hypergeometric functions, Bessel functions, etc.
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.
This is also interesting: https://core.ac.uk/download/pdf/187723002.pdf
Since your answer I read many papers dedicated to the computation of the incomplete beta:
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).
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.
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 :
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
Larger computation can be computed as fast :
This continued fraction expansion is listed on https://functions.wolfram.com/GammaBetaErf/BetaRegularized/10/ but its validity domain is not given.