mirror of
https://github.com/vale981/arb
synced 2025-03-05 17:31:38 -05:00
Merge pull request #290 from p15-git-acc/approx_ws_derivative
interpolation derivative for zeros refinement under limited precision
This commit is contained in:
Fredrik Johansson
GitHub
commit
e46c12c066
9 changed files with 447 additions and 49 deletions
|
|
@ -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);
|
||||
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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))
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
||||
|
|
|
|||
80
arb_poly/sinc_pi_series.c
Normal file
80
arb_poly/sinc_pi_series.c
Normal file
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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);
|
||||
}
|
||||
|
||||
104
arb_poly/test/t-sinc_pi_series.c
Normal file
104
arb_poly/test/t-sinc_pi_series.c
Normal file
|
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#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;
|
||||
}
|
||||
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue