diff --git a/acb_dirichlet.h b/acb_dirichlet.h index f2eca095..0a6e15e7 100644 --- a/acb_dirichlet.h +++ b/acb_dirichlet.h @@ -45,6 +45,8 @@ 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 diff --git a/acb_dirichlet/riemann_xi.c b/acb_dirichlet/riemann_xi.c new file mode 100644 index 00000000..de3e9245 --- /dev/null +++ b/acb_dirichlet/riemann_xi.c @@ -0,0 +1,51 @@ +/* + 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-riemann_xi.c new file mode 100644 index 00000000..da2a5fe1 --- /dev/null +++ b/acb_dirichlet/test/t-riemann_xi.c @@ -0,0 +1,74 @@ +/* + 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" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("riemann_xi...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 3000 * arb_test_multiplier(); iter++) + { + acb_t a, b, c; + slong prec1, prec2; + + prec1 = 2 + n_randint(state, 1000); + prec2 = prec1 + 30; + + acb_init(a); + acb_init(b); + acb_init(c); + + 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); + + if (!acb_overlaps(b, c)) + { + flint_printf("FAIL: overlap\n\n"); + flint_printf("a = "); acb_print(a); flint_printf("\n\n"); + flint_printf("b = "); acb_print(b); flint_printf("\n\n"); + flint_printf("c = "); acb_print(c); flint_printf("\n\n"); + flint_abort(); + } + + /* check riemann_xi(s) = riemann_xi(1-s) */ + acb_sub_ui(c, a, 1, prec1); + acb_neg(c, c); + acb_riemann_xi(c, c, prec1); + + if (!acb_overlaps(b, c)) + { + flint_printf("FAIL: functional equation\n\n"); + flint_printf("a = "); acb_print(a); flint_printf("\n\n"); + flint_printf("b = "); acb_print(b); flint_printf("\n\n"); + flint_printf("c = "); acb_print(c); flint_printf("\n\n"); + flint_abort(); + } + + acb_clear(a); + acb_clear(b); + acb_clear(c); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +}