arb_hypgeom_gamma_upper(357, 356) gives nan +/- inf

[p15-git-acc]

According to Pari/GP and Mathematica it should be 1.2397102993749311592885574337209916963 E755

#include "arb_hypgeom.h"

int main()
{
    slong prec = 300;
    arb_t g, s, x;
    arb_init(g);
    arb_init(s);
    arb_init(x);
    arb_set_ui(s, 357);
    arb_set_ui(x, 356);
    arb_hypgeom_gamma_upper(g, s, x, 0, prec);
    arb_printd(g, 30);
    flint_printf("\n");
    arb_clear(g);
    arb_clear(s);
    arb_clear(x);
    flint_cleanup();
    return 0;
}

Edit: It gives the expected result if prec is increased from 300 to 400. This is probably because of conditionals in acb_hypgeom_u_asymp. For my application I think I need only a few bits of precision in the output, not anywhere near 300 or 400 bits.

[p15-git-acc]

After some more testing I think the problem is in the determination of whether or not the asymptotic expansion should be used. Here's an example with more generic input (not exact or integers or tiny intervals around integers):

#include "arb_hypgeom.h"
#include "acb_hypgeom.h"

static void
check_incomplete_gamma(const char *sstr, const char *xstr, slong prec)
{
    arb_t g, s, x;
    acb_t gc, sc, xc;

    arb_init(g);
    arb_init(s);
    arb_init(x);
    acb_init(gc);
    acb_init(sc);
    acb_init(xc);

    arb_set_str(s, sstr, prec);
    arb_set_str(x, xstr, prec);
    arb_add_error_2exp_si(s, -200);
    arb_add_error_2exp_si(x, -200);
    acb_set_arb(sc, s);
    acb_set_arb(xc, x);
    flint_printf("s: "); arb_printd(s, 20); flint_printf("\n");
    flint_printf("x: "); arb_printd(x, 20); flint_printf("\n");
    arb_hypgeom_gamma_upper(g, s, x, 0, prec);
    flint_printf("gamma_upper(s, x): ");
    arb_printd(g, 20); flint_printf("\n");
    acb_hypgeom_gamma_upper_1f1b(gc, sc, xc, 0, prec);
    flint_printf("gamma_upper_1f1b(s, x): ");
    acb_printd(gc, 20); flint_printf("\n");
    flint_printf("\n");

    arb_clear(g);
    arb_clear(s);
    arb_clear(x);
    acb_clear(gc);
    acb_clear(sc);
    acb_clear(xc);
}

int main()
{
    slong prec = 300;
    check_incomplete_gamma("356.123", "356.456", prec);
    check_incomplete_gamma("357.123", "356.456", prec);
    check_incomplete_gamma("357.456", "356.123", prec);
    flint_cleanup();
    return 0;
}

s: 356.123 +/- 6.223e-61
x: 356.456 +/- 6.223e-61
gamma_upper(s, x): 6.7788573374431732297e+752 +/- 3.1745e+58542                
gamma_upper_1f1b(s, x): (6.7788573374431732297e+752 + 0j)  +/-  (1.11e+694, 0j)

s: 357.123 +/- 6.223e-61
x: 356.456 +/- 6.223e-61
gamma_upper(s, x): 2.5190935761679531971e+755 +/- 4.7286e+3732762              
gamma_upper_1f1b(s, x): (2.5190935761679531971e+755 + 0j)  +/-  (3.89e+696, 0j)

s: 357.456 +/- 6.223e-61
x: 356.123 +/- 6.223e-61
gamma_upper(s, x): nan +/- inf
gamma_upper_1f1b(s, x): (1.832655185502837768e+756 + 0j)  +/-  (2.72e+697, 0j) 

[fredrik-johansson]

Same issue as #166

Note that increasing the precision adaptively is a (not necessarily convenient) workaround.

[p15-git-acc]

From latest DLMF: "Sharp error bounds and an exponentially-improved extension for (8.11.7) can be found in Nemes (2016)".