From b924a3afcd05f2a021d27f85823027c93e0c0336 Mon Sep 17 00:00:00 2001 From: fredrik Date: Sat, 27 Nov 2021 13:29:53 +0100 Subject: [PATCH] acb_dirichlet_l_fmpq_afe: handle the singular series --- acb_dirichlet/l_fmpq_afe.c | 103 +++++++++++++------- acb_dirichlet/test/t-l_fmpq_afe.c | 6 +- arb_hypgeom.h | 2 + arb_hypgeom/gamma_upper_fmpq.c | 134 ++++++++++++++++++++++++++ arb_hypgeom/test/t-gamma_upper_fmpq.c | 56 +++++++++++ doc/source/arb_hypgeom.rst | 12 +++ 6 files changed, 277 insertions(+), 36 deletions(-) diff --git a/acb_dirichlet/l_fmpq_afe.c b/acb_dirichlet/l_fmpq_afe.c index ee937c63..45a972bb 100644 --- a/acb_dirichlet/l_fmpq_afe.c +++ b/acb_dirichlet/l_fmpq_afe.c @@ -135,13 +135,12 @@ acb_dirichlet_fmpq_sum_afe(acb_t res, const fmpq_t s, const dirichlet_group_t G, arb_t ns, t, u, v, z, z0, z1, x, x2, Ga, Gz1, Gz0, expmz0, z0_prevn, Gz0_prevn, expmz0_prevn; acb_t c; fmpq_t s2; - int parity; + int parity, gamma_singular; ulong q; double abs_tol_mag; double gammainc_mag, gamma_mag, ns_mag; double aa, zz; - double cancellation; #if VERBOSE double t1, t2, t3; @@ -188,7 +187,11 @@ acb_dirichlet_fmpq_sum_afe(acb_t res, const fmpq_t s, const dirichlet_group_t G, /* s2 = (s+parity)/2 */ fmpq_add_ui(s2, s, parity); fmpq_div_2exp(s2, s2, 1); - arb_gamma_fmpq(Ga, s2, gamma_cached_prec); + + gamma_singular = (fmpz_is_one(fmpq_denref(s2)) && fmpz_sgn(fmpq_numref(s2)) <= 0); + + if (!gamma_singular) + arb_gamma_fmpq(Ga, s2, gamma_cached_prec); for (n = 1; ; n += 1) { @@ -222,28 +225,34 @@ acb_dirichlet_fmpq_sum_afe(acb_t res, const fmpq_t s, const dirichlet_group_t G, /* Gamma((s+parity)/2, z) (want lower bound, to estimate cancellation) */ gammainc_mag = log_gamma_upper_approx(aa, zz) / log(2); - /* Gamma((s+parity)/2) */ - gamma_mag = ARF_EXP(arb_midref(Ga)); /* n^-s */ ns_mag = -fmpq_get_d(s) * log(n) / log(2); - cancellation = FLINT_MAX(gamma_mag - gammainc_mag, 0); - cancellation = cancellation; /* suppress warning when VERBOSE == 0 */ /* Want Gamma(a,z) n^-s with abs_tol --> want Gamma(a,z) with abs_tol * n^s */ mag_set_ui_2exp_si(abs_tol_gamma, 1, abs_tol_mag - ns_mag); - /* Precision needed sans cancellation. */ + /* wp = Precision needed sans cancellation. */ wp = gammainc_mag + ns_mag - abs_tol_mag + 5; wp = FLINT_MAX(wp, 30); - /* Precision needed with cancellation. */ - wp2 = FLINT_MAX(gamma_mag, gammainc_mag) + ns_mag - abs_tol_mag + 5; - wp2 = FLINT_MAX(wp2, 30); + /* wp2 = Precision needed with cancellation. */ + if (gamma_singular) + { + /* Max term is roughly n^(-s) * z^((s+parity)/2) * exp(z) */ + wp2 = -abs_tol_mag + ns_mag + aa * log(zz) + zz / log(2) + 5; + wp2 = FLINT_MAX(wp2, 30); + } + else + { + /* Estimate of Gamma((s+parity)/2) */ + gamma_mag = ARF_EXP(arb_midref(Ga)); + wp2 = FLINT_MAX(gamma_mag, gammainc_mag) + ns_mag - abs_tol_mag + 5; + wp2 = FLINT_MAX(wp2, 30); + } #if VERBOSE printf(" abs_tol_gamma = "); mag_printd(abs_tol_gamma, 5); printf("\n"); printf(" gamma(a) = "); arb_printd(Ga, 10); printf("\n"); - printf(" cancellation = %f\n", cancellation); printf(" wp = %ld wp2 = %ld\n", wp, wp2); #endif @@ -286,32 +295,54 @@ acb_dirichlet_fmpq_sum_afe(acb_t res, const fmpq_t s, const dirichlet_group_t G, else { /* Otherwise fallback to the series at 0. */ - NN = _arb_hypgeom_gamma_lower_fmpq_0_choose_N(AE, s2, z0, abs_tol_gamma); - -#if VERBOSE - flint_printf(" lower series with N = %wd, z0 = ", NN); arb_printd(z0, 10); printf(" "); - mag_printd(AE, 10); printf("\n"); -#endif - - _arb_hypgeom_gamma_lower_fmpq_0_bsplit(Gz0, s2, z0, NN, wp2); - arb_add_error_mag(Gz0, AE); - - while (mag_cmp(arb_radref(Ga), abs_tol_gamma) > 0) + if (gamma_singular) { - gamma_cached_prec *= 2; - arb_gamma_fmpq(Ga, s2, gamma_cached_prec); + slong nn; + + nn = *fmpq_numref(s2); + + if (COEFF_IS_MPZ(nn)) + { + arb_indeterminate(Gz0); + } + else + { + nn = -nn; + + NN = _arb_hypgeom_gamma_upper_singular_si_choose_N(AE, nn, z0, abs_tol_gamma); + _arb_hypgeom_gamma_upper_singular_si_bsplit(Gz0, nn, z0, NN, wp2); + arb_add_error_mag(Gz0, AE); + +#if VERBOSE + flint_printf(" singular series with N = %wd, z0 = ", NN); arb_printd(z0, 10); printf(" "); + mag_printd(AE, 10); printf("\n"); +#endif + } } - + else + { + NN = _arb_hypgeom_gamma_lower_fmpq_0_choose_N(AE, s2, z0, abs_tol_gamma); + _arb_hypgeom_gamma_lower_fmpq_0_bsplit(Gz0, s2, z0, NN, wp2); + arb_add_error_mag(Gz0, AE); #if VERBOSE - flint_printf(" lower series with N = %wd: ", NN); arb_printd(Gz0, 10); printf("\n"); + flint_printf(" lower series with N = %wd, z0 = ", NN); arb_printd(z0, 10); printf(" "); + mag_printd(AE, 10); printf("\n"); #endif - arb_sub(Gz0, Ga, Gz0, wp); - + while (mag_cmp(arb_radref(Ga), abs_tol_gamma) > 0) + { + gamma_cached_prec *= 2; + arb_gamma_fmpq(Ga, s2, gamma_cached_prec); + } #if VERBOSE - flint_printf(" G(a) - lower series: "); arb_printd(Gz0, 10); printf("\n"); + flint_printf(" lower series with N = %wd: ", NN); arb_printd(Gz0, 10); printf("\n"); #endif + arb_sub(Gz0, Ga, Gz0, wp); +#if VERBOSE + flint_printf(" G(a) - lower series: "); arb_printd(Gz0, 10); printf("\n"); +#endif + } } } else @@ -437,17 +468,19 @@ acb_dirichlet_l_fmpq_afe(acb_t res, const fmpq_t s, const dirichlet_group_t G, c q = (G == NULL) ? 1 : G->q; parity = (G == NULL) ? 0 : dirichlet_parity_char(G, chi); - /* Limit at gamma(-n) is not implemented. */ + /* Division by gamma((s+parity)/2) at a pole. */ if (fmpz_is_one(fmpq_denref(s))) { const fmpz * n = fmpq_numref(s); if ((parity == 0 && fmpz_sgn(n) <= 0 && fmpz_is_even(n)) || - (parity == 0 && fmpz_sgn(n) > 0 && fmpz_is_odd(n)) || - (parity == 1 && fmpz_sgn(n) < 0 && fmpz_is_odd(n)) || - (parity == 1 && fmpz_sgn(n) > 0 && fmpz_is_even(n))) + (parity == 1 && fmpz_sgn(n) < 0 && fmpz_is_odd(n))) { - acb_indeterminate(res); + /* Special case for zeta */ + if (q == 1 && fmpz_is_zero(n)) + acb_set_d(res, -0.5); + else + acb_zero(res); return; } } diff --git a/acb_dirichlet/test/t-l_fmpq_afe.c b/acb_dirichlet/test/t-l_fmpq_afe.c index e5100d66..910bc0d4 100644 --- a/acb_dirichlet/test/t-l_fmpq_afe.c +++ b/acb_dirichlet/test/t-l_fmpq_afe.c @@ -45,7 +45,11 @@ int main() if (dirichlet_char_is_primitive(G, chi)) { - fmpq_randtest(x, state, 2 + n_randint(state, 8)); + if (n_randint(state, 2)) + fmpq_set_si(x, -8 + (slong) n_randint(state, 8), 1); + else + fmpq_randtest(x, state, 2 + n_randint(state, 8)); + acb_set_fmpq(s, x, prec); if (n_randint(state, 2)) diff --git a/arb_hypgeom.h b/arb_hypgeom.h index 833b6689..615425fb 100644 --- a/arb_hypgeom.h +++ b/arb_hypgeom.h @@ -177,6 +177,8 @@ void _arb_hypgeom_gamma_upper_fmpq_inf_bsplit(arb_t res, const fmpq_t a, const a slong _arb_hypgeom_gamma_lower_fmpq_0_choose_N(mag_t err, const fmpq_t a, const arb_t z, const mag_t abs_tol); void _arb_hypgeom_gamma_lower_fmpq_0_bsplit(arb_t res, const fmpq_t a, const arb_t z, slong N, slong prec); void _arb_gamma_upper_fmpq_step_bsplit(arb_t Gz1, const fmpq_t a, const arb_t z0, const arb_t z1, const arb_t Gz0, const arb_t expmz0, const mag_t abs_tol, slong prec); +slong _arb_hypgeom_gamma_upper_singular_si_choose_N(mag_t err, slong n, const arb_t z, const mag_t abs_tol); +void _arb_hypgeom_gamma_upper_singular_si_bsplit(arb_t res, slong n, const arb_t z, slong N, slong prec); #ifdef __cplusplus } diff --git a/arb_hypgeom/gamma_upper_fmpq.c b/arb_hypgeom/gamma_upper_fmpq.c index bf420c46..659fffef 100644 --- a/arb_hypgeom/gamma_upper_fmpq.c +++ b/arb_hypgeom/gamma_upper_fmpq.c @@ -351,3 +351,137 @@ _arb_hypgeom_gamma_lower_fmpq_0_bsplit(arb_t res, const fmpq_t a, const arb_t z, arb_clear(S); arb_clear(Q); } + +/* bounded by 1/z^n * sum z^k / k! */ +slong +_arb_hypgeom_gamma_upper_singular_si_choose_N(mag_t err, slong n, const arb_t z, const mag_t abs_tol) +{ + slong k; + mag_t t, u, zm; + + mag_init(t); + mag_init(u); + mag_init(zm); + + arb_get_mag(zm, z); + + arb_get_mag_lower(t, z); + mag_inv(t, t); + mag_pow_ui(t, t, n); + + for (k = 1; ; k++) + { + mag_mul(t, t, zm); + mag_div_ui(t, t, k); + + if (mag_cmp(t, abs_tol) < 0) + { + mag_div_ui(u, zm, k); + mag_geom_series(u, u, 0); + mag_mul(u, t, u); + + if (mag_cmp(u, abs_tol) < 0) + { + mag_swap(err, t); + break; + } + } + } + + mag_clear(t); + mag_clear(u); + mag_clear(zm); + + return k; +} + +static void +singular_bsplit(arb_t M, arb_t S, arb_t Q, slong n, const arb_t z, slong na, slong nb, int cont, slong prec) +{ + if (nb - na == 0) + { + arb_zero(M); + arb_zero(S); + arb_one(Q); + } + else if (nb - na == 1) + { + fmpz_t t; + slong k; + + k = na; + + if (k == n) + arb_neg(M, z); + else + arb_mul_si(M, z, n - k, prec); + + arb_set_si(S, (k != n) ? (k + 1) : 0); + + fmpz_init_set_si(t, k + 1); + if (k != n) + fmpz_mul_si(t, t, k - n); + arb_set_fmpz(Q, t); + fmpz_clear(t); + } + else + { + slong m; + arb_t M2, S2, Q2; + + m = na + (nb - na) / 2; + + arb_init(M2); + arb_init(S2); + arb_init(Q2); + + singular_bsplit(M, S, Q, n, z, na, m, 1, prec); + singular_bsplit(M2, S2, Q2, n, z, m, nb, cont, prec); + + arb_mul(S, S, Q2, prec); + arb_addmul(S, M, S2, prec); + + if (cont) + arb_mul(M, M, M2, prec); + + arb_mul(Q, Q, Q2, prec); + + arb_clear(M2); + arb_clear(S2); + arb_clear(Q2); + } +} + +void +_arb_hypgeom_gamma_upper_singular_si_bsplit(arb_t res, slong n, const arb_t z, slong N, slong prec) +{ + arb_t M, S, Q; + + arb_init(M); + arb_init(S); + arb_init(Q); + + N = FLINT_MAX(N, 0); + + singular_bsplit(M, S, Q, n, z, 0, N, 0, prec); + + /* (-1)**n/fac(n) * (digamma(n+1) - ln(z)) - (S/Q)/z**n */ + arb_pow_ui(M, z, n, prec); + arb_mul(Q, Q, M, prec); + arb_div(S, S, Q, prec); + + arb_set_ui(M, n + 1); + arb_digamma(M, M, prec); + arb_log(Q, z, prec); + arb_sub(M, M, Q, prec); + arb_fac_ui(Q, n, prec); + arb_div(M, M, Q, prec); + if (n & 1) + arb_neg(M, M); + + arb_sub(res, M, S, prec); + + arb_clear(M); + arb_clear(S); + arb_clear(Q); +} diff --git a/arb_hypgeom/test/t-gamma_upper_fmpq.c b/arb_hypgeom/test/t-gamma_upper_fmpq.c index b09e5b65..9c2988c6 100644 --- a/arb_hypgeom/test/t-gamma_upper_fmpq.c +++ b/arb_hypgeom/test/t-gamma_upper_fmpq.c @@ -141,6 +141,62 @@ int main() mag_clear(err2); } + for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++) + { + slong na; + arb_t z, r1, r2; + mag_t tol1, tol2, err1, err2; + slong N1, N2, prec1, prec2; + + mag_init(tol1); + mag_init(tol2); + mag_init(err1); + mag_init(err2); + + arb_init(z); + arb_init(r1); + arb_init(r2); + + prec1 = 2 + n_randint(state, 100); + prec2 = 2 + n_randint(state, 100); + + na = n_randint(state, 50); + + arb_set_ui(z, 1 + n_randint(state, 100)); + arb_div_ui(z, z, 1 + n_randint(state, 100), 200); + + mag_set_ui_2exp_si(tol1, 1, -(slong) n_randint(state, 100)); + mag_set_ui_2exp_si(tol2, 1, -(slong) n_randint(state, 100)); + + N1 = _arb_hypgeom_gamma_upper_singular_si_choose_N(err1, na, z, tol1); + N2 = _arb_hypgeom_gamma_upper_singular_si_choose_N(err2, na, z, tol2); + + _arb_hypgeom_gamma_upper_singular_si_bsplit(r1, na, z, N1, prec1); + arb_add_error_mag(r1, err1); + + _arb_hypgeom_gamma_upper_singular_si_bsplit(r2, na, z, N2, prec2); + arb_add_error_mag(r2, err2); + + if (!arb_overlaps(r1, r2)) + { + flint_printf("FAIL: overlap\n\n"); + flint_printf("na = %wd", na); flint_printf("\n\n"); + flint_printf("z = "); arb_printn(z, 100, 0); flint_printf("\n\n"); + flint_printf("r1 = "); arb_printn(r1, 100, 0); flint_printf("\n\n"); + flint_printf("r2 = "); arb_printn(r2, 100, 0); flint_printf("\n\n"); + flint_abort(); + } + + arb_clear(z); + arb_clear(r1); + arb_clear(r2); + + mag_clear(tol1); + mag_clear(tol2); + mag_clear(err1); + mag_clear(err2); + } + flint_randclear(state); flint_cleanup(); flint_printf("PASS\n"); diff --git a/doc/source/arb_hypgeom.rst b/doc/source/arb_hypgeom.rst index 9a773ef3..ee83dd6c 100644 --- a/doc/source/arb_hypgeom.rst +++ b/doc/source/arb_hypgeom.rst @@ -313,6 +313,18 @@ Internal evaluation functions the Taylor series at zero before term *N*. The truncation error bound has to be added separately. +.. function:: slong _arb_hypgeom_gamma_upper_singular_si_choose_N(mag_t err, slong n, const arb_t z, const mag_t abs_tol) + + Returns number of terms *N* and sets *err* to the truncation error for evaluating + `\Gamma(-n,z)` using the Taylor series at zero, targeting an absolute + tolerance of *abs_tol*. + +.. function:: void _arb_hypgeom_gamma_upper_singular_si_bsplit(arb_t res, slong n, const arb_t z, slong N, slong prec) + + Sets *res* to the approximation of `\Gamma(-n,z)` obtained by truncating + the Taylor series at zero before term *N*. + The truncation error bound has to be added separately. + .. function:: void _arb_gamma_upper_fmpq_step_bsplit(arb_t Gz1, const fmpq_t a, const arb_t z0, const arb_t z1, const arb_t Gz0, const arb_t expmz0, const mag_t abs_tol, slong prec) Given *Gz0* and *expmz0* representing the values `\Gamma(a,z_0)` and `\exp(-z_0)`,