From eee8b433a1d943ff720c55d771d31beb37cdcd74 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Mon, 19 Jun 2017 11:39:34 +0200 Subject: [PATCH] add test code for acb_dirichlet_hurwitz; more special case handling --- acb_dirichlet/hurwitz.c | 16 +++++++ acb_dirichlet/test/t-hurwitz.c | 84 ++++++++++++++++++++++++++++++++++ acb_dirichlet/zeta.c | 6 +++ 3 files changed, 106 insertions(+) create mode 100644 acb_dirichlet/test/t-hurwitz.c diff --git a/acb_dirichlet/hurwitz.c b/acb_dirichlet/hurwitz.c index 400746be..5b1614a7 100644 --- a/acb_dirichlet/hurwitz.c +++ b/acb_dirichlet/hurwitz.c @@ -29,6 +29,22 @@ acb_dirichlet_hurwitz(acb_t res, const acb_t s, const acb_t a, slong prec) return; } + if (acb_is_int(s) && arf_sgn(arb_midref(acb_realref(s))) < 0 && + arf_cmpabs_ui(arb_midref(acb_realref(s)), prec / 2) < 0) + { + slong n = arf_get_si(arb_midref(acb_realref(s)), ARF_RND_FLOOR); + + acb_bernoulli_poly_ui(res, 1 - n, a, prec); + acb_div_si(res, res, n - 1, prec); + return; + } + + if (arb_contains_zero(acb_imagref(s)) && arb_contains_si(acb_realref(s), 1)) + { + acb_indeterminate(res); + return; + } + _acb_poly_zeta_cpx_series(res, s, a, 0, 1, prec); } diff --git a/acb_dirichlet/test/t-hurwitz.c b/acb_dirichlet/test/t-hurwitz.c new file mode 100644 index 00000000..c3787a81 --- /dev/null +++ b/acb_dirichlet/test/t-hurwitz.c @@ -0,0 +1,84 @@ +/* + Copyright (C) 2017 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("hurwitz...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 3000 * arb_test_multiplier(); iter++) + { + acb_t s, a, b, z1, z2; + slong prec1, prec2; + + acb_init(s); + acb_init(a); + acb_init(b); + acb_init(z1); + acb_init(z2); + + prec1 = 2 + n_randint(state, 200); + prec2 = 2 + n_randint(state, 200); + + acb_randtest(s, state, 1 + n_randint(state, 500), 3); + acb_randtest(a, state, 1 + n_randint(state, 500), 3); + acb_randtest(b, state, 1 + n_randint(state, 500), 10); + acb_randtest(z1, state, 1 + n_randint(state, 500), 10); + acb_randtest(z2, state, 1 + n_randint(state, 500), 10); + + if (n_randint(state, 2)) + acb_set_si(a, n_randint(state, 100) - 50); + if (n_randint(state, 2)) + acb_set_si(s, n_randint(state, 100) - 50); + + /* test zeta(s,a) = 2^(-s) (zeta(s,a/2) + zeta(s,(a+1)/2)) */ + acb_mul_2exp_si(b, a, -1); + acb_dirichlet_hurwitz(z2, s, b, prec2); + acb_add_ui(b, a, 1, prec2); + acb_mul_2exp_si(b, b, -1); + acb_dirichlet_hurwitz(b, s, b, prec2); + acb_add(z2, z2, b, prec2); + acb_neg(b, s); + acb_set_ui(z1, 2); + acb_pow(b, z1, b, prec2); + acb_mul(z2, z2, b, prec2); + + acb_dirichlet_hurwitz(z1, s, a, prec1); + + if (!acb_overlaps(z1, z2)) + { + flint_printf("FAIL: overlap\n\n"); + flint_printf("s = "); acb_printd(s, 30); flint_printf("\n\n"); + flint_printf("a = "); acb_printd(a, 30); flint_printf("\n\n"); + flint_printf("z1 = "); acb_printd(z1, 30); flint_printf("\n\n"); + flint_printf("z2 = "); acb_printd(z2, 30); flint_printf("\n\n"); + flint_abort(); + } + + acb_clear(s); + acb_clear(a); + acb_clear(b); + acb_clear(z1); + acb_clear(z2); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/acb_dirichlet/zeta.c b/acb_dirichlet/zeta.c index 67b4c916..26046229 100644 --- a/acb_dirichlet/zeta.c +++ b/acb_dirichlet/zeta.c @@ -37,6 +37,12 @@ acb_dirichlet_zeta(acb_t res, const acb_t s, slong prec) return; } + if (arb_contains_zero(acb_imagref(s)) && arb_contains_si(acb_realref(s), 1)) + { + acb_indeterminate(res); + return; + } + cutoff = 24.0 * prec * sqrt(prec); if (arf_cmpabs_d(arb_midref(acb_imagref(s)), cutoff) >= 0 &&