diff --git a/acb_dirichlet/platt_c_bound.c b/acb_dirichlet/platt_c_bound.c index ffe9e1e5..bfea8874 100644 --- a/acb_dirichlet/platt_c_bound.c +++ b/acb_dirichlet/platt_c_bound.c @@ -12,33 +12,6 @@ #include "acb_dirichlet.h" #include "arb_hypgeom.h" -/* Increase precision adaptively. */ -static void -_gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z, - int regularized, slong prec) -{ - if (!arb_is_finite(s) || !arb_is_finite(z)) - { - arb_indeterminate(res); - } - else - { - arb_t x; - slong i; - arb_init(x); - for (i = 0; i < 5; i++) - { - arb_hypgeom_gamma_upper(x, s, z, regularized, prec << i); - if (arb_rel_accuracy_bits(x) > 1) - { - break; - } - } - arb_swap(res, x); - arb_clear(x); - } -} - static void _arb_pow_si(arb_t res, const arb_t x, slong y, slong prec) { @@ -175,7 +148,7 @@ _pre_c_p(arb_ptr res, slong sigma, const arb_t h, ulong k, slong prec) arb_set_si(x, k + l + 1); arb_mul_2exp_si(x, x, -1); - _gamma_upper_workaround(x, x, u, 0, prec); + arb_hypgeom_gamma_upper(x, x, u, 0, prec); arb_mul(res + l, res + l, x, prec); } diff --git a/acb_dirichlet/platt_lemma_A11.c b/acb_dirichlet/platt_lemma_A11.c index bad7fa8b..f69009d1 100644 --- a/acb_dirichlet/platt_lemma_A11.c +++ b/acb_dirichlet/platt_lemma_A11.c @@ -12,33 +12,6 @@ #include "acb_dirichlet.h" #include "arb_hypgeom.h" -/* Increase precision adaptively. */ -static void -_gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z, - int regularized, slong prec) -{ - if (!arb_is_finite(s) || !arb_is_finite(z)) - { - arb_indeterminate(res); - } - else - { - arb_t x; - slong i; - arb_init(x); - for (i = 0; i < 5; i++) - { - arb_hypgeom_gamma_upper(x, s, z, regularized, prec << i); - if (arb_rel_accuracy_bits(x) > 1) - { - break; - } - } - arb_swap(res, x); - arb_clear(x); - } -} - static void _platt_lemma_A11_X(arb_t out, const arb_t t0, const arb_t h, slong B, const arb_t beta, slong prec) @@ -90,7 +63,7 @@ _platt_lemma_A11_Y(arb_t out, arb_sqr(x4, x4, prec); arb_mul_2exp_si(x4, x4, -3); - _gamma_upper_workaround(x5, x3, x4, 0, prec); + arb_hypgeom_gamma_upper(x5, x3, x4, 0, prec); arb_mul(out, x1, x2, prec); arb_mul(out, out, x5, prec); @@ -131,7 +104,7 @@ _platt_lemma_A11_Z(arb_t out, arb_sqr(x4, x4, prec); arb_mul_2exp_si(x4, x4, -1); - _gamma_upper_workaround(x5, x3, x4, 0, prec); + arb_hypgeom_gamma_upper(x5, x3, x4, 0, prec); arb_mul(out, x1, x2, prec); arb_mul(out, out, x5, prec); diff --git a/acb_dirichlet/platt_lemma_A5.c b/acb_dirichlet/platt_lemma_A5.c index c8f30cf4..6e17d16c 100644 --- a/acb_dirichlet/platt_lemma_A5.c +++ b/acb_dirichlet/platt_lemma_A5.c @@ -12,34 +12,6 @@ #include "acb_dirichlet.h" #include "arb_hypgeom.h" - -/* Increase precision adaptively. */ -static void -_gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z, - int regularized, slong prec) -{ - if (!arb_is_finite(s) || !arb_is_finite(z)) - { - arb_indeterminate(res); - } - else - { - arb_t x; - slong i; - arb_init(x); - for (i = 0; i < 5; i++) - { - arb_hypgeom_gamma_upper(x, s, z, regularized, prec << i); - if (arb_rel_accuracy_bits(x) > 1) - { - break; - } - } - arb_swap(res, x); - arb_clear(x); - } -} - /* Lemma A5 requires B > h*sqrt(k) */ static int _platt_lemma_A5_constraint(slong B, const arb_t h, slong k, slong prec) @@ -106,7 +78,7 @@ acb_dirichlet_platt_lemma_A5(arb_t out, slong B, const arb_t h, arb_mul_2exp_si(x6, b, -1); - _gamma_upper_workaround(x6, x6, a, 0, prec); + arb_hypgeom_gamma_upper(x6, x6, a, 0, prec); arb_mul(out, x4, x5, prec); arb_mul(out, out, x6, prec); diff --git a/acb_dirichlet/platt_ws_interpolation.c b/acb_dirichlet/platt_ws_interpolation.c index d93aa993..be28c285 100644 --- a/acb_dirichlet/platt_ws_interpolation.c +++ b/acb_dirichlet/platt_ws_interpolation.c @@ -13,33 +13,6 @@ #include "arb_hypgeom.h" #include "arb_poly.h" -/* Increase precision adaptively. */ -static void -_gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z, - int regularized, slong prec) -{ - if (!arb_is_finite(s) || !arb_is_finite(z)) - { - arb_indeterminate(res); - } - else - { - arb_t x; - slong i; - arb_init(x); - for (i = 0; i < 5; i++) - { - arb_hypgeom_gamma_upper(x, s, z, regularized, prec << i); - if (arb_rel_accuracy_bits(x) > 1) - { - break; - } - } - arb_swap(res, x); - arb_clear(x); - } -} - static void _arb_div_si_si(arb_t res, slong x, slong y, slong prec) { @@ -218,7 +191,7 @@ _platt_bound_C3_Y(arb_t res, const arb_t t0, slong A, const arb_t H, arb_div(g2, g2, H, prec); arb_sqr(g2, g2, prec); arb_mul_2exp_si(g2, g2, -1); - _gamma_upper_workaround(g, g1, g2, 0, prec); + arb_hypgeom_gamma_upper(g, g1, g2, 0, prec); /* res = a*b*A*H*g */ arb_mul_si(res, H, A, prec); @@ -261,7 +234,7 @@ _platt_bound_C3_Z(arb_t res, const arb_t t0, slong A, const arb_t H, arb_div(g2, t0, H, prec); arb_sqr(g2, g2, prec); arb_mul_2exp_si(g2, g2, -1); - _gamma_upper_workaround(g, g1, g2, 0, prec); + arb_hypgeom_gamma_upper(g, g1, g2, 0, prec); /* res = a*b*A*g */ arb_mul_si(res, g, A, prec); diff --git a/acb_hypgeom/gamma_upper.c b/acb_hypgeom/gamma_upper.c index 5f182bbd..68f3eaef 100644 --- a/acb_hypgeom/gamma_upper.c +++ b/acb_hypgeom/gamma_upper.c @@ -398,7 +398,9 @@ acb_hypgeom_gamma_upper(acb_t res, const acb_t s, const acb_t z, int regularized c = 2.391097143; ds = arf_get_d(arb_midref(acb_realref(s)), ARF_RND_DOWN); dx = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN); - if (pow(dx - c, 8) > pow(a*ds, 8) + pow(b*prec, 8)) + if (dx - c > a*ds + b*prec || + (dx - c > FLINT_MAX(a*ds, b*prec) && + pow(dx - c, 8) > pow(a*ds, 8) + pow(b*prec, 8))) { acb_hypgeom_gamma_upper_asymp(res, s, z, regularized, prec); return; diff --git a/acb_hypgeom/test/t-gamma_upper.c b/acb_hypgeom/test/t-gamma_upper.c index 0dcbcb00..fef392b8 100644 --- a/acb_hypgeom/test/t-gamma_upper.c +++ b/acb_hypgeom/test/t-gamma_upper.c @@ -11,6 +11,28 @@ #include "acb_hypgeom.h" + +static void +_accuracy_regression_test(const acb_t s, const acb_t z, + int regularized, slong prec, slong issue, slong accuracy) +{ + acb_t g; + acb_init(g); + acb_hypgeom_gamma_upper(g, s, z, regularized, prec); + if (acb_rel_accuracy_bits(g) < accuracy) + { + flint_printf("FAIL: accuracy regression in issue #%ld\n\n", issue); + flint_printf("prec = %d\n\n", prec); + flint_printf("regularized = %d\n\n", regularized); + flint_printf("s = "); acb_printd(s, 30); flint_printf("\n\n"); + flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n"); + flint_printf("g = "); acb_printd(g, 30); flint_printf("\n\n"); + flint_abort(); + } + acb_clear(g); +} + + int main() { slong iter; @@ -211,6 +233,40 @@ int main() acb_clear(u); } + /* Accuracy regression tests. */ + { + acb_t s, z; + slong prec, issue, accuracy; + acb_init(s); + acb_init(z); + + issue = 166; + prec = 165; + accuracy = 100; + acb_zero(s); + acb_set_si(z, 110); + _accuracy_regression_test(s, z, 2, prec, issue, accuracy); + + issue = 276; + prec = 300; + accuracy = 100; + acb_set_ui(s, 357); + acb_set_ui(z, 356); + _accuracy_regression_test(s, z, 0, prec, issue, accuracy); + arb_set_str(acb_realref(s), "356.123", prec); + arb_set_str(acb_realref(z), "356.456", prec); + _accuracy_regression_test(s, z, 0, prec, issue, accuracy); + arb_set_str(acb_realref(s), "357.123", prec); + arb_set_str(acb_realref(z), "356.456", prec); + _accuracy_regression_test(s, z, 0, prec, issue, accuracy); + arb_set_str(acb_realref(s), "357.456", prec); + arb_set_str(acb_realref(z), "356.123", prec); + _accuracy_regression_test(s, z, 0, prec, issue, accuracy); + + acb_clear(s); + acb_clear(z); + } + flint_randclear(state); flint_cleanup(); flint_printf("PASS\n");