From 3c219ee2e31aec9f2ea135a7542e299c2330bcdf Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Tue, 7 Jan 2025 15:38:02 +0100 Subject: [PATCH] allow evaluating deflated Riemann zeta function on small intervals containing 1 --- src/acb_poly/test/t-zeta_cpx_series.c | 57 +++++++++++++++++++++++++ src/acb_poly/zeta_series.c | 60 +++++++++++++++++++++++++++ 2 files changed, 117 insertions(+) diff --git a/src/acb_poly/test/t-zeta_cpx_series.c b/src/acb_poly/test/t-zeta_cpx_series.c index b074df0354..639d9349a1 100644 --- a/src/acb_poly/test/t-zeta_cpx_series.c +++ b/src/acb_poly/test/t-zeta_cpx_series.c @@ -83,5 +83,62 @@ TEST_FUNCTION_START(acb_poly_zeta_cpx_series, state) _acb_vec_clear(z2, len); } + /* check deflations around 1 */ + for (iter = 0; iter < 100 * 0.1 * flint_test_multiplier(); iter++) + { + acb_t s1, s2, a; + acb_ptr z1, z2; + slong i, len, prec; + int deflate; + + acb_init(s1); + acb_init(s2); + acb_init(a); + + acb_one(a); + + prec = 2 + n_randint(state, 300); + len = 1 + n_randint(state, 20); + deflate = 1; + + acb_randtest(s1, state, 300, 2); + acb_mul_2exp_si(s1, s1, -(slong) n_randint(state, 100)); + + acb_zero(s2); + arb_get_mag(arb_radref(acb_realref(s2)), acb_realref(s1)); + arb_get_mag(arb_radref(acb_imagref(s2)), acb_imagref(s1)); + + acb_add_ui(s1, s1, 1, prec); + acb_add_ui(s2, s2, 1, prec); + + z1 = _acb_vec_init(len); + z2 = _acb_vec_init(len); + + _acb_poly_zeta_cpx_series(z1, s1, a, deflate, len, prec); + _acb_poly_zeta_cpx_series(z2, s2, a, deflate, len, prec); + + for (i = 0; i < len; i++) + { + if (!acb_overlaps(z1 + i, z2 + i)) + { + flint_printf("FAIL: overlap (deflation)\n\n"); + flint_printf("iter = %wd\n", iter); + flint_printf("deflate = %d, len = %wd, i = %wd\n\n", deflate, len, i); + flint_printf("s1 = "); acb_printd(s1, prec / 3.33); flint_printf("\n\n"); + flint_printf("s2 = "); acb_printd(s2, prec / 3.33); flint_printf("\n\n"); + flint_printf("a = "); acb_printd(a, prec / 3.33); flint_printf("\n\n"); + flint_printf("z1 = "); acb_printd(z1 + i, prec / 3.33); flint_printf("\n\n"); + flint_printf("z2 = "); acb_printd(z2 + i, prec / 3.33); flint_printf("\n\n"); + flint_abort(); + } + } + + acb_clear(a); + acb_clear(s1); + acb_clear(s2); + _acb_vec_clear(z1, len); + _acb_vec_clear(z2, len); + } + TEST_FUNCTION_END(state); } diff --git a/src/acb_poly/zeta_series.c b/src/acb_poly/zeta_series.c index 59ebedec31..fcc188bf83 100644 --- a/src/acb_poly/zeta_series.c +++ b/src/acb_poly/zeta_series.c @@ -39,6 +39,66 @@ _acb_poly_zeta_cpx_series(acb_ptr z, const acb_t s, const acb_t a, int deflate, return; } + /* F(s) = zeta(s) - 1/(s-1) on small intervals containing 1 */ + if (deflate && acb_is_one(a) && !acb_is_one(s) && + (arb_contains_zero(acb_imagref(s)) && arb_contains_si(acb_realref(s), 1))) + { + acb_t t; + mag_t r, u, v; + slong n; + + /* An extremely crude bound is |gamma_n| <= 8^(-n) n!. + + (F^{(n)}(1+x) - F^{(n)}(1)) / n! = sum_{k>=1} F^{(n+k)} / (k! n!) + |(F^{(n)}(1+x) - F^{(n)}(1)) / n!| <= sum_{k>=1} x^k |gamma_{n+k}| / (k! n!) + <= sum_{k>=1} x^k (n+k)! / (k! n!) / 8^(n+k) + <= sum_{k>=1} x^k (n+1)^k / 8^(n+k) + = 8^(-n) (n+1) x / (8 - (n+1) x) + */ + + mag_init(r); + mag_init(u); + mag_init(v); + acb_init(t); + + is_real = acb_is_real(s); + + acb_sub_ui(t, s, 1, prec); + acb_get_mag(r, t); + mag_mul_ui(u, r, d); + + if (mag_cmp_2exp_si(u, 3) < 0) + { + acb_one(t); + _acb_poly_zeta_cpx_series(z, t, a, 1, d, prec); + + for (n = 0; n < d; n++) + { + mag_mul_ui(u, r, n + 1); + mag_one(v); + mag_mul_2exp_si(v, v, 3); + mag_sub_lower(v, v, u); + mag_div(v, u, v); + mag_mul_2exp_si(v, v, -3 * n); + + arb_add_error_mag(acb_realref(z + n), v); + if (!is_real) + arb_add_error_mag(acb_imagref(z + n), v); + } + } + else + { + _acb_vec_indeterminate(z, d); + } + + mag_clear(r); + mag_clear(u); + mag_clear(v); + acb_clear(t); + + return; + } + is_real = const_is_real = 0; if (acb_is_real(s) && acb_is_real(a))