incomplete gamma inaccurate in asymptotic range

[pascalmolin]

Currently, Gamma(0, 110) returns a 10^-4 ball (test code below).

A precise result is needed for L functions (see commit acc874e from #165) I guess one needs asymptotic expansion like http://dlmf.nist.gov/8.11#i

example code

#include "acb_hypgeom.h"

int main()
{
    slong prec = 165;
    acb_t g, s, x;
    acb_init(g);
    acb_init(s);
    acb_init(x);
    /* s = 0 */
    acb_zero(s);
    /* x = 110 */
    acb_set_si(x, 110);
    acb_hypgeom_gamma_upper(g, s, x, 2, prec);
    acb_printd(g, 30);
    flint_printf("\n");
    acb_clear(g);
    acb_clear(s);
    acb_clear(x);
    flint_cleanup();
    return EXIT_SUCCESS;
}

[fredrik-johansson]

The asymptotic 2F0 series is used when |z| is large enough, but the cutoff is not optimized. Right now, it uses the 1F1 series when the 2F0 is estimated to give less than prec bits of accuracy, even if the accuracy with the 1F1 series is much less than that (because of cancellation). The algorithm should really estimate which of the 1F1 and 2F0 formulas gives the most accuracy; something similar is done in the 1F1 code right now.

However, even with a more intelligent cutoff, there will be some loss of accuracy in the intermediate range.

For L-functions, we don't need to improve the incomplete gamma function in general; it's enough to make it more precise in this special evaluation mode.

It's very easy if s is exact: the L-functions just needs to wrap acb_hypgeom_gamma_upper in a loop that increases the working precision used to call this function until the error is less than 2^-prec (z is not exact here, but since it's just a rational multiple of pi, we can recompute it to the working precision).

If s is not exact, we could evaluate Gamma(s,z) at the midpoint and bound the propagated error.

I think we should get a good bound on |Gamma(S,z)| on a complex ball S surrounding s (to use the Cauchy integral formula) from \int_z^{\infty} |t^{S-1} e^{-t}| dt, or directly for the s-derivative from \int_z^{\infty} |t^{s-1} log(t) e^{-t}| dt.

[fredrik-johansson]

There is now a trivial solution with the new integration code.

What is less trivial is to select algorithms automatically.