Merge pull request flintlib#2139 from fredrik-johansson/sinc

sinc_series: allow intervals containing 0

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 <https://www.gnu.org/licenses/>.
 */
 
+#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);

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);
 }

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);