From 90ac72b841acd232c13a3b54f520dfdb12b941f8 Mon Sep 17 00:00:00 2001
From: p15-git-acc <37548430+p15-git-acc@users.noreply.github.com>
Date: Mon, 10 Jun 2019 16:13:14 -0500
Subject: [PATCH] upper incomplete gamma workaround

---
 acb_dirichlet/platt_c_bound.c          | 29 +++++++++++++++++++++++-
 acb_dirichlet/platt_lemma_A11.c        | 31 ++++++++++++++++++++++++--
 acb_dirichlet/platt_lemma_A5.c         | 29 +++++++++++++++++++++++-
 acb_dirichlet/platt_ws_interpolation.c | 31 ++++++++++++++++++++++++--
 4 files changed, 114 insertions(+), 6 deletions(-)

diff --git a/acb_dirichlet/platt_c_bound.c b/acb_dirichlet/platt_c_bound.c
index bfea8874..107aca1b 100644
--- a/acb_dirichlet/platt_c_bound.c
+++ b/acb_dirichlet/platt_c_bound.c
@@ -12,6 +12,33 @@
 #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 * (1 << 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)
 {
@@ -148,7 +175,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);
-        arb_hypgeom_gamma_upper(x, x, u, 0, prec);
+        _gamma_upper_workaround(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 f69009d1..18eaaa4b 100644
--- a/acb_dirichlet/platt_lemma_A11.c
+++ b/acb_dirichlet/platt_lemma_A11.c
@@ -12,6 +12,33 @@
 #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 * (1 << 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)
@@ -63,7 +90,7 @@ _platt_lemma_A11_Y(arb_t out,
     arb_sqr(x4, x4, prec);
     arb_mul_2exp_si(x4, x4, -3);
 
-    arb_hypgeom_gamma_upper(x5, x3, x4, 0, prec);
+    _gamma_upper_workaround(x5, x3, x4, 0, prec);
 
     arb_mul(out, x1, x2, prec);
     arb_mul(out, out, x5, prec);
@@ -104,7 +131,7 @@ _platt_lemma_A11_Z(arb_t out,
     arb_sqr(x4, x4, prec);
     arb_mul_2exp_si(x4, x4, -1);
 
-    arb_hypgeom_gamma_upper(x5, x3, x4, 0, prec);
+    _gamma_upper_workaround(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 54308aca..81f515b9 100644
--- a/acb_dirichlet/platt_lemma_A5.c
+++ b/acb_dirichlet/platt_lemma_A5.c
@@ -13,6 +13,33 @@
 #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 * (1 << 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)
@@ -79,7 +106,7 @@ acb_dirichlet_platt_lemma_A5(arb_t out, slong B, const arb_t h,
 
     arb_mul_2exp_si(x6, b, -1);
 
-    arb_hypgeom_gamma_upper(x6, x6, a, 0, prec);
+    _gamma_upper_workaround(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 79ff15cf..0c14aa12 100644
--- a/acb_dirichlet/platt_ws_interpolation.c
+++ b/acb_dirichlet/platt_ws_interpolation.c
@@ -12,6 +12,33 @@
 #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 * (1 << 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)
 {
@@ -188,7 +215,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);
-    arb_hypgeom_gamma_upper(g, g1, g2, 0, prec);
+    _gamma_upper_workaround(g, g1, g2, 0, prec);
 
     /* res = a*b*A*H*g */
     arb_mul_si(res, H, A, prec);
@@ -231,7 +258,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);
-    arb_hypgeom_gamma_upper(g, g1, g2, 0, prec);
+    _gamma_upper_workaround(g, g1, g2, 0, prec);
 
     /* res = a*b*A*g */
     arb_mul_si(res, g, A, prec);