From 66741901ad5a76fe598d11eb39a46f021e8e0c01 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Wed, 16 Nov 2016 22:49:45 +0100 Subject: [PATCH] implement Hardy Z-function for Dirichlet L-functions --- acb_dirichlet.h | 9 +++ acb_dirichlet/hardy_theta.c | 126 +++++++++++++++++++++++++++++++++ acb_dirichlet/hardy_z.c | 62 ++++++++++++++++ acb_dirichlet/test/t-hardy_z.c | 122 +++++++++++++++++++++++++++++++ doc/source/acb_dirichlet.rst | 25 ++++++- doc/source/dirichlet.rst | 5 +- examples/real_roots.c | 89 +++++++++++++++++------ 7 files changed, 415 insertions(+), 23 deletions(-) create mode 100644 acb_dirichlet/hardy_theta.c create mode 100644 acb_dirichlet/hardy_z.c create mode 100644 acb_dirichlet/test/t-hardy_z.c diff --git a/acb_dirichlet.h b/acb_dirichlet.h index e7eb7c67..2dddc0fb 100644 --- a/acb_dirichlet.h +++ b/acb_dirichlet.h @@ -126,6 +126,15 @@ void acb_dirichlet_l(acb_t res, const acb_t s, const dirichlet_group_t G, const void acb_dirichlet_l_jet(acb_ptr res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, int deflate, slong len, slong prec); +void acb_dirichlet_hardy_theta(acb_ptr res, const acb_t t, + const dirichlet_group_t G, const dirichlet_char_t chi, + slong len, slong prec); + +void acb_dirichlet_hardy_z(acb_t res, const acb_t t, + const dirichlet_group_t G, const dirichlet_char_t chi, + slong len, slong prec); + + /* utils */ ACB_DIRICHLET_INLINE void diff --git a/acb_dirichlet/hardy_theta.c b/acb_dirichlet/hardy_theta.c new file mode 100644 index 00000000..72f1651b --- /dev/null +++ b/acb_dirichlet/hardy_theta.c @@ -0,0 +1,126 @@ +/* + 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 . +*/ + +#include "acb_dirichlet.h" +#include "acb_poly.h" + +void +acb_dirichlet_hardy_theta(acb_ptr res, const acb_t t, + const dirichlet_group_t G, const dirichlet_char_t chi, + slong len, slong prec) +{ + acb_struct y[2]; + arb_t c; + slong k; + int parity; + ulong q; + + if (len <= 0) + return; + + if (t == res) + { + acb_init(y); + acb_set(y, t); + acb_dirichlet_hardy_theta(res, y, G, chi, len, prec); + acb_clear(y); + return; + } + + if (G == NULL) + { + parity = 0; + q = 1; + } + else + { + parity = dirichlet_parity_char(G, chi); + q = G->q; + + if (q != dirichlet_conductor_char(G, chi)) + { + flint_printf("hardy theta: need primitive character\n"); + abort(); + } + } + + if (!acb_is_finite(t)) + { + _acb_vec_indeterminate(res, len); + return; + } + + acb_init(y + 0); + acb_init(y + 1); + arb_init(c); + + /* res = log gamma((s+parity)/2), s = 0.5+i(t+x) */ + acb_mul_onei(y, t); + arb_set_d(c, 0.5 + parity); + arb_add(acb_realref(y), acb_realref(y), c, prec); + acb_mul_2exp_si(y, y, -1); + acb_onei(y + 1); + acb_mul_2exp_si(y + 1, y + 1, -1); + _acb_poly_lgamma_series(res, y, FLINT_MIN(len, 2), len, prec); + + if (acb_is_real(t)) + { + for (k = 0; k < len; k++) + { + arb_set(acb_realref(res + k), acb_imagref(res + k)); + arb_zero(acb_imagref(res + k)); + } + } + else + { + acb_ptr v = _acb_vec_init(len); + + /* v = log gamma(((1-s)+parity)/2), s = 0.5+i(t+x) */ + acb_div_onei(y, t); + arb_set_d(c, 0.5 + parity); + arb_add(acb_realref(y), acb_realref(y), c, prec); + acb_mul_2exp_si(y, y, -1); + acb_neg(y + 1, y + 1); + _acb_poly_lgamma_series(v, y, FLINT_MIN(len, 2), len, prec); + + _acb_vec_sub(res, res, v, len, prec); + for (k = 0; k < len; k++) + { + acb_div_onei(res + k, res + k); + acb_mul_2exp_si(res + k, res + k, -1); + } + + _acb_vec_clear(v, len); + } + + /* (t+x) [-(1/2) log(pi/q)] */ + arb_const_pi(c, prec); + arb_div_ui(c, c, q, prec); + arb_log(c, c, prec); + arb_mul_2exp_si(c, c, -1); + acb_submul_arb(res, t, c, prec); + if (len > 1) + acb_sub_arb(res + 1, res + 1, c, prec); + + /* i log(eps) / 2 */ + if (q > 1) + { + acb_dirichlet_root_number(y, G, chi, prec); + acb_arg(c, y, prec); + arb_mul_2exp_si(c, c, -1); + arb_sub(acb_realref(res), acb_realref(res), c, prec); + } + + acb_clear(y + 0); + acb_clear(y + 1); + arb_clear(c); +} + diff --git a/acb_dirichlet/hardy_z.c b/acb_dirichlet/hardy_z.c new file mode 100644 index 00000000..cb0bc5d9 --- /dev/null +++ b/acb_dirichlet/hardy_z.c @@ -0,0 +1,62 @@ +/* + 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 . +*/ + +#include "acb_dirichlet.h" +#include "acb_poly.h" + +void +acb_dirichlet_hardy_z(acb_t res, const acb_t t, + const dirichlet_group_t G, const dirichlet_char_t chi, + slong len, slong prec) +{ + acb_ptr v, w; + slong k; + int is_real; + + if (len <= 0) + return; + + v = _acb_vec_init(len); + w = _acb_vec_init(len); + + is_real = acb_is_real(t); + + /* v = exp(i theta(t+x)) */ + acb_dirichlet_hardy_theta(v, t, G, chi, len, prec); + _acb_vec_scalar_mul_onei(v, v, len); + _acb_poly_exp_series(v, v, len, len, prec); + + /* w = L(1/2 + i (t+x)) */ + acb_one(w); + acb_mul_2exp_si(w, w, -1); + arb_sub(acb_realref(w), acb_realref(w), acb_imagref(t), prec); + arb_set(acb_imagref(w), acb_realref(t)); + acb_dirichlet_l_jet(w, w, G, chi, 0, len, prec); + for (k = 0; k < len; k++) + { + if (k % 4 == 1) + acb_mul_onei(w + k, w + k); + else if (k % 4 == 2) + acb_neg(w + k, w + k); + else if (k % 4 == 3) + acb_div_onei(w + k, w + k); + } + + _acb_poly_mullow(res, v, len, w, len, len, prec); + + if (is_real) + for (k = 0; k < len; k++) + arb_zero(acb_imagref(res + k)); + + _acb_vec_clear(v, len); + _acb_vec_clear(w, len); +} + diff --git a/acb_dirichlet/test/t-hardy_z.c b/acb_dirichlet/test/t-hardy_z.c new file mode 100644 index 00000000..626afbcc --- /dev/null +++ b/acb_dirichlet/test/t-hardy_z.c @@ -0,0 +1,122 @@ +/* + 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 . +*/ + +#include "acb_dirichlet.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("hardy_z...."); + fflush(stdout); + + flint_randinit(state); + + /* test self-consistency */ + for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) + { + acb_t s, s2; + dirichlet_group_t G; + dirichlet_char_t chi; + acb_ptr vec1, vec2; + slong len1, len2; + slong prec1, prec2; + ulong q, k; + slong i; + + len1 = n_randint(state, 6); + len2 = n_randint(state, 6); + prec1 = 2 + n_randint(state, 100); + prec2 = 2 + n_randint(state, 100); + + do { + q = 1 + n_randint(state, 30); + } while (q % 4 == 2); + + dirichlet_group_init(G, q); + dirichlet_char_init(chi, G); + + do { + k = n_randint(state, n_euler_phi(q)); + dirichlet_char_index(chi, G, k); + } while (dirichlet_conductor_char(G, chi) != q); + + acb_init(s); + acb_init(s2); + vec1 = _acb_vec_init(len1); + vec2 = _acb_vec_init(len2); + + acb_randtest(s, state, 2 + n_randint(state, 200), 2); + acb_randtest(s2, state, 2 + n_randint(state, 200), 2); + acb_sub(s2, s2, s2, 200); + acb_add(s2, s, s2, 200); + + acb_dirichlet_hardy_z(vec1, s, G, chi, len1, prec1); + acb_dirichlet_hardy_z(vec2, s2, G, chi, len2, prec2); + + for (i = 0; i < FLINT_MIN(len1, len2); i++) + { + if (!acb_overlaps(vec1 + i, vec2 + i)) + { + flint_printf("FAIL: overlap\n\n"); + flint_printf("iter = %wd q = %wu k = %wu i = %wd\n\n", iter, q, k, i); + flint_printf("s = "); acb_printn(s, 50, 0); flint_printf("\n\n"); + flint_printf("r1 = "); acb_printn(vec1 + i, 50, 0); flint_printf("\n\n"); + flint_printf("r2 = "); acb_printn(vec2 + i, 50, 0); flint_printf("\n\n"); + abort(); + } + } + + if (arb_contains_zero(acb_imagref(s))) + { + for (i = 0; i < len1; i++) + { + if (!arb_contains_zero(acb_imagref(vec1 + i))) + { + flint_printf("FAIL: real 1\n\n"); + flint_printf("iter = %wd q = %wu k = %wu i = %wd\n\n", iter, q, k, i); + flint_printf("s = "); acb_printn(s, 50, 0); flint_printf("\n\n"); + flint_printf("r1 = "); acb_printn(vec1 + i, 50, 0); flint_printf("\n\n"); + abort(); + } + } + } + + if (arb_contains_zero(acb_imagref(s2))) + { + for (i = 0; i < len2; i++) + { + if (!arb_contains_zero(acb_imagref(vec2 + i))) + { + flint_printf("FAIL: real 1\n\n"); + flint_printf("iter = %wd q = %wu k = %wu i = %wd\n\n", iter, q, k, i); + flint_printf("s = "); acb_printn(s, 50, 0); flint_printf("\n\n"); + flint_printf("r1 = "); acb_printn(vec2 + i, 50, 0); flint_printf("\n\n"); + abort(); + } + } + } + + dirichlet_char_clear(chi); + dirichlet_group_clear(G); + acb_clear(s); + acb_clear(s2); + _acb_vec_clear(vec1, len1); + _acb_vec_clear(vec2, len2); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/doc/source/acb_dirichlet.rst b/doc/source/acb_dirichlet.rst index 397b3063..ed2fb986 100644 --- a/doc/source/acb_dirichlet.rst +++ b/doc/source/acb_dirichlet.rst @@ -89,7 +89,7 @@ The coefficients `C_k(p)` in the asymptotic part of the expansion are expressed in terms of certain auxiliary coefficients `d_j^{(k)}` and `F^{(j)}(p)`. Because of artificial discontinuities, *s* should be exact inside -the evaluation (automatic reduction to the exact case is not yet implemented). +the evaluation. .. function:: void acb_dirichlet_zeta_rs_f_coeffs(acb_ptr f, const arb_t p, slong n, slong prec) @@ -373,3 +373,26 @@ L-functions gives the regular part of the Laurent expansion. When *chi* is non-principal, *deflate* has no effect. +Hardy Z-functions +------------------------------------------------------------------------------- + +.. function:: void acb_dirichlet_hardy_theta(acb_ptr res, const acb_t t, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec) + + Computes the phase function used to construct the Z-function. + We have + + .. math :: + + \theta(t) = -\frac{t}{2} \log(\pi/q) - \frac{\epsilon}{2} + + \frac{\log \Gamma((s+\delta)/2) - \log \Gamma((1-s+\delta)/2)}{2i} + + where `s = 1/2+it`, `\delta` is the parity of *chi*, and `\epsilon` + is the root number as computed by :func:`acb_dirichlet_root_number`. + The first *len* terms in the Taylor expansion are written to the output. + +.. function:: void acb_dirichlet_hardy_z(acb_t res, const acb_t t, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec) + + Computes the Hardy Z-function, also known as the Riemann-Siegel Z-function + `Z(t) = e^{i \theta(t)} L(1/2+it)`, which is real-valued for real *t*. + The first *len* terms in the Taylor expansion are written to the output. + diff --git a/doc/source/dirichlet.rst b/doc/source/dirichlet.rst index b630b820..723934f6 100644 --- a/doc/source/dirichlet.rst +++ b/doc/source/dirichlet.rst @@ -243,7 +243,7 @@ encoding the zero value. .. function:: ulong dirichlet_chi(const dirichlet_group_t G, const dirichlet_char_t chi, ulong n) - Compute that value `\chi(n)` as the exponent modulo *G->expo*. + Compute the value `\chi(n)` as the exponent modulo *G->expo*. .. function:: void dirichlet_chi_vec(ulong * v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv) @@ -264,7 +264,7 @@ Character operations .. function:: void dirichlet_char_pow(dirichlet_char_t c, const dirichlet_group_t G, const dirichlet_char_t a, ulong n) - Take the power of some character. + Take the power of a character. .. function:: void dirichlet_char_lift(dirichlet_char_t chi_G, const dirichlet_group_t G, const dirichlet_char_t chi_H, const dirichlet_group_t H) @@ -278,3 +278,4 @@ Character operations This requires `c(\chi_G)\mid q_H\mid q_G`, where `c(\chi_G)` is the conductor of `\chi_G` and `q_G, q_H` are the moduli of G and H. + diff --git a/examples/real_roots.c b/examples/real_roots.c index a7f37cdf..dacc71ac 100644 --- a/examples/real_roots.c +++ b/examples/real_roots.c @@ -3,25 +3,51 @@ #include #include "arb_calc.h" #include "acb_hypgeom.h" +#include "acb_dirichlet.h" #include "flint/profiler.h" slong eval_count = 0; +typedef struct +{ + dirichlet_group_t * G; + dirichlet_char_t * chi; +} +z_param_struct; + int z_function(arb_ptr out, const arb_t inp, void * params, slong order, slong prec) { - arb_struct x[2]; + z_param_struct * par = params; - arb_init(x); - arb_init(x + 1); + if (par->G == NULL) + { + arb_struct x[2]; - arb_set(x, inp); - arb_one(x + 1); + arb_init(x); + arb_init(x + 1); - _arb_poly_riemann_siegel_z_series(out, x, FLINT_MIN(2, order), order, prec); + arb_set(x, inp); + arb_one(x + 1); - arb_clear(x); - arb_clear(x + 1); + _arb_poly_riemann_siegel_z_series(out, x, FLINT_MIN(2, order), order, prec); + + arb_clear(x); + arb_clear(x + 1); + } + else + { + acb_ptr tmp; + slong k; + + tmp = _acb_vec_init(order); + acb_set_arb(tmp, inp); + acb_dirichlet_hardy_z(tmp, tmp, *(par->G), *(par->chi), order, prec); + for (k = 0; k < order; k++) + arb_set(out + k, acb_realref(tmp + k)); + + _acb_vec_clear(tmp, order); + } eval_count++; return 0; @@ -152,7 +178,11 @@ int main(int argc, char *argv[]) arf_interval_ptr blocks; arb_calc_func_t function; int * info; - int params; + void * params; + int param1; + z_param_struct param2; + dirichlet_group_t G; + dirichlet_char_t chi; slong digits, low_prec, high_prec, i, num, found_roots, found_unknown; slong maxdepth, maxeval, maxfound; int refine; @@ -166,7 +196,7 @@ int main(int argc, char *argv[]) flint_printf("real_roots function a b [-refine d] [-verbose] " "[-maxdepth n] [-maxeval n] [-maxfound n] [-prec n]\n"); flint_printf("available functions:\n"); - flint_printf(" 0 Z(x), Riemann-Siegel Z-function\n"); + flint_printf(" 0 Z(x), Z-function (Riemann zeta or Dirichlet L-function)\n"); flint_printf(" 1 sin(x)\n"); flint_printf(" 2 sin(x^2)\n"); flint_printf(" 3 sin(1/x)\n"); @@ -174,15 +204,20 @@ int main(int argc, char *argv[]) flint_printf(" 5 Ai'(x), Airy function\n"); flint_printf(" 6 Bi(x), Airy function\n"); flint_printf(" 7 Bi'(x), Airy function\n"); + flint_printf("With 0, specify optional Dirichlet character with [-character q n]\n"); return 1; } - params = 0; + param1 = 0; + param2.G = NULL; + param2.chi = NULL; + params = ¶m1; switch (atoi(argv[1])) { case 0: function = z_function; + params = ¶m2; break; case 1: function = sin_x; @@ -195,19 +230,19 @@ int main(int argc, char *argv[]) break; case 4: function = airy; - params = 0; + param1 = 0; break; case 5: function = airy; - params = 1; + param1 = 1; break; case 6: function = airy; - params = 2; + param1 = 2; break; case 7: function = airy; - params = 3; + param1 = 3; break; default: flint_printf("require a function 0-7\n"); @@ -257,6 +292,14 @@ int main(int argc, char *argv[]) { low_prec = atol(argv[i+1]); } + else if (!strcmp(argv[i], "-character")) + { + dirichlet_group_init(G, atol(argv[i+1])); + dirichlet_char_init(chi, G); + dirichlet_char_log(chi, G, atol(argv[i+2])); + param2.G = &G; + param2.chi = χ + } } high_prec = digits * 3.32192809488736 + 10; @@ -280,7 +323,7 @@ int main(int argc, char *argv[]) TIMEIT_ONCE_START num = arb_calc_isolate_roots(&blocks, &info, function, - ¶ms, interval, maxdepth, maxeval, maxfound, low_prec); + params, interval, maxdepth, maxeval, maxfound, low_prec); for (i = 0; i < num; i++) { @@ -302,7 +345,7 @@ int main(int argc, char *argv[]) continue; if (arb_calc_refine_root_bisect(t, - function, ¶ms, blocks + i, 5, low_prec) + function, params, blocks + i, 5, low_prec) != ARB_CALC_SUCCESS) { flint_printf("warning: some bisection steps failed!\n"); @@ -316,7 +359,7 @@ int main(int argc, char *argv[]) } if (arb_calc_refine_root_bisect(blocks + i, - function, ¶ms, t, 5, low_prec) + function, params, t, 5, low_prec) != ARB_CALC_SUCCESS) { flint_printf("warning: some bisection steps failed!\n"); @@ -330,10 +373,10 @@ int main(int argc, char *argv[]) } arf_interval_get_arb(v, t, high_prec); - arb_calc_newton_conv_factor(C, function, ¶ms, v, low_prec); + arb_calc_newton_conv_factor(C, function, params, v, low_prec); arf_interval_get_arb(w, blocks + i, high_prec); - if (arb_calc_refine_root_newton(z, function, ¶ms, + if (arb_calc_refine_root_newton(z, function, params, w, v, C, 10, high_prec) != ARB_CALC_SUCCESS) { flint_printf("warning: some newton steps failed!\n"); @@ -357,6 +400,12 @@ int main(int argc, char *argv[]) flint_free(blocks); flint_free(info); + if (param2.G != NULL) + { + dirichlet_group_clear(G); + dirichlet_char_clear(chi); + } + arf_interval_clear(t); arf_interval_clear(interval); arf_clear(C);