partially refactor and move power series code to the acb_dirichlet module

acb_dirichlet.h

@@ -27,6 +27,11 @@
 extern "C" {
 #endif
 
+void acb_dirichlet_powsum_term(acb_ptr res, arb_t log_prev, ulong * prev,
+    const acb_t s, ulong k, int integer, int critical_line, slong len, slong prec);
+
+void acb_dirichlet_powsum_sieved(acb_ptr z, const acb_t s, ulong n, slong len, slong prec);
+
 typedef struct
 {
     acb_struct s;

acb_dirichlet/powsum_sieved.c

@@ -0,0 +1,119 @@
+/*
+    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"
+#include "acb_poly.h"
+
+#define POWER(_k) (powers + (((_k)-1)/2) * (len))
+#define DIVISOR(_k) (divisors[((_k)-1)/2])
+
+void
+acb_dirichlet_powsum_sieved(acb_ptr z, const acb_t s, ulong n, slong len, slong prec)
+{
+    slong * divisors;
+    slong powers_alloc;
+    slong i, j, k, ibound, power_of_two, horner_point;
+    ulong kprev;
+    int critical_line, integer;
+
+    acb_ptr powers;
+    acb_ptr t, u, x;
+    acb_ptr p1, p2;
+    arb_t logk, v, w;
+
+    if (n <= 1)
+    {
+        acb_set_ui(z, n);
+        _acb_vec_zero(z + 1, len - 1);
+        return;
+    }
+
+    critical_line = arb_is_exact(acb_realref(s)) &&
+        (arf_cmp_2exp_si(arb_midref(acb_realref(s)), -1) == 0);
+
+    integer = arb_is_zero(acb_imagref(s)) && arb_is_int(acb_realref(s));
+
+    divisors = flint_calloc(n / 2 + 1, sizeof(slong));
+    powers_alloc = (n / 6 + 1) * len;
+    powers = _acb_vec_init(powers_alloc);
+
+    ibound = n_sqrt(n);
+    for (i = 3; i <= ibound; i += 2)
+        if (DIVISOR(i) == 0)
+            for (j = i * i; j <= n; j += 2 * i)
+                DIVISOR(j) = i;
+
+    t = _acb_vec_init(len);
+    u = _acb_vec_init(len);
+    x = _acb_vec_init(len);
+    arb_init(logk);
+    arb_init(v);
+    arb_init(w);
+
+    power_of_two = 1;
+    while (power_of_two * 2 <= n)
+        power_of_two *= 2;
+    horner_point = n / power_of_two;
+
+    _acb_vec_zero(z, len);
+
+    kprev = 1;
+    acb_dirichlet_powsum_term(x, logk, &kprev, s, 2,
+        integer, critical_line, len, prec);
+
+    for (k = 1; k <= n; k += 2)
+    {
+        /* t = k^(-s) */
+        if (DIVISOR(k) == 0)
+        {
+            acb_dirichlet_powsum_term(t, logk, &kprev, s, k,
+                integer, critical_line, len, prec);
+        }
+        else
+        {
+            p1 = POWER(DIVISOR(k));
+            p2 = POWER(k / DIVISOR(k));
+
+            if (len == 1)
+                acb_mul(t, p1, p2, prec);
+            else
+                _acb_poly_mullow(t, p1, len, p2, len, len, prec);
+        }
+
+        if (k * 3 <= n)
+            _acb_vec_set(POWER(k), t, len);
+
+        _acb_vec_add(u, u, t, len, prec);
+
+        while (k == horner_point && power_of_two != 1)
+        {
+            _acb_poly_mullow(t, z, len, x, len, len, prec);
+            _acb_vec_add(z, t, u, len, prec);
+
+            power_of_two /= 2;
+            horner_point = n / power_of_two;
+            horner_point -= (horner_point % 2 == 0);
+        }
+    }
+
+    _acb_poly_mullow(t, z, len, x, len, len, prec);
+    _acb_vec_add(z, t, u, len, prec);
+
+    flint_free(divisors);
+    _acb_vec_clear(powers, powers_alloc);
+    _acb_vec_clear(t, len);
+    _acb_vec_clear(u, len);
+    _acb_vec_clear(x, len);
+    arb_clear(logk);
+    arb_clear(v);
+    arb_clear(w);
+}
+

acb_dirichlet/powsum_term.c

@@ -0,0 +1,73 @@
+/*
+    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_dirichlet_powsum_term(acb_ptr res, arb_t log_prev, ulong * prev,
+    const acb_t s, ulong k, int integer, int critical_line, slong len, slong prec)
+{
+    slong i;
+
+    if (integer)
+    {
+        arb_neg(acb_realref(res), acb_realref(s));
+        arb_set_ui(acb_imagref(res), k);
+        arb_pow(acb_realref(res), acb_imagref(res), acb_realref(res), prec);
+        arb_zero(acb_imagref(res));
+
+        if (len != 1)
+        {
+            arb_log_ui_from_prev(log_prev, k, log_prev, *prev, prec);
+            *prev = k;
+        }
+    }
+    else
+    {
+        arb_t w;
+        arb_init(w);
+
+        arb_log_ui_from_prev(log_prev, k, log_prev, *prev, prec);
+        *prev = k;
+        arb_mul(w, log_prev, acb_imagref(s), prec);
+        arb_sin_cos(acb_imagref(res), acb_realref(res), w, prec);
+        arb_neg(acb_imagref(res), acb_imagref(res));
+
+        if (critical_line)
+        {
+            arb_rsqrt_ui(w, k, prec);
+            acb_mul_arb(res, res, w, prec);
+        }
+        else
+        {
+            arb_mul(w, acb_realref(s), log_prev, prec);
+            arb_neg(w, w);
+            arb_exp(w, w, prec);
+            acb_mul_arb(res, res, w, prec);
+        }
+
+        arb_clear(w);
+    }
+
+    if (len > 1)
+    {
+        arb_neg(log_prev, log_prev);
+
+        for (i = 1; i < len; i++)
+        {
+            acb_mul_arb(res + i, res + i - 1, log_prev, prec);
+            acb_div_ui(res + i, res + i, i, prec);
+        }
+
+        arb_neg(log_prev, log_prev);
+    }
+}
+

acb_poly/powsum_one_series_sieved.c

@@ -10,141 +10,11 @@
 */
 
 #include "acb_poly.h"
-
+#include "acb_dirichlet.h"
-#define POWER(_k) (powers + (((_k)-1)/2) * (len))
-#define DIVISOR(_k) (divisors[((_k)-1)/2])
-
-#define COMPUTE_POWER(t, k, kprev) \
-  do { \
-    if (integer) \
-    { \
-        arb_neg(w, acb_realref(s)); \
-        arb_set_ui(v, k); \
-        arb_pow(acb_realref(t), v, w, prec); \
-        arb_zero(acb_imagref(t)); \
-        if (len != 1) \
-        { \
-            arb_log_ui_from_prev(logk, k, logk, kprev, prec); \
-            kprev = k; \
-            arb_neg(logk, logk); \
-        } \
-    } \
-    else \
-    { \
-        arb_log_ui_from_prev(logk, k, logk, kprev, prec); \
-        kprev = k; \
-        arb_neg(logk, logk); \
-        arb_mul(w, logk, acb_imagref(s), prec); \
-        arb_sin_cos(acb_imagref(t), acb_realref(t), w, prec); \
-        if (critical_line) \
-        { \
-            arb_rsqrt_ui(w, k, prec); \
-            acb_mul_arb(t, t, w, prec); \
-        } \
-        else \
-        { \
-            arb_mul(w, acb_realref(s), logk, prec); \
-            arb_exp(w, w, prec); \
-            acb_mul_arb(t, t, w, prec); \
-        } \
-    } \
-    for (i = 1; i < len; i++) \
-    { \
-        acb_mul_arb(t + i, t + i - 1, logk, prec); \
-        acb_div_ui(t + i, t + i, i, prec); \
-    } \
-    arb_neg(logk, logk); \
-  } while (0); \
 
 void
 _acb_poly_powsum_one_series_sieved(acb_ptr z, const acb_t s, slong n, slong len, slong prec)
 {
-    slong * divisors;
+    acb_dirichlet_powsum_sieved(z, s, n, len, prec);
-    slong powers_alloc;
-    slong i, j, k, ibound, kprev, power_of_two, horner_point;
-    int critical_line, integer;
-
-    acb_ptr powers;
-    acb_ptr t, u, x;
-    acb_ptr p1, p2;
-    arb_t logk, v, w;
-
-    critical_line = arb_is_exact(acb_realref(s)) &&
-        (arf_cmp_2exp_si(arb_midref(acb_realref(s)), -1) == 0);
-
-    integer = arb_is_zero(acb_imagref(s)) && arb_is_int(acb_realref(s));
-
-    divisors = flint_calloc(n / 2 + 1, sizeof(slong));
-    powers_alloc = (n / 6 + 1) * len;
-    powers = _acb_vec_init(powers_alloc);
-
-    ibound = n_sqrt(n);
-    for (i = 3; i <= ibound; i += 2)
-        if (DIVISOR(i) == 0)
-            for (j = i * i; j <= n; j += 2 * i)
-                DIVISOR(j) = i;
-
-    t = _acb_vec_init(len);
-    u = _acb_vec_init(len);
-    x = _acb_vec_init(len);
-    arb_init(logk);
-    arb_init(v);
-    arb_init(w);
-
-    power_of_two = 1;
-    while (power_of_two * 2 <= n)
-        power_of_two *= 2;
-    horner_point = n / power_of_two;
-
-    _acb_vec_zero(z, len);
-
-    kprev = 0;
-    COMPUTE_POWER(x, 2, kprev);
-
-    for (k = 1; k <= n; k += 2)
-    {
-        /* t = k^(-s) */
-        if (DIVISOR(k) == 0)
-        {
-            COMPUTE_POWER(t, k, kprev);
-        }
-        else
-        {
-            p1 = POWER(DIVISOR(k));
-            p2 = POWER(k / DIVISOR(k));
-
-            if (len == 1)
-                acb_mul(t, p1, p2, prec);
-            else
-                _acb_poly_mullow(t, p1, len, p2, len, len, prec);
-        }
-
-        if (k * 3 <= n)
-            _acb_vec_set(POWER(k), t, len);
-
-        _acb_vec_add(u, u, t, len, prec);
-
-        while (k == horner_point && power_of_two != 1)
-        {
-            _acb_poly_mullow(t, z, len, x, len, len, prec);
-            _acb_vec_add(z, t, u, len, prec);
-
-            power_of_two /= 2;
-            horner_point = n / power_of_two;
-            horner_point -= (horner_point % 2 == 0);
-        }
-    }
-
-    _acb_poly_mullow(t, z, len, x, len, len, prec);
-    _acb_vec_add(z, t, u, len, prec);
-
-    flint_free(divisors);
-    _acb_vec_clear(powers, powers_alloc);
-    _acb_vec_clear(t, len);
-    _acb_vec_clear(u, len);
-    _acb_vec_clear(x, len);
-    arb_clear(logk);
-    arb_clear(v);
-    arb_clear(w);
 }
 

doc/source/acb_dirichlet.rst

@@ -22,6 +22,35 @@ The code in other modules for computing the Riemann zeta function,
 Hurwitz zeta function and polylogarithm will possibly be migrated to this
 module in the future.
 
+Truncated L-series and power sums
+-------------------------------------------------------------------------------
+
+.. function:: void acb_dirichlet_powsum_term(acb_ptr res, arb_t log_prev, ulong * prev, const acb_t s, ulong k, int integer, int critical_line, slong len, slong prec)
+
+    Sets *res* to `k^{-(s+x)}` as a power series in *x* truncated to length *len*.
+    The flags *integer* and *critical_line* respectively specify optimizing
+    for *s* being an integer or having real part 1/2.
+
+    On input *log_prev* should contain the natural logarithm of the integer
+    at *prev*. If *prev* is close to *k*, this can be used to speed up
+    computations. If `\log(k)` is computed internally by this function, then
+    *log_prev* is overwritten by this value, and the integer at *prev* is
+    overwritten by *k*, allowing *log_prev* to be recycled for the next
+    term when evaluating a power sum.
+
+.. function:: void acb_dirichlet_powsum_sieved(acb_ptr res, const acb_t s, ulong n, slong len, slong prec)
+
+    Sets *res* to `\sum_{k=1}^n k^{-(s+x)}`
+    as a power series in *x* truncated to length *len*.
+    This function stores a table of powers that have already been calculated,
+    computing `(ij)^r` as `i^r j^r` whenever `k = ij` is
+    composite. As a further optimization, it groups all even `k` and
+    evaluates the sum as a polynomial in `2^{-(s+x)}`.
+    This scheme requires about `n / \log n` powers, `n / 2` multiplications,
+    and temporary storage of `n / 6` power series. Due to the extra
+    power series multiplications, it is only faster than the naive
+    algorithm when *len* is small.
+
 Hurwitz zeta function
 -------------------------------------------------------------------------------