interpolation derivative for zeros refinement under limited precision

acb_dirichlet.h

@@ -229,8 +229,8 @@ void acb_dirichlet_platt_c_precomp_init(acb_dirichlet_platt_c_precomp_t pre,
     slong sigma, const arb_t h, ulong k, slong prec);
 void acb_dirichlet_platt_c_precomp_clear(acb_dirichlet_platt_c_precomp_t pre);
 void acb_dirichlet_platt_c_bound_precomp(arb_t res,
-        const acb_dirichlet_platt_c_precomp_t pre, slong sigma, const arb_t t0,
+    const acb_dirichlet_platt_c_precomp_t pre, slong sigma, const arb_t t0,
-        const arb_t h, slong k, slong prec);
+    const arb_t h, slong k, slong prec);
 void acb_dirichlet_platt_c_bound(arb_t res,
     slong sigma, const arb_t t0, const arb_t h, slong k, slong prec);
 
@@ -249,15 +249,15 @@ void acb_dirichlet_platt_i_bound(arb_t res,
 /* Platt Gaussian-windowed Whittaker-Shannon interpolation */
 
 void acb_dirichlet_platt_ws_precomp_init(acb_dirichlet_platt_ws_precomp_t pre,
-        slong A, const arb_t H, slong sigma, slong prec);
+    slong A, const arb_t H, slong sigma, slong prec);
 void acb_dirichlet_platt_ws_precomp_clear(acb_dirichlet_platt_ws_precomp_t pre);
-void acb_dirichlet_platt_ws_interpolation_precomp(arb_t res,
+void acb_dirichlet_platt_ws_interpolation_precomp(arb_t res, arf_t deriv,
     const acb_dirichlet_platt_ws_precomp_t pre, const arb_t t0,
     arb_srcptr p, const fmpz_t T, slong A, slong B, slong Ns_max,
     const arb_t H, slong sigma, slong prec);
-void acb_dirichlet_platt_ws_interpolation(arb_t res, const arb_t t0,
+void acb_dirichlet_platt_ws_interpolation(arb_t res, arf_t deriv,
-    arb_srcptr p, const fmpz_t T, slong A, slong B, slong Ns_max,
+    const arb_t t0, arb_srcptr p, const fmpz_t T, slong A, slong B,
-    const arb_t H, slong sigma, slong prec);
+    slong Ns_max, const arb_t H, slong sigma, slong prec);
 void acb_dirichlet_platt_bound_C3(arb_t res, const arb_t t0, slong A,
     const arb_t H, slong Ns, slong prec);
 

acb_dirichlet/platt_local_hardy_z_zeros.c

@@ -82,21 +82,22 @@ platt_ctx_clear(platt_ctx_t ctx)
 }
 
 static void
-platt_ctx_interpolate(arb_t res,
+platt_ctx_interpolate(arb_t res, arf_t deriv,
         const platt_ctx_t ctx, const arb_t t0, slong prec)
 {
-    acb_dirichlet_platt_ws_interpolation_precomp(res, &ctx->pre, t0, ctx->p,
+    acb_dirichlet_platt_ws_interpolation_precomp(res, deriv,
-            &ctx->T, ctx->A, ctx->B, ctx->Ns_max, &ctx->H, ctx->sigma, prec);
+        &ctx->pre, t0, ctx->p, &ctx->T, ctx->A, ctx->B, ctx->Ns_max,
+        &ctx->H, ctx->sigma, prec);
 }
 
 static void
-platt_ctx_interpolate_arf(arb_t res,
+platt_ctx_interpolate_arf(arb_t res, arf_t deriv,
         const platt_ctx_t ctx, const arf_t t0, slong prec)
 {
     arb_t t;
     arb_init(t);
     arb_set_arf(t, t0);
-    platt_ctx_interpolate(res, ctx, t, prec);
+    platt_ctx_interpolate(res, deriv, ctx, t, prec);
     arb_clear(t);
 }
 
@@ -228,7 +229,7 @@ create_non_gram_node(const arf_t t, const platt_ctx_t ctx, slong prec)
     zz_node_ptr p = flint_malloc(sizeof(zz_node_struct));
     zz_node_init(p);
     arf_set(&p->t, t);
-    platt_ctx_interpolate_arf(&p->v, ctx, t, prec);
+    platt_ctx_interpolate_arf(&p->v, NULL, ctx, t, prec);
     if (arb_contains_zero(&p->v))
     {
         zz_node_clear(p);
@@ -258,7 +259,7 @@ create_gram_node(const fmpz_t n, const platt_ctx_t ctx, slong prec)
 
     acb_dirichlet_gram_point(t, n, NULL, NULL, prec + fmpz_sizeinbase(n, 2));
     acb_set_arb(z, t);
-    platt_ctx_interpolate(v, ctx, t, prec);
+    platt_ctx_interpolate(v, NULL, ctx, t, prec);
     if (!arb_contains_zero(v))
     {
         /* t contains g(n) and does not contain a zero of the f function */
@@ -1136,10 +1137,10 @@ _refine_local_hardy_z_zero_illinois(arb_t res,
     abs_tol = nmag - prec - 4;
 
     wp = prec + nmag + 8;
-    platt_ctx_interpolate_arf(z, ctx, a, wp);
+    platt_ctx_interpolate_arf(z, NULL, ctx, a, wp);
     asign = arb_sgn_nonzero(z);
     arf_set(fa, arb_midref(z));
-    platt_ctx_interpolate_arf(z, ctx, b, wp);
+    platt_ctx_interpolate_arf(z, NULL, ctx, b, wp);
     bsign = arb_sgn_nonzero(z);
     arf_set(fb, arb_midref(z));
 
@@ -1164,22 +1165,90 @@ _refine_local_hardy_z_zero_illinois(arb_t res,
         arf_mul(c, c, fa, wp, ARF_RND_NEAR);
         arf_sub(c, a, c, wp, ARF_RND_NEAR);
 
-        /* if c is not sandwiched between a and b, improve precision
+        /* if c is not sandwiched between a and b,
-           and fall back to one bisection step */
+           fall back to one bisection step */
         if (!arf_is_finite(c) ||
             !((arf_cmp(a, c) < 0 && arf_cmp(c, b) < 0) ||
               (arf_cmp(b, c) < 0 && arf_cmp(c, a) < 0)))
         {
             /* flint_printf("no sandwich (k = %wd)\n", k); */
-            wp += 32;
             arf_add(c, a, b, ARF_PREC_EXACT, ARF_RND_DOWN);
             arf_mul_2exp_si(c, c, -1);
         }
 
-        platt_ctx_interpolate_arf(z, ctx, c, wp);
+        platt_ctx_interpolate_arf(z, NULL, ctx, c, wp);
         csign = arb_sgn_nonzero(z);
+
+        /* If the guess is close enough to a zero that the sign
+         * cannot be determined, then use the derivative to
+         * make an appropriately small interval around the guess. */
         if (!csign)
+        {
+            arf_t deriv, aprime, bprime, faprime, fbprime, err, delta;
+            slong i, aprimesign, bprimesign;
+
+            arf_init(deriv);
+            arf_init(aprime);
+            arf_init(bprime);
+            arf_init(faprime);
+            arf_init(fbprime);
+            arf_init(err);
+            arf_init(delta);
+
+            arf_set_mag(err, arb_radref(z));
+            platt_ctx_interpolate_arf(NULL, deriv, ctx, c, wp);
+            arf_div(delta, err, deriv, wp, ARF_RND_NEAR);
+            arf_mul_si(delta, delta, 3, wp, ARF_RND_NEAR);
+            arf_mul_2exp_si(delta, delta, -1);
+            arf_set(aprime, c);
+            arf_set(bprime, c);
+
+            /* When the context allows the interval endpoints to
+             * be evaluated to relatively high precision,
+             * this should not require more than one or two iterations. */
+            for (i = 0; i < 5; i++)
+            {
+                arf_sub(aprime, aprime, delta, wp, ARF_RND_DOWN);
+                arf_add(bprime, bprime, delta, wp, ARF_RND_UP);
+                if (arf_cmp(a, b) < 0)
+                {
+                    if (arf_cmp(aprime, a) < 0)
+                        arf_set(aprime, a);
+                    if (arf_cmp(b, bprime) < 0)
+                        arf_set(bprime, b);
+                }
+                else
+                {
+                    if (arf_cmp(aprime, b) < 0)
+                        arf_set(aprime, b);
+                    if (arf_cmp(a, bprime) < 0)
+                        arf_set(bprime, a);
+                }
+                platt_ctx_interpolate_arf(z, NULL, ctx, aprime, wp);
+                arf_set(faprime, arb_midref(z));
+                aprimesign = arb_sgn_nonzero(z);
+                platt_ctx_interpolate_arf(z, NULL, ctx, bprime, wp);
+                arf_set(fbprime, arb_midref(z));
+                bprimesign = arb_sgn_nonzero(z);
+                if (aprimesign && bprimesign && aprimesign != bprimesign)
+                {
+                    arf_set(a, aprime);
+                    arf_set(b, bprime);
+                    arf_set(fa, faprime);
+                    arf_set(fb, fbprime);
+                    break;
+                }
+            }
+
+            arf_clear(deriv);
+            arf_clear(aprime);
+            arf_clear(bprime);
+            arf_clear(faprime);
+            arf_clear(fbprime);
+            arf_clear(err);
+            arf_clear(delta);
             break;
+        }
         arf_set(fc, arb_midref(z));
 
         if (csign != bsign)

acb_dirichlet/platt_ws_interpolation.c

@@ -335,29 +335,27 @@ finish:
     arb_clear(rhs);
 }
 
-
+/* Does not account for limited resolution and supporting points. */
 static void
-_interpolation_helper(arb_t res, const acb_dirichlet_platt_ws_precomp_t pre,
+_interpolation_helper_raw(arb_t res,
         const arb_t t0, arb_srcptr p, const fmpz_t T, slong A, slong B,
-        slong i0, slong Ns, const arb_t H, slong sigma, slong prec)
+        slong i0, slong Ns, const arb_t H, slong prec)
 {
     mag_t m;
     arb_t accum1; /* sum of terms where the argument of sinc is small */
     arb_t accum2; /* sum of terms where the argument of sinc is large */
-    arb_t total, dt0, dt, a, b, s, err, pi, g, x, c;
+    arb_t dt0, dt, a, b, s, pi, g, x, c;
     slong i;
     slong N = A*B;
 
     mag_init(m);
     arb_init(accum1);
     arb_init(accum2);
-    arb_init(total);
     arb_init(dt0);
     arb_init(dt);
     arb_init(a);
     arb_init(b);
     arb_init(s);
-    arb_init(err);
     arb_init(pi);
     arb_init(g);
     arb_init(x);
@@ -398,31 +396,191 @@ _interpolation_helper(arb_t res, const acb_dirichlet_platt_ws_precomp_t pre,
             arb_add(accum2, accum2, b, prec);
         }
     }
-    arb_set(total, accum1);
+    arb_set(res, accum1);
-    arb_addmul(total, accum2, c, prec);
+    arb_addmul(res, accum2, c, prec);
-    acb_dirichlet_platt_bound_C3(err, t0, A, H, Ns, prec);
-    arb_add_error(total, err);
-    acb_dirichlet_platt_i_bound_precomp(
-            err, &pre->pre_i, &pre->pre_c, t0, A, H, sigma, prec);
-    arb_add_error(total, err);
-    arb_set(res, total);
 
     mag_clear(m);
     arb_clear(accum1);
     arb_clear(accum2);
-    arb_clear(total);
     arb_clear(dt0);
     arb_clear(dt);
     arb_clear(a);
     arb_clear(b);
     arb_clear(s);
-    arb_clear(err);
     arb_clear(pi);
     arb_clear(g);
     arb_clear(x);
     arb_clear(c);
 }
 
+static void
+_arb_poly_sinc_pi_series(arb_ptr g, arb_srcptr h, slong hlen, slong n, slong prec)
+{
+    hlen = FLINT_MIN(hlen, n);
+    if (hlen == 1)
+    {
+        arb_sinc_pi(g, h, prec);
+        _arb_vec_zero(g + 1, n - 1);
+    }
+    else
+    {
+        arb_t pi;
+        arb_ptr t, u;
+
+        arb_init(pi);
+        t = _arb_vec_init(n + 1);
+        u = _arb_vec_init(hlen);
+
+        arb_const_pi(pi, prec);
+        _arb_vec_set(u, h, hlen);
+
+        if (arb_is_zero(h))
+        {
+            _arb_poly_sin_pi_series(t, u, hlen, n + 1, prec);
+            _arb_poly_div_series(g, t + 1, n, u + 1, hlen - 1, n, prec);
+        }
+        else
+        {
+            _arb_poly_sin_pi_series(t, u, hlen, n, prec);
+            _arb_poly_div_series(g, t, n, u, hlen, n, prec);
+        }
+        _arb_vec_scalar_div(g, g, n, pi, prec);
+
+        arb_clear(pi);
+        _arb_vec_clear(t, n + 1);
+        _arb_vec_clear(u, hlen);
+    }
+}
+
+/* Sets res to the function (a * exp(-(b-h)^2 / c)) * sinc_pi(d*(b-h)))
+ * of the power series h, for the purpose of computing derivatives
+ * of the Gaussian-windowed Whittaker-Shannon interpolation.
+ * Supports aliasing. */
+static void
+_arb_poly_gwws_series(arb_ptr res, arb_srcptr h, slong hlen,
+        const arb_t a, const arb_t b, const arb_t c, const arb_t d,
+        slong len, slong prec)
+{
+    arb_ptr u, u2, v, w;
+    hlen = FLINT_MIN(hlen, len);
+
+    u = _arb_vec_init(hlen);
+    u2 = _arb_vec_init(len);
+    v = _arb_vec_init(len);
+    w = _arb_vec_init(len);
+
+    /* u = b-h; u2 = (b-h)^2 */
+    _arb_vec_neg(u, h, hlen);
+    arb_add(u, u, b, prec);
+    _arb_poly_mullow(u2, u, hlen, u, hlen, len, prec);
+
+    /* v = exp(-(b-h)^2 / c) */
+    _arb_vec_scalar_div(v, u2, len, c, prec);
+    _arb_vec_neg(v, v, len);
+    _arb_poly_exp_series(v, v, len, len, prec);
+
+    /* w = sinc_pi(d*(b-h)) */
+    _arb_vec_scalar_mul(w, u, hlen, d, prec);
+    _arb_poly_sinc_pi_series(w, w, hlen, len, prec);
+
+    /* res = a * exp(-(b-h)^2 / c)) * sinc_pi(d*(b-h)) */
+    _arb_poly_mullow(res, v, len, w, len, len, prec);
+    _arb_vec_scalar_mul(res, res, len, a, prec);
+
+    _arb_vec_clear(u, hlen);
+    _arb_vec_clear(u2, len);
+    _arb_vec_clear(v, len);
+    _arb_vec_clear(w, len);
+}
+
+/* Does not account for limited resolution and supporting points. */
+static void
+_interpolation_helper_raw_series(arb_ptr res, arb_srcptr t0, slong t0len,
+        arb_srcptr p, const fmpz_t T, slong A, slong B, slong i0,
+        slong Ns, const arb_t H, slong trunc, slong prec)
+{
+    t0len = FLINT_MIN(t0len, trunc);
+    if (t0len == 1)
+    {
+        _interpolation_helper_raw(res, t0, p, T, A, B, i0, Ns, H, prec);
+        _arb_vec_zero(res + 1, trunc - 1);
+    }
+    else
+    {
+        arb_ptr h, g, accum;
+        arb_t b, c, d;
+        slong N = A*B;
+        slong i;
+
+        arb_init(b);
+        arb_init(c);
+        arb_init(d);
+        h = _arb_vec_init(t0len);
+        g = _arb_vec_init(trunc);
+        accum = _arb_vec_init(trunc);
+
+        arb_sqr(c, H, prec);
+        arb_mul_2exp_si(c, c, 1);
+        arb_set_si(d, A);
+        _arb_vec_set(h, t0, t0len);
+        arb_sub_fmpz(h, t0, T, prec + fmpz_clog_ui(T, 2));
+
+        for (i = i0; i < i0 + 2*Ns; i++)
+        {
+            slong n = i - N/2;
+            _arb_div_si_si(b, n, A, prec);
+            _arb_poly_gwws_series(g, h, t0len, p + i, b, c, d, trunc, prec);
+            _arb_vec_add(accum, accum, g, trunc, prec);
+        }
+        _arb_vec_set(res, accum, trunc);
+
+        arb_clear(b);
+        arb_clear(c);
+        arb_clear(d);
+        _arb_vec_clear(h, t0len);
+        _arb_vec_clear(g, trunc);
+        _arb_vec_clear(accum, trunc);
+    }
+}
+
+static void
+_interpolation_deriv_helper(arf_t res, const arb_t t0,
+        arb_srcptr p, const fmpz_t T, slong A, slong B, slong i0,
+        slong Ns, const arb_t H, slong prec)
+{
+    arb_ptr t, h;
+    t = _arb_vec_init(2);
+    h = _arb_vec_init(2);
+    arb_set(t+0, t0);
+    arb_one(t+1);
+    _interpolation_helper_raw_series(
+        h, t, 2, p, T, A, B, i0, Ns, H, 2, prec);
+    arf_set(res, arb_midref(h+1));
+    _arb_vec_clear(t, 2);
+    _arb_vec_clear(h, 2);
+}
+
+/* Accounts for limited resolution and supporting points. */
+static void
+_interpolation_helper(arb_t res, const acb_dirichlet_platt_ws_precomp_t pre,
+        const arb_t t0, arb_srcptr p, const fmpz_t T, slong A, slong B,
+        slong i0, slong Ns, const arb_t H, slong sigma, slong prec)
+{
+    arb_t total, err;
+    arb_init(total);
+    arb_init(err);
+    _interpolation_helper_raw(
+        total, t0, p, T, A, B, i0, Ns, H, prec);
+    acb_dirichlet_platt_bound_C3(err, t0, A, H, Ns, prec);
+    arb_add_error(total, err);
+    acb_dirichlet_platt_i_bound_precomp(
+        err, &pre->pre_i, &pre->pre_c, t0, A, H, sigma, prec);
+    arb_add_error(total, err);
+    arb_set(res, total);
+    arb_clear(total);
+    arb_clear(err);
+}
+
 
 void
 acb_dirichlet_platt_ws_precomp_init(acb_dirichlet_platt_ws_precomp_t pre,
@@ -439,7 +597,7 @@ acb_dirichlet_platt_ws_precomp_clear(acb_dirichlet_platt_ws_precomp_t pre)
     acb_dirichlet_platt_i_precomp_clear(&pre->pre_i);
 }
 
-void acb_dirichlet_platt_ws_interpolation_precomp(arb_t res,
+void acb_dirichlet_platt_ws_interpolation_precomp(arb_t res, arf_t deriv,
     const acb_dirichlet_platt_ws_precomp_t pre, const arb_t t0,
     arb_srcptr p, const fmpz_t T, slong A, slong B, slong Ns_max,
     const arb_t H, slong sigma, slong prec)
@@ -466,6 +624,10 @@ void acb_dirichlet_platt_ws_interpolation_precomp(arb_t res,
         arb_mul_si(dt0A, dt0, A, prec);
         arb_get_lbound_arf(lower_f, dt0A, prec);
         lower_n = arf_get_si(lower_f, ARF_RND_FLOOR);
+        if (deriv)
+        {
+            arf_zero(deriv);
+        }
 
         /*
          * More than one iteration is needed only when the set of
@@ -483,20 +645,31 @@ void acb_dirichlet_platt_ws_interpolation_precomp(arb_t res,
             else
             {
                 slong i0 = N/2 + n - (Ns - 1);
-                _interpolation_helper(
+                if (res)
-                        x, pre, t0, p, T, A, B, i0, Ns, H, sigma, prec);
-                if (n == lower_n)
                 {
-                    arb_set(total, x);
+                    _interpolation_helper(
+                            x, pre, t0, p, T, A, B, i0, Ns, H, sigma, prec);
+                    if (n == lower_n)
+                    {
+                        arb_set(total, x);
+                    }
+                    else
+                    {
+                        arb_union(total, total, x, prec);
+                    }
                 }
-                else
+                if (deriv)
                 {
-                    arb_union(total, total, x, prec);
+                    _interpolation_deriv_helper(
+                            deriv, t0, p, T, A, B, i0, Ns, H, prec);
                 }
             }
         }
 
-        arb_set(res, total);
+        if (res)
+        {
+            arb_set(res, total);
+        }
 
         arb_clear(x);
         arb_clear(dt0);
@@ -507,13 +680,13 @@ void acb_dirichlet_platt_ws_interpolation_precomp(arb_t res,
 }
 
 void
-acb_dirichlet_platt_ws_interpolation(arb_t res, const arb_t t0,
+acb_dirichlet_platt_ws_interpolation(arb_t res, arf_t deriv, const arb_t t0,
         arb_srcptr p, const fmpz_t T, slong A, slong B,
         slong Ns_max, const arb_t H, slong sigma, slong prec)
 {
     acb_dirichlet_platt_ws_precomp_t pre;
     acb_dirichlet_platt_ws_precomp_init(pre, A, H, sigma, prec);
     acb_dirichlet_platt_ws_interpolation_precomp(
-            res, pre, t0, p, T, A, B, Ns_max, H, sigma, prec);
+            res, deriv, pre, t0, p, T, A, B, Ns_max, H, sigma, prec);
     acb_dirichlet_platt_ws_precomp_clear(pre);
 }

acb_dirichlet/test/t-platt_ws_interpolation.c

@@ -65,7 +65,7 @@ int main()
             arb_abs(H, H);
 
             acb_dirichlet_platt_scaled_lambda(expected, t0, prec);
-            acb_dirichlet_platt_ws_interpolation(observed, t0, vec,
+            acb_dirichlet_platt_ws_interpolation(observed, NULL, t0, vec,
                     T, A, B, Ns_max, H, sigma, prec);
 
             if (!arb_overlaps(expected, observed))

doc/source/acb_dirichlet.rst

@@ -766,11 +766,12 @@ and formulas described by David J. Platt in [Pla2017]_.
     discrete Fourier transforms, and it requires the four additional tuning
     parameters *h*, *J*, *K*, and *sigma*.
 
-.. function:: void acb_dirichlet_platt_ws_interpolation(arb_t res, const arb_t t0, arb_srcptr p, const fmpz_t T, slong A, slong B, slong Ns_max, const arb_t H, slong sigma, slong prec)
+.. function:: void acb_dirichlet_platt_ws_interpolation(arb_t res, arf_t deriv, const arb_t t0, arb_srcptr p, const fmpz_t T, slong A, slong B, slong Ns_max, const arb_t H, slong sigma, slong prec)
 
     Compute :func:`acb_dirichlet_platt_scaled_lambda` at *t0* by
     Gaussian-windowed Whittaker-Shannon interpolation of points evaluated by
-    :func:`acb_dirichlet_platt_scaled_lambda_vec`.
+    :func:`acb_dirichlet_platt_scaled_lambda_vec`. The derivative is
+    also approximated if the output parameter *deriv* is not *NULL*.
     *Ns_max* defines the maximum number of supporting points to be used in
     the interpolation on either side of *t0*. *H* is the standard deviation
     of the Gaussian window centered on *t0* to be applied before the