From 06df4d6dcb3c266500184d5b6cc2e52108242bcc Mon Sep 17 00:00:00 2001
From: Fredrik Johansson <fredrik.johansson@gmail.com>
Date: Tue, 25 Oct 2016 21:57:37 +0200
Subject: [PATCH] wrapper for zeta algorithms

---
 acb/zeta.c                   | 52 +---------------------
 acb_dirichlet.h              |  1 +
 acb_dirichlet/zeta.c         | 83 ++++++++++++++++++++++++++++++++++++
 arf.h                        | 12 ++----
 arf/cmpabs.c                 | 15 +++++++
 doc/source/acb_dirichlet.rst |  4 ++
 doc/source/arf.rst           |  2 +
 7 files changed, 111 insertions(+), 58 deletions(-)
 create mode 100644 acb_dirichlet/zeta.c

diff --git a/acb/zeta.c b/acb/zeta.c
index 68f18d88..7be9c33b 100644
--- a/acb/zeta.c
+++ b/acb/zeta.c
@@ -11,6 +11,7 @@
 
 #include "acb.h"
 #include "acb_poly.h"
+#include "acb_dirichlet.h"
 
 void
 acb_zeta_si(acb_t z, slong s, slong prec)
@@ -46,55 +47,6 @@ acb_hurwitz_zeta(acb_t z, const acb_t s, const acb_t a, slong prec)
 void
 acb_zeta(acb_t z, const acb_t s, slong prec)
 {
-    acb_t a;
-    acb_init(a);
-    acb_one(a);
-
-    if (acb_is_int(s) &&
-        arf_cmpabs_2exp_si(arb_midref(acb_realref(s)), FLINT_BITS - 1) < 0)
-    {
-        acb_zeta_si(z, arf_get_si(arb_midref(acb_realref(s)), ARF_RND_DOWN), prec);
-        return;
-    }
-
-    if (arf_sgn(arb_midref(acb_realref(s))) < 0)
-    {
-        acb_t t, u, v;
-        slong wp = prec + 6;
-
-        acb_init(t);
-        acb_init(u);
-        acb_init(v);
-
-        acb_sub_ui(t, s, 1, wp);
-
-        /* 2 * (2pi)^(s-1) */
-        arb_const_pi(acb_realref(u), wp);
-        acb_mul_2exp_si(u, u, 1);
-        acb_pow(u, u, t, wp);
-        acb_mul_2exp_si(u, u, 1);
-
-        /* sin(pi*s/2) */
-        acb_mul_2exp_si(v, s, -1);
-        acb_sin_pi(v, v, wp);
-        acb_mul(u, u, v, wp);
-
-        /* gamma(1-s) zeta(1-s) */
-        acb_neg(t, t);
-        acb_gamma(v, t, wp);
-        acb_mul(u, u, v, wp);
-        acb_hurwitz_zeta(v, t, a, wp);
-        acb_mul(z, u, v, prec);
-
-        acb_clear(t);
-        acb_clear(u);
-        acb_clear(v);
-    }
-    else
-    {
-        acb_hurwitz_zeta(z, s, a, prec);
-    }
-
-    acb_clear(a);
+    acb_dirichlet_zeta(z, s, prec);
 }
 
diff --git a/acb_dirichlet.h b/acb_dirichlet.h
index 8b106e60..f0f84274 100644
--- a/acb_dirichlet.h
+++ b/acb_dirichlet.h
@@ -38,6 +38,7 @@ void acb_dirichlet_zeta_rs_d_coeffs(arb_ptr d, const arb_t sigma, slong k, slong
 void acb_dirichlet_zeta_rs_bound(mag_t err, const acb_t s, slong K);
 void acb_dirichlet_zeta_rs_r(acb_t res, const acb_t s, slong K, slong prec);
 void acb_dirichlet_zeta_rs(acb_t res, const acb_t s, slong K, slong prec);
+void acb_dirichlet_zeta(acb_t res, const acb_t s, slong prec);
 
 typedef struct
 {
diff --git a/acb_dirichlet/zeta.c b/acb_dirichlet/zeta.c
new file mode 100644
index 00000000..8452df69
--- /dev/null
+++ b/acb_dirichlet/zeta.c
@@ -0,0 +1,83 @@
+/*
+    Copyright (C) 2016 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_dirichlet.h"
+
+void acb_zeta_si(acb_t z, slong s, slong prec);
+
+void
+acb_dirichlet_zeta(acb_t res, const acb_t s, slong prec)
+{
+    acb_t a;
+    double cutoff;
+
+    if (acb_is_int(s) &&
+        arf_cmpabs_2exp_si(arb_midref(acb_realref(s)), FLINT_BITS - 1) < 0)
+    {
+        acb_zeta_si(res, arf_get_si(arb_midref(acb_realref(s)), ARF_RND_DOWN), prec);
+        return;
+    }
+
+    cutoff = 24.0 * prec * sqrt(prec);
+
+    /* todo: support non-exact s properly in RS */
+    if (acb_is_exact(s) &&
+        arf_cmpabs_d(arb_midref(acb_imagref(s)), cutoff) >= 0 &&
+        arf_cmpabs_d(arb_midref(acb_realref(s)), 10 + prec * 0.1) <= 0)
+    {
+        acb_dirichlet_zeta_rs(res, s, 0, prec);
+        return;
+    }
+
+    acb_init(a);
+    acb_one(a);
+
+    if (arf_sgn(arb_midref(acb_realref(s))) < 0)
+    {
+        acb_t t, u, v;
+        slong wp = prec + 6;
+
+        acb_init(t);
+        acb_init(u);
+        acb_init(v);
+
+        acb_sub_ui(t, s, 1, wp);
+
+        /* 2 * (2pi)^(s-1) */
+        arb_const_pi(acb_realref(u), wp);
+        acb_mul_2exp_si(u, u, 1);
+        acb_pow(u, u, t, wp);
+        acb_mul_2exp_si(u, u, 1);
+
+        /* sin(pi*s/2) */
+        acb_mul_2exp_si(v, s, -1);
+        acb_sin_pi(v, v, wp);
+        acb_mul(u, u, v, wp);
+
+        /* gamma(1-s) zeta(1-s) */
+        acb_neg(t, t);
+        acb_gamma(v, t, wp);
+        acb_mul(u, u, v, wp);
+        acb_hurwitz_zeta(v, t, a, wp);
+        acb_mul(res, u, v, prec);
+
+        acb_clear(t);
+        acb_clear(u);
+        acb_clear(v);
+    }
+    else
+    {
+        acb_hurwitz_zeta(res, s, a, prec);
+    }
+
+    acb_clear(a);
+}
+
diff --git a/arf.h b/arf.h
index b4780b82..1253622e 100644
--- a/arf.h
+++ b/arf.h
@@ -343,6 +343,10 @@ int arf_cmp(const arf_t x, const arf_t y);
 
 int arf_cmpabs(const arf_t x, const arf_t y);
 
+int arf_cmpabs_ui(const arf_t x, ulong y);
+
+int arf_cmpabs_d(const arf_t x, double y);
+
 int arf_cmp_si(const arf_t x, slong y);
 
 int arf_cmp_ui(const arf_t x, ulong y);
@@ -466,14 +470,6 @@ arf_set_si(arf_t x, slong v)
         ARF_NEG(x);
 }
 
-ARF_INLINE int
-arf_cmpabs_ui(const arf_t x, ulong y)
-{
-    arf_t t;
-    arf_init_set_ui(t, y);  /* no need to free */
-    return arf_cmpabs(x, t);
-}
-
 ARF_INLINE void
 arf_init_set_shallow(arf_t z, const arf_t x)
 {
diff --git a/arf/cmpabs.c b/arf/cmpabs.c
index 7d4540a2..81019415 100644
--- a/arf/cmpabs.c
+++ b/arf/cmpabs.c
@@ -53,3 +53,18 @@ arf_cmpabs(const arf_t x, const arf_t y)
     return 0;
 }
 
+int arf_cmpabs_ui(const arf_t x, ulong y)
+{
+    arf_t t;
+    arf_init_set_ui(t, y); /* no need to free */
+    return arf_cmpabs(x, t);
+}
+
+int arf_cmpabs_d(const arf_t x, double y)
+{
+    arf_t t;
+    arf_init(t); /* no need to free */
+    arf_set_d(t, y);
+    return arf_cmpabs(x, t);
+}
+
diff --git a/doc/source/acb_dirichlet.rst b/doc/source/acb_dirichlet.rst
index 7a067221..4a7df33d 100644
--- a/doc/source/acb_dirichlet.rst
+++ b/doc/source/acb_dirichlet.rst
@@ -121,6 +121,10 @@ the evaluation (automatic reduction to the exact case is not yet implemented).
     otherwise chooses the number of terms automatically based on *s* and the
     precision.
 
+.. function:: void acb_dirichlet_zeta(acb_t res, const acb_t s, slong prec)
+
+    Computes `\zeta(s)` using an automatic choice of algorithm.
+
 Hurwitz zeta function
 -------------------------------------------------------------------------------
 
diff --git a/doc/source/arf.rst b/doc/source/arf.rst
index a34ab6f1..9529ebee 100644
--- a/doc/source/arf.rst
+++ b/doc/source/arf.rst
@@ -344,6 +344,8 @@ Comparisons and bounds
 
 .. function:: int arf_cmpabs_ui(const arf_t x, ulong y)
 
+.. function:: int arf_cmpabs_d(const arf_t x, ulong y)
+
 .. function:: int arf_cmpabs_mag(const arf_t x, const mag_t y)
 
     Compares the absolute values of *x* and *y*.