diff --git a/acb_dirichlet.h b/acb_dirichlet.h
index 733b4753..893c1a68 100644
--- a/acb_dirichlet.h
+++ b/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 arb_t h, slong k, slong prec);
+ const acb_dirichlet_platt_c_precomp_t pre, slong sigma, const arb_t t0,
+ 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,
- 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, 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);
void acb_dirichlet_platt_bound_C3(arb_t res, const arb_t t0, slong A,
const arb_t H, slong Ns, slong prec);
diff --git a/acb_dirichlet/platt_local_hardy_z_zeros.c b/acb_dirichlet/platt_local_hardy_z_zeros.c
index 451c70dc..e3c80646 100644
--- a/acb_dirichlet/platt_local_hardy_z_zeros.c
+++ b/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,
- &ctx->T, ctx->A, ctx->B, ctx->Ns_max, &ctx->H, ctx->sigma, prec);
+ acb_dirichlet_platt_ws_interpolation_precomp(res, deriv,
+ &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
- and fall back to one bisection step */
+ /* if c is not sandwiched between a and b,
+ 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)
diff --git a/acb_dirichlet/platt_ws_interpolation.c b/acb_dirichlet/platt_ws_interpolation.c
index 11d7c482..d93aa993 100644
--- a/acb_dirichlet/platt_ws_interpolation.c
+++ b/acb_dirichlet/platt_ws_interpolation.c
@@ -11,6 +11,7 @@
#include "acb_dirichlet.h"
#include "arb_hypgeom.h"
+#include "arb_poly.h"
/* Increase precision adaptively. */
static void
@@ -335,29 +336,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 +397,152 @@ _interpolation_helper(arb_t res, const acb_dirichlet_platt_ws_precomp_t pre,
arb_add(accum2, accum2, b, prec);
}
}
- arb_set(total, accum1);
- arb_addmul(total, 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);
+ arb_set(res, accum1);
+ arb_addmul(res, accum2, c, prec);
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);
}
+/* 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 +559,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 +586,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 +607,31 @@ void acb_dirichlet_platt_ws_interpolation_precomp(arb_t res,
else
{
slong i0 = N/2 + n - (Ns - 1);
- _interpolation_helper(
- x, pre, t0, p, T, A, B, i0, Ns, H, sigma, prec);
- if (n == lower_n)
+ if (res)
{
- 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 +642,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);
}
diff --git a/acb_dirichlet/test/t-platt_ws_interpolation.c b/acb_dirichlet/test/t-platt_ws_interpolation.c
index bbee32f6..d682bfef 100644
--- a/acb_dirichlet/test/t-platt_ws_interpolation.c
+++ b/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))
diff --git a/arb_poly.h b/arb_poly.h
index 03464c29..60da6ad0 100644
--- a/arb_poly.h
+++ b/arb_poly.h
@@ -636,6 +636,9 @@ void arb_poly_tan_series(arb_poly_t g, const arb_poly_t h, slong n, slong prec);
void _arb_poly_sinc_series(arb_ptr g, arb_srcptr h, slong hlen, slong n, slong prec);
void arb_poly_sinc_series(arb_poly_t g, const arb_poly_t h, slong n, slong prec);
+void _arb_poly_sinc_pi_series(arb_ptr g, arb_srcptr h, slong hlen, slong n, slong prec);
+void arb_poly_sinc_pi_series(arb_poly_t g, const arb_poly_t h, slong n, slong prec);
+
void _arb_poly_compose_series_brent_kung(arb_ptr res, arb_srcptr poly1, slong len1,
arb_srcptr poly2, slong len2, slong n, slong prec);
diff --git a/arb_poly/sinc_pi_series.c b/arb_poly/sinc_pi_series.c
new file mode 100644
index 00000000..c91f8cc8
--- /dev/null
+++ b/arb_poly/sinc_pi_series.c
@@ -0,0 +1,80 @@
+/*
+ Copyright (C) 2016 Fredrik Johansson
+ Copyright (C) 2019 D.H.J. Polymath
+
+ This file is part of Arb.
+
+ Arb is free software: you can redistribute it and/or modify it under
+ the terms of the GNU Lesser General Public License (LGPL) as published
+ by the Free Software Foundation; either version 2.1 of the License, or
+ (at your option) any later version. See .
+*/
+
+#include "arb_poly.h"
+
+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);
+ }
+}
+
+void
+arb_poly_sinc_pi_series(arb_poly_t g, const arb_poly_t h, slong n, slong prec)
+{
+ slong hlen = h->length;
+
+ if (n == 0)
+ {
+ arb_poly_zero(g);
+ return;
+ }
+
+ if (hlen == 0)
+ {
+ arb_poly_one(g);
+ return;
+ }
+
+ if (hlen == 1)
+ n = 1;
+
+ arb_poly_fit_length(g, n);
+ _arb_poly_sinc_pi_series(g->coeffs, h->coeffs, hlen, n, prec);
+ _arb_poly_set_length(g, n);
+ _arb_poly_normalise(g);
+}
+
diff --git a/arb_poly/test/t-sinc_pi_series.c b/arb_poly/test/t-sinc_pi_series.c
new file mode 100644
index 00000000..1beb075a
--- /dev/null
+++ b/arb_poly/test/t-sinc_pi_series.c
@@ -0,0 +1,104 @@
+/*
+ Copyright (C) 2012, 2013 Fredrik Johansson
+ Copyright (C) 2019 D.H.J. Polymath
+
+ This file is part of Arb.
+
+ Arb is free software: you can redistribute it and/or modify it under
+ the terms of the GNU Lesser General Public License (LGPL) as published
+ by the Free Software Foundation; either version 2.1 of the License, or
+ (at your option) any later version. See .
+*/
+
+#include "arb_poly.h"
+
+int main()
+{
+ slong iter;
+ flint_rand_t state;
+
+ flint_printf("sinc_pi_series....");
+ fflush(stdout);
+
+ flint_randinit(state);
+
+ for (iter = 0; iter < 200 * arb_test_multiplier(); iter++)
+ {
+ slong m, n1, n2, rbits1, rbits2, rbits3, rbits4;
+ arb_poly_t a, b, c, d;
+ arb_t pi;
+
+ rbits1 = 2 + n_randint(state, 300);
+ rbits2 = 2 + n_randint(state, 300);
+ rbits3 = 2 + n_randint(state, 300);
+ rbits4 = 2 + n_randint(state, 300);
+
+ m = n_randint(state, 15);
+ n1 = n_randint(state, 15);
+ n2 = n_randint(state, 15);
+
+ arb_poly_init(a);
+ arb_poly_init(b);
+ arb_poly_init(c);
+ arb_poly_init(d);
+ arb_init(pi);
+
+ arb_poly_randtest(a, state, m, rbits1, 10);
+ arb_poly_randtest(b, state, 10, rbits1, 10);
+ arb_poly_randtest(c, state, 10, rbits1, 10);
+
+ arb_poly_sinc_pi_series(b, a, n1, rbits2);
+ arb_poly_sinc_pi_series(c, a, n2, rbits3);
+
+ arb_poly_set(d, b);
+ arb_poly_truncate(d, FLINT_MIN(n1, n2));
+ arb_poly_truncate(c, FLINT_MIN(n1, n2));
+
+ arb_const_pi(pi, rbits4);
+
+ if (!arb_poly_overlaps(c, d))
+ {
+ flint_printf("FAIL\n\n");
+ flint_printf("n1 = %wd, n2 = %wd, bits2 = %wd, bits3 = %wd, bits4 = %wd\n", n1, n2, rbits2, rbits3, rbits4);
+ flint_printf("a = "); arb_poly_printd(a, 50); flint_printf("\n\n");
+ flint_printf("b = "); arb_poly_printd(b, 50); flint_printf("\n\n");
+ flint_printf("c = "); arb_poly_printd(c, 50); flint_printf("\n\n");
+ flint_abort();
+ }
+
+ /* check pi x sinc_pi(x) = sin_pi(x) */
+ arb_poly_mullow(c, b, a, n1, rbits2);
+ arb_poly_scalar_mul(c, c, pi, rbits2);
+ arb_poly_sin_pi_series(d, a, n1, rbits2);
+
+ if (!arb_poly_overlaps(c, d))
+ {
+ flint_printf("FAIL (functional equation)\n\n");
+ flint_printf("a = "); arb_poly_printd(a, 15); flint_printf("\n\n");
+ flint_printf("b = "); arb_poly_printd(b, 15); flint_printf("\n\n");
+ flint_printf("c = "); arb_poly_printd(c, 15); flint_printf("\n\n");
+ flint_printf("d = "); arb_poly_printd(d, 15); flint_printf("\n\n");
+ flint_abort();
+ }
+
+ arb_poly_sinc_pi_series(a, a, n1, rbits2);
+
+ if (!arb_poly_overlaps(a, b))
+ {
+ flint_printf("FAIL (aliasing)\n\n");
+ flint_abort();
+ }
+
+ arb_poly_clear(a);
+ arb_poly_clear(b);
+ arb_poly_clear(c);
+ arb_poly_clear(d);
+ arb_clear(pi);
+ }
+
+ flint_randclear(state);
+ flint_cleanup();
+ flint_printf("PASS\n");
+ return EXIT_SUCCESS;
+}
+
diff --git a/doc/source/acb_dirichlet.rst b/doc/source/acb_dirichlet.rst
index ea620626..dbc5271b 100644
--- a/doc/source/acb_dirichlet.rst
+++ b/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
diff --git a/doc/source/arb_poly.rst b/doc/source/arb_poly.rst
index a1ec6e33..74708b82 100644
--- a/doc/source/arb_poly.rst
+++ b/doc/source/arb_poly.rst
@@ -983,6 +983,12 @@ Powers and elementary functions
Sets *c* to the sinc function of the power series *h*, truncated
to length *n*.
+.. function:: void _arb_poly_sinc_pi_series(arb_ptr s, arb_srcptr h, slong hlen, slong n, slong prec)
+
+.. function:: void arb_poly_sinc_pi_series(arb_poly_t s, const arb_poly_t h, slong n, slong prec)
+
+ Compute the sinc function of the input multiplied by `\pi`.
+
Lambert W function
-------------------------------------------------------------------------------