simplified init call for hurwitz_precomp

acb_dirichlet.h

@@ -59,6 +59,7 @@ acb_dirichlet_hurwitz_precomp_struct;
 typedef acb_dirichlet_hurwitz_precomp_struct acb_dirichlet_hurwitz_precomp_t[1];
 
 void acb_dirichlet_hurwitz_precomp_init(acb_dirichlet_hurwitz_precomp_t pre, const acb_t s, int deflate, slong A, slong K, slong N, slong prec);
+void acb_dirichlet_hurwitz_precomp_init_num(acb_dirichlet_hurwitz_precomp_t pre, const acb_t s, int deflate, double num_eval, slong prec);
 void acb_dirichlet_hurwitz_precomp_clear(acb_dirichlet_hurwitz_precomp_t pre);
 void acb_dirichlet_hurwitz_precomp_bound(mag_t res, const acb_t s, slong A, slong K, slong N);
 void acb_dirichlet_hurwitz_precomp_eval(acb_t res, const acb_dirichlet_hurwitz_precomp_t pre, ulong p, ulong q, slong prec);

acb_dirichlet/hurwitz_precomp_clear.c

@@ -15,7 +15,11 @@ void
 acb_dirichlet_hurwitz_precomp_clear(acb_dirichlet_hurwitz_precomp_t pre)
 {
     acb_clear(&pre->s);
+
+    if (pre->A != 0)
+    {
         mag_clear(&pre->err);
         _acb_vec_clear(pre->coeffs, pre->N * pre->K);
+    }
 }
 

acb_dirichlet/hurwitz_precomp_eval.c

@@ -20,7 +20,25 @@ acb_dirichlet_hurwitz_precomp_eval(acb_t res,
     acb_t a, t;
 
     if (p > q)
+    {
+        flint_printf("hurwitz_precomp_eval: require p <= n\n");
         flint_abort();
+    }
+
+    if (pre->A == 0)
+    {
+        acb_init(a);
+        acb_set_ui(a, p);
+        acb_div_ui(a, a, q, prec);
+
+        if (pre->deflate == 0)
+            acb_hurwitz_zeta(res, &pre->s, a, prec);
+        else
+            _acb_poly_zeta_cpx_series(res, &pre->s, a, 1, 1, prec);
+
+        acb_clear(a);
+        return;
+    }
 
     acb_init(a);
     acb_init(t);

acb_dirichlet/hurwitz_precomp_init.c

@@ -18,18 +18,26 @@ acb_dirichlet_hurwitz_precomp_init(acb_dirichlet_hurwitz_precomp_t pre,
 {
     slong i, k;
 
-    if (A < 1 || K < 1 || N < 1)
-        flint_abort();
-
     pre->deflate = deflate;
     pre->A = A;
     pre->K = K;
     pre->N = N;
+
+    acb_init(&pre->s);
+    acb_set(&pre->s, s);
+
+    if (A == 0)
+        return;
+
+    if (A < 1 || K < 1 || N < 1)
+    {
+        flint_printf("hurwitz_precomp_init: require A, K, N >= 1 (unless A == 0)\n");
+        flint_abort();
+    }
+
     pre->coeffs = _acb_vec_init(N * K);
 
     mag_init(&pre->err);
-    acb_init(&pre->s);
-    acb_set(&pre->s, s);
 
     acb_dirichlet_hurwitz_precomp_bound(&pre->err, s, A, K, N);
 
@@ -79,3 +87,12 @@ acb_dirichlet_hurwitz_precomp_init(acb_dirichlet_hurwitz_precomp_t pre,
     }
 }
 
+void
+acb_dirichlet_hurwitz_precomp_init_num(acb_dirichlet_hurwitz_precomp_t pre,
+        const acb_t s, int deflate, double num_eval, slong prec)
+{
+    ulong A, K, N;
+    acb_dirichlet_hurwitz_precomp_choose_param(&A, &K, &N, s, num_eval, prec);
+    acb_dirichlet_hurwitz_precomp_init(pre, s, deflate, A, K, N, prec);
+}
+

acb_dirichlet/l.c

@@ -23,24 +23,12 @@ acb_dirichlet_l_general(acb_t res, const acb_t s,
     }
     else
     {
-        ulong A, K, N;
-        slong wp;
-
-        wp = prec + n_clog(G->phi_q, 2);
-        acb_dirichlet_hurwitz_precomp_choose_param(&A, &K, &N, s, G->phi_q, wp);
-
-        if (A == 0)
-        {
-            acb_dirichlet_l_hurwitz(res, s, NULL, G, chi, prec);
-        }
-        else
-        {
+        slong wp = prec + n_clog(G->phi_q, 2);
         acb_dirichlet_hurwitz_precomp_t pre;
-            acb_dirichlet_hurwitz_precomp_init(pre, s, acb_is_one(s), A, K, N, wp);
+        acb_dirichlet_hurwitz_precomp_init_num(pre, s, acb_is_one(s), G->phi_q, wp);
         acb_dirichlet_l_hurwitz(res, s, pre, G, chi, prec);
         acb_dirichlet_hurwitz_precomp_clear(pre);
     }
-    }
 }
 
 void

acb_dirichlet/test/t-hurwitz_precomp.c

@@ -43,6 +43,7 @@ int main()
         acb_init(z2);
 
         acb_randtest(s, state, 1 + n_randint(state, 200), 2);
+
         acb_dirichlet_hurwitz_precomp_init(pre, s, deflate, A, K, N, prec1);
 
         for (i = 0; i < 10; i++)

acb_dirichlet/test/t-l_vec_hurwitz.c

@@ -45,13 +45,18 @@ int main()
 
         v = _acb_vec_init(G->phi_q);
 
+        if (n_randint(state, 2))
+        {
             acb_dirichlet_hurwitz_precomp_choose_param(&A, &K, &N, s, G->phi_q, prec);
-
-        if (A != 0)
             acb_dirichlet_hurwitz_precomp_init(pre, s, acb_is_one(s), A, K, N, prec);
+        }
+        else
+        {
+            acb_dirichlet_hurwitz_precomp_init_num(pre, s, acb_is_one(s), G->phi_q, prec);
+        }
 
         /* all at once */
-        if (A != 0 && n_randint(state, 2))
+        if (n_randint(state, 2))
             acb_dirichlet_l_vec_hurwitz(v, s, pre, G, prec);
         else
             acb_dirichlet_l_vec_hurwitz(v, s, NULL, G, prec);
@@ -103,7 +108,6 @@ int main()
         dirichlet_char_clear(chi);
         dirichlet_group_clear(G);
 
-        if (A != 0)
         acb_dirichlet_hurwitz_precomp_clear(pre);
     }
 

doc/source/acb_dirichlet.rst

@@ -199,7 +199,11 @@ Hurwitz zeta function precomputation
     requires *2K* simple arithmetic operations (polynomial evaluation) plus
     *A* powers. As *K* grows, the error is at most `O(1/(2AN)^K)`.
 
-    We require that *A*, *K* and *N* are all positive. Moreover, for a finite
+    This function can be called with *A* set to zero, in which case
+    no Taylor series precomputation is performed. This means that evaluation
+    will be identical to calling :func:`acb_dirichlet_hurwitz` directly.
+
+    Otherwise, we require that *A*, *K* and *N* are all positive. For a finite
     error bound, we require `K+\operatorname{re}(s) > 1`.
     To avoid an initial "bump" that steals precision
     and slows convergence, *AN* should be at least roughly as large as `|s|`,
@@ -208,6 +212,12 @@ Hurwitz zeta function precomputation
     If *deflate* is set, the deflated Hurwitz zeta function is used,
     removing the pole at `s = 1`.
 
+.. function:: void acb_dirichlet_hurwitz_precomp_init_num(acb_dirichlet_hurwitz_precomp_t pre, const acb_t s, int deflate, double num_eval, slong prec)
+
+    Initializes *pre*, choosing the parameters *A*, *K*, and *N*
+    automatically to minimize the cost of *num_eval* evaluations of the
+    Hurwitz zeta function at argument *s* to precision *prec*.
+
 .. function:: void acb_dirichlet_hurwitz_precomp_clear(acb_dirichlet_hurwitz_precomp_t pre)
 
     Clears the precomputed data.

examples/lcentral.c

@@ -18,7 +18,6 @@ int main(int argc, char *argv[])
     slong prec = 100, digits = 30;
     ulong qmin, qmax, q;
     acb_t s;
-    ulong A, K, N;
     acb_dirichlet_hurwitz_precomp_t pre;
 
     if (argc < 3)
@@ -53,10 +52,8 @@ int main(int argc, char *argv[])
     acb_one(s);
     acb_div_si(s, s, 2, prec);
 
-    acb_dirichlet_hurwitz_precomp_choose_param(&A, &K, &N, s,
+    acb_dirichlet_hurwitz_precomp_init_num(pre, s, 0,
         (qmax - qmin + 1) * 0.5 * qmax, prec);
-    if (A != 0)
-        acb_dirichlet_hurwitz_precomp_init(pre, s, 0, A, K, N, prec);
 
     for (q = qmin; q <= qmax; q++)
     {
@@ -73,7 +70,7 @@ int main(int argc, char *argv[])
 
         z =  _acb_vec_init(G->phi_q);
 
-        acb_dirichlet_l_vec_hurwitz(z, s, (A == 0) ? NULL : pre, G, prec);
+        acb_dirichlet_l_vec_hurwitz(z, s, pre, G, prec);
 
         if (out || check)
         {
@@ -107,7 +104,6 @@ int main(int argc, char *argv[])
         dirichlet_group_clear(G);
     }
 
-    if (A != 0)
     acb_dirichlet_hurwitz_precomp_clear(pre);
 
     acb_clear(s);