diff --git a/acb_dirichlet.h b/acb_dirichlet.h index 0a6e15e7..7e2338dc 100644 --- a/acb_dirichlet.h +++ b/acb_dirichlet.h @@ -45,8 +45,6 @@ void acb_dirichlet_zeta(acb_t res, const acb_t s, slong prec); void acb_dirichlet_zeta_jet_rs(acb_ptr res, const acb_t s, slong len, slong prec); void acb_dirichlet_zeta_jet(acb_t res, const acb_t s, int deflate, slong len, slong prec); -void acb_riemann_xi(acb_t res, const acb_t s, slong prec); - void acb_dirichlet_hurwitz(acb_t res, const acb_t s, const acb_t a, slong prec); typedef struct @@ -75,6 +73,8 @@ void _acb_dirichlet_euler_product_real_ui(arb_t res, ulong s, void acb_dirichlet_eta(acb_t res, const acb_t s, slong prec); +void acb_dirichlet_xi(acb_t res, const acb_t s, slong prec); + void acb_dirichlet_pairing(acb_t res, const dirichlet_group_t G, ulong m, ulong n, slong prec); void acb_dirichlet_pairing_char(acb_t res, const dirichlet_group_t G, const dirichlet_char_t a, const dirichlet_char_t b, slong prec); diff --git a/acb_dirichlet/riemann_xi.c b/acb_dirichlet/riemann_xi.c deleted file mode 100644 index de3e9245..00000000 --- a/acb_dirichlet/riemann_xi.c +++ /dev/null @@ -1,51 +0,0 @@ -/* - Copyright (C) 2018 D.H.J Polymath - - 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" - -void acb_riemann_xi(acb_t res, const acb_t s, slong prec) -{ - acb_t pi, z1, z2, z3, z4; - - acb_init(pi); - acb_const_pi(pi, prec); - - /* s*(s-1)/2 */ - acb_init(z1); - acb_sub_ui(z1, s, 1, prec); - acb_mul(z1, z1, s, prec); - acb_mul_2exp_si(z1, z1, -1); - - /* pi^(-s/2) */ - acb_init(z2); - acb_mul_2exp_si(z2, s, -1); - acb_neg(z2, z2); - acb_pow(z2, pi, z2, prec); - - /* gamma(s/2) */ - acb_init(z3); - acb_mul_2exp_si(z3, s, -1); - acb_gamma(z3, z3, prec); - - /* zeta(s) */ - acb_init(z4); - acb_zeta(z4, s, prec); - - acb_mul(res, z1, z2, prec); - acb_mul(res, res, z3, prec); - acb_mul(res, res, z4, prec); - - acb_clear(pi); - acb_clear(z1); - acb_clear(z2); - acb_clear(z3); - acb_clear(z4); -} diff --git a/acb_dirichlet/test/t-riemann_xi.c b/acb_dirichlet/test/t-xi.c similarity index 75% rename from acb_dirichlet/test/t-riemann_xi.c rename to acb_dirichlet/test/t-xi.c index da2a5fe1..a3819631 100644 --- a/acb_dirichlet/test/t-riemann_xi.c +++ b/acb_dirichlet/test/t-xi.c @@ -16,7 +16,7 @@ int main() slong iter; flint_rand_t state; - flint_printf("riemann_xi...."); + flint_printf("xi...."); fflush(stdout); flint_randinit(state); @@ -36,8 +36,8 @@ int main() arb_randtest_precise(acb_realref(a), state, 1 + n_randint(state, 1000), 3); arb_randtest_precise(acb_imagref(a), state, 1 + n_randint(state, 1000), 3); - acb_riemann_xi(b, a, prec1); - acb_riemann_xi(c, a, prec2); + acb_dirichlet_xi(b, a, prec1); + acb_dirichlet_xi(c, a, prec2); if (!acb_overlaps(b, c)) { @@ -48,10 +48,22 @@ int main() flint_abort(); } - /* check riemann_xi(s) = riemann_xi(1-s) */ + acb_set(c, a); + acb_dirichlet_xi(c, c, prec1); + + if (!acb_equal(b, c)) + { + flint_printf("FAIL: aliasing\n\n"); + flint_printf("a = "); acb_printd(a, 30); flint_printf("\n\n"); + flint_printf("b = "); acb_printd(b, 30); flint_printf("\n\n"); + flint_printf("c = "); acb_printd(c, 30); flint_printf("\n\n"); + flint_abort(); + } + + /* check xi(s) = xi(1-s) */ acb_sub_ui(c, a, 1, prec1); acb_neg(c, c); - acb_riemann_xi(c, c, prec1); + acb_dirichlet_xi(c, c, prec1); if (!acb_overlaps(b, c)) { diff --git a/acb_dirichlet/xi.c b/acb_dirichlet/xi.c new file mode 100644 index 00000000..26789be2 --- /dev/null +++ b/acb_dirichlet/xi.c @@ -0,0 +1,68 @@ +/* + Copyright (C) 2018 D.H.J Polymath + + 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" + +static void +_acb_dirichlet_xi(acb_t res, const acb_t s, slong prec) +{ + acb_t pi, z1, z2, z3; + + acb_init(pi); + acb_init(z1); + acb_init(z2); + acb_init(z3); + + /* (s-1) * pi^(-s/2) * gamma(1 + s/2) * zeta(s) */ + acb_const_pi(pi, prec); + acb_mul_2exp_si(z1, s, -1); + acb_neg(z1, z1); + acb_pow(z1, pi, z1, prec); + acb_mul_2exp_si(z2, s, -1); + acb_add_ui(z2, z2, 1, prec); + acb_gamma(z2, z2, prec); + acb_zeta(z3, s, prec); + acb_sub_ui(res, s, 1, prec); + acb_mul(res, res, z1, prec); + acb_mul(res, res, z2, prec); + acb_mul(res, res, z3, prec); + + acb_clear(pi); + acb_clear(z1); + acb_clear(z2); + acb_clear(z3); +} + +/* todo: interval containing 1 */ +void acb_dirichlet_xi(acb_t res, const acb_t s, slong prec) +{ + if (!acb_is_finite(s)) + { + acb_indeterminate(res); + } + else if (acb_is_one(s)) + { + acb_one(res); + acb_mul_2exp_si(res, res, -1); + } + else if (arf_sgn(arb_midref(acb_realref(s))) < 0) + { + acb_t t; + acb_sub_ui(t, s, 1, prec); + acb_neg(t, t); + _acb_dirichlet_xi(res, t, prec); + acb_clear(t); + } + else + { + _acb_dirichlet_xi(res, s, prec); + } +} diff --git a/doc/source/acb_dirichlet.rst b/doc/source/acb_dirichlet.rst index 24cf9aa3..b5e25916 100644 --- a/doc/source/acb_dirichlet.rst +++ b/doc/source/acb_dirichlet.rst @@ -123,6 +123,12 @@ Riemann zeta function Note that the alternating character `\{1,-1\}` is not itself a Dirichlet character. +.. function:: void acb_dirichlet_xi(acb_t res, const acb_t s, slong prec) + + Sets *res* to the Riemann xi function + `\xi(s) = \frac{1}{2} s (s-1) \pi^{-s/2} \Gamma(\frac{1}{2} s) \zeta(s)`. + The functional equation for xi is `\xi(1-s) = \xi(s)`. + Riemann-Siegel formula -------------------------------------------------------------------------------