diff --git a/src/acb_poly/sinc_series.c b/src/acb_poly/sinc_series.c index af5d667c93..48081357ad 100644 --- a/src/acb_poly/sinc_series.c +++ b/src/acb_poly/sinc_series.c @@ -1,5 +1,5 @@ /* - Copyright (C) 2016 Fredrik Johansson + Copyright (C) 2016, 2025 Fredrik Johansson This file is part of FLINT. @@ -9,8 +9,118 @@ (at your option) any later version. See . */ +#include "mag.h" #include "acb_poly.h" +void +_acb_sinc_jet_zero(acb_ptr res, const acb_t z, slong len, slong prec) +{ + mag_ptr D; + slong n; + int is_real, is_imag; + mag_t zm, err, fac; + arb_t wide; + + is_real = acb_is_real(z); + is_imag = arb_is_zero(acb_realref(z)); + + mag_init(zm); + mag_init(err); + mag_init(fac); + arb_init(wide); + + acb_get_mag(zm, z); + + D = _mag_vec_init(len + 1); + + /* Compute sequence of derivative bounds + + sinc^{(n)}(z) = int_0^1 t^n cos(t z + n pi / 2) dt + + |sinc^{(n)}(z)| <= cosh(im(z)) min(1/|z|, 1/(n+1)) + |sinc^{(n)}(z)| <= cosh(im(z)) min(1/|z|, 1/(n+1), |z|/(n+2)) (odd n) + */ + { + mag_t C, b1, b2, b3; + + mag_init(C); + mag_init(b1); + mag_init(b2); + mag_init(b3); + + /* C = cosh(im(z)) */ + arb_get_mag(C, acb_imagref(z)); + mag_cosh(C, C); + + /* b1 = 1/|z| */ + acb_get_mag_lower(b1, z); + mag_inv(b1, b1); + + for (n = 0; n <= len; n++) + { + /* b2 = 1/(n+1) */ + mag_one(b2); + mag_div_ui(b2, b2, n + 1); + mag_min(b2, b1, b2); + + /* b3 = |z|/(n+2) */ + if (n % 2 == 1) + { + mag_div_ui(b3, zm, n + 2); + mag_min(b2, b2, b3); + } + + mag_mul(D + n, C, b2); + } + + mag_clear(C); + mag_clear(b1); + mag_clear(b2); + mag_clear(b3); + } + + mag_one(fac); + + for (n = 0; n < len; n++) + { + /* sinc^{(n)}(0) / n! */ + if (n == 0) + acb_one(res + n); + else if (n % 2 == 1) + acb_zero(res + n); + else + acb_div_si(res + n, res + n - 2, -n * (n + 1), prec); + + if (n > 1) + mag_div_ui(fac, fac, n); + + /* |sinc^{(n)}(0 + eps) - sinc^{(n)}(0)| / n! <= |eps| * |sinc^{(n+1)}(z)| / n!, z = (+/- eps) */ + mag_mul(err, zm, D + n + 1); + mag_mul(err, err, fac); + acb_add_error_mag(res + n, err); + + /* sinc^{(n)}(+/- eps) is a possibly better enclosure */ + arb_zero(wide); + mag_mul(err, D + n, fac); + arb_add_error_mag(wide, err); + arb_intersection(acb_realref(res + n), acb_realref(res + n), wide, prec); + arb_intersection(acb_imagref(res + n), acb_imagref(res + n), wide, prec); + + if (is_real || (is_imag && (n % 2 == 0))) + arb_zero(acb_imagref(res + n)); + + if (is_imag && (n % 2 == 1)) + arb_zero(acb_realref(res + n)); + } + + _mag_vec_clear(D, len + 1); + + mag_clear(zm); + mag_clear(err); + mag_clear(fac); + arb_clear(wide); +} + void _acb_poly_sinc_series(acb_ptr g, acb_srcptr h, slong hlen, slong n, slong prec) { @@ -35,6 +145,13 @@ _acb_poly_sinc_series(acb_ptr g, acb_srcptr h, slong hlen, slong n, slong prec) _acb_poly_sin_series(t, u, hlen, n + 1, prec); _acb_poly_div_series(g, t + 1, n, u + 1, hlen - 1, n, prec); } + else if (acb_contains_zero(h)) + { + _acb_sinc_jet_zero(t, h, n, prec); + /* compose with nonconstant part */ + acb_zero(u); + _acb_poly_compose_series(g, t, n, u, hlen, n, prec); + } else { _acb_poly_sin_series(t, u, hlen, n, prec); diff --git a/src/acb_poly/test/t-sinc_series.c b/src/acb_poly/test/t-sinc_series.c index 2d204279c7..aa554d7532 100644 --- a/src/acb_poly/test/t-sinc_series.c +++ b/src/acb_poly/test/t-sinc_series.c @@ -83,5 +83,59 @@ TEST_FUNCTION_START(acb_poly_sinc_series, state) acb_poly_clear(d); } + /* check intervals containing 0 */ + for (iter = 0; iter < 100 * flint_test_multiplier(); iter++) + { + acb_t z1, z2; + acb_poly_t f1, f2, g1, g2; + slong len, prec; + + acb_init(z1); + acb_init(z2); + + acb_poly_init(f1); + acb_poly_init(f2); + acb_poly_init(g1); + acb_poly_init(g2); + + prec = 2 + n_randint(state, 300); + len = 1 + n_randint(state, 5); + + acb_randtest(z1, state, 300, 5); + acb_mul_2exp_si(z1, z1, -(slong) n_randint(state, 100)); + if (n_randint(state, 2)) + arb_zero(acb_realref(z1)); + + arb_get_mag(arb_radref(acb_realref(z2)), acb_realref(z1)); + arb_get_mag(arb_radref(acb_imagref(z2)), acb_imagref(z1)); + + acb_poly_set_coeff_si(f1, 1, n_randint(state, 2) ? 1 : -1); + acb_poly_set(f2, f1); + acb_poly_set_coeff_acb(f1, 0, z1); + acb_poly_set_coeff_acb(f2, 0, z2); + + acb_poly_sinc_series(g1, f1, len, prec); + acb_poly_sinc_series(g2, f2, len, prec); + + if (!acb_poly_overlaps(g1, g2)) + { + flint_printf("FAIL: overlap (0)\n\n"); + flint_printf("iter = %wd\n", iter); + flint_printf("f1 = "); acb_poly_printd(f1, prec / 3.33); flint_printf("\n\n"); + flint_printf("f2 = "); acb_poly_printd(f2, prec / 3.33); flint_printf("\n\n"); + flint_printf("g1 = "); acb_poly_printd(g1, prec / 3.33); flint_printf("\n\n"); + flint_printf("g2 = "); acb_poly_printd(g2, prec / 3.33); flint_printf("\n\n"); + flint_abort(); + } + + acb_clear(z1); + acb_clear(z2); + + acb_poly_clear(f1); + acb_poly_clear(f2); + acb_poly_clear(g1); + acb_poly_clear(g2); + } + TEST_FUNCTION_END(state); } diff --git a/src/arb_poly/sinc_series.c b/src/arb_poly/sinc_series.c index 62d0891cf0..6d2060e252 100644 --- a/src/arb_poly/sinc_series.c +++ b/src/arb_poly/sinc_series.c @@ -10,6 +10,23 @@ */ #include "arb_poly.h" +#include "acb_poly.h" + +void _acb_sinc_jet_zero(acb_ptr res, const acb_t z, slong len, slong prec); + +void +_arb_sinc_jet_zero(arb_ptr res, const arb_t z, slong len, slong prec) +{ + acb_ptr t; + slong i; + t = _acb_vec_init(len + 1); + acb_set_arb(t + len, z); + _acb_sinc_jet_zero(t, t + len, len, prec); + for (i = 0; i < len; i++) + arb_swap(res + i, acb_realref(t + i)); + _acb_vec_clear(t, len + 1); + +} void _arb_poly_sinc_series(arb_ptr g, arb_srcptr h, slong hlen, slong n, slong prec) @@ -35,6 +52,13 @@ _arb_poly_sinc_series(arb_ptr g, arb_srcptr h, slong hlen, slong n, slong prec) _arb_poly_sin_series(t, u, hlen, n + 1, prec); _arb_poly_div_series(g, t + 1, n, u + 1, hlen - 1, n, prec); } + else if (arb_contains_zero(h)) + { + _arb_sinc_jet_zero(t, h, n, prec); + /* compose with nonconstant part */ + arb_zero(u); + _arb_poly_compose_series(g, t, n, u, hlen, n, prec); + } else { _arb_poly_sin_series(t, u, hlen, n, prec);