implement Euler product for L-functions (acb_dirichlet_l_euler_product)

acb_dirichlet.h

@@ -87,6 +87,8 @@ void acb_dirichlet_jacobi_sum_gauss(acb_t res, const dirichlet_group_t G, const
 void acb_dirichlet_jacobi_sum(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec);
 void acb_dirichlet_jacobi_sum_ui(acb_t res, const dirichlet_group_t G, ulong a, ulong b, slong prec);
 
+void acb_dirichlet_l_euler_product(acb_t res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
+
 void acb_dirichlet_l_hurwitz(acb_t res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
 void acb_dirichlet_l(acb_t res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
 

acb_dirichlet/l.c

@@ -16,6 +16,14 @@ void
 acb_dirichlet_l(acb_t res, const acb_t s,
     const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
 {
-    acb_dirichlet_l_hurwitz(res, s, G, chi, prec);
+    /* this cutoff is probably too conservative when q is large */
+    if (arf_cmp_d(arb_midref(acb_realref(s)), 8 + 0.5 * prec / log(prec)) >= 0)
+    {
+        acb_dirichlet_l_euler_product(res, s, G, chi, prec);
+    }
+    else
+    {
+        acb_dirichlet_l_hurwitz(res, s, G, chi, prec);
+    }
 }
 

acb_dirichlet/l_euler_product.c

@@ -0,0 +1,128 @@
+/*
+    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"
+
+#define ONE_OVER_LOG2 1.4426950408889634
+
+void
+acb_dirichlet_l_euler_product(acb_t res, const acb_t s,
+    const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
+{
+    arf_t left;
+    slong wp, powprec, left_s;
+    ulong p, p_limit;
+    double p_needed_approx, powmag, logp, errmag;
+    int is_real;
+    acb_t t, u, v, c;
+    mag_t err;
+
+    if (!acb_is_finite(s))
+    {
+        acb_indeterminate(res);
+        return;
+    }
+
+    arf_init(left);
+    arf_set_mag(left, arb_radref(acb_realref(s)));
+    arf_sub(left, arb_midref(acb_realref(s)), left, 2 * MAG_BITS, ARF_RND_FLOOR);
+
+    /* Require re(s) >= 2. */
+    if (arf_cmp_si(left, 2) < 0)
+    {
+        acb_indeterminate(res);
+        arf_clear(left);
+        return;
+    }
+
+    is_real = acb_is_real(s) && dirichlet_char_is_real(G, chi);
+
+    /* L(s) ~= 1. */
+    if (arf_cmp_si(left, prec) > 0)
+    {
+        acb_one(res);
+        mag_hurwitz_zeta_uiui(arb_radref(acb_realref(res)), prec, 2);
+        if (!is_real)
+            mag_set(arb_radref(acb_imagref(res)), arb_radref(acb_realref(res)));
+        acb_inv(res, res, prec);
+        arf_clear(left);
+        return;
+    }
+
+    left_s = arf_get_si(left, ARF_RND_FLOOR);
+
+    /* Heuristic. */
+    wp = prec + FLINT_BIT_COUNT(prec) + (prec / left_s) + 4;
+
+    /* Don't work too hard if a small s was passed as input. */
+    p_limit = 100 + prec * sqrt(prec);
+
+    /* Truncating at p ~= 2^(prec/s) gives an error of 2^-prec */
+    if (((double) prec) / left_s > 50.0)
+        p_needed_approx = pow(2.0, 50.0);
+    else
+        p_needed_approx = pow(2.0, ((double) prec) / left_s);
+
+    p_needed_approx = FLINT_MIN(p_limit, p_needed_approx);
+    /* todo: use this to determine whether to precompute chi */
+
+    acb_init(t);
+    acb_init(u);
+    acb_init(v);
+    acb_init(c);
+
+    acb_one(v);
+
+    for (p = 2; p < p_limit; p = n_nextprime(p, 1))
+    {
+        /* p^s */
+        logp = log(p);
+        powmag = left_s * logp * ONE_OVER_LOG2;
+
+        /* zeta(s,p) ~= 1/p^s + 1/((s-1) p^(s-1)) */
+        errmag = (log(left_s - 1.0) + (left_s - 1.0) * logp) * ONE_OVER_LOG2;
+        errmag = FLINT_MIN(powmag, errmag);
+
+        if (errmag > prec + 2)
+            break;
+
+        powprec = FLINT_MAX(wp - powmag, 8);
+
+        acb_dirichlet_chi(c, G, chi, p, powprec);
+
+        if (!acb_is_zero(c))
+        {
+            acb_set_ui(t, p);
+            acb_pow(t, t, s, powprec);
+            acb_set_round(u, v, powprec);
+            acb_div(t, u, t, powprec);
+            acb_mul(t, t, c, powprec);
+            acb_sub(v, v, t, wp);
+        }
+    }
+
+    mag_init(err);
+    mag_hurwitz_zeta_uiui(err, left_s, p);
+    if (is_real)
+        arb_add_error_mag(acb_realref(v), err);
+    else
+        acb_add_error_mag(v, err);
+    mag_clear(err);
+
+    acb_inv(res, v, prec);
+
+    acb_clear(t);
+    acb_clear(u);
+    acb_clear(v);
+    acb_clear(c);
+    arf_clear(left);
+}
+

acb_dirichlet/test/t-l_euler_product.c

@@ -0,0 +1,85 @@
+/*
+    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"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("l_euler_product....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 500 * arb_test_multiplier(); iter++)
+    {
+        acb_t s, t, u;
+        dirichlet_group_t G;
+        dirichlet_char_t chi;
+        ulong q, k;
+        slong prec;
+
+        acb_init(s);
+        acb_init(t);
+        acb_init(u);
+
+        q = 1 + n_randint(state, 50);
+        prec = 2 + n_randint(state, 500);
+        k = n_randint(state, n_euler_phi(q));
+
+        dirichlet_group_init(G, q);
+        dirichlet_char_init(chi, G);
+        dirichlet_char_index(chi, G, k);
+
+        if (n_randint(state, 2))
+        {
+            acb_set_ui(s, 2 + n_randint(state, prec + 50));
+        }
+        else
+        {
+            acb_randtest(s, state, 2 + n_randint(state, 200), 2);
+            arb_abs(acb_realref(s), acb_realref(s));
+            if (n_randint(state, 2))
+                acb_add_ui(s, s, n_randtest(state), prec);
+        }
+
+        if (n_randint(state, 2))
+            acb_dirichlet_l_hurwitz(t, s, G, chi, prec);
+        else
+            acb_dirichlet_l_euler_product(t, s, G, chi, prec * 1.5);
+
+        acb_dirichlet_l_euler_product(u, s, G, chi, prec);
+
+        if (!acb_overlaps(t, u))
+        {
+            flint_printf("FAIL: overlap\n\n");
+            flint_printf("iter = %ld  q = %lu  k = %lu  prec = %ld\n\n", iter, q, k, prec);
+            flint_printf("s = "); acb_printn(s, 100, 0); flint_printf("\n\n");
+            flint_printf("t = "); acb_printn(t, 100, 0); flint_printf("\n\n");
+            flint_printf("u = "); acb_printn(u, 100, 0); flint_printf("\n\n");
+            abort();
+        }
+
+        dirichlet_char_clear(chi);
+        dirichlet_group_clear(G);
+        acb_clear(s);
+        acb_clear(t);
+        acb_clear(u);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

doc/source/acb_dirichlet.rst

@@ -206,8 +206,13 @@ L-functions
 
     This formula is slow for large *q*.
 
+.. function:: void acb_dirichlet_l_euler_product(acb_t res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
+
+    Computes `L(s,\chi)` directly using the Euler product. This is
+    efficient if *s* has large positive real part. As implemented, this
+    function only gives a finite result if `\operatorname{re}(s) \ge 2`.
+
 .. function:: void acb_dirichlet_l(acb_t res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
 
     Computes `L(s,\chi)` using a default choice of algorithm.
-    Currently, only the Hurwitz zeta formula is used.