hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-04 17:01:40 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

platt whittaker shannon interpolation

Browse source
This commit is contained in:
p15-git-acc 2019-04-19 00:05:14 -05:00
parent 56ce687ea1
commit 45c91a0051
7 changed files with 1092 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

70
acb_dirichlet.h
View file

@ -195,6 +195,76 @@ acb_dirichlet_zeta_zero(acb_t res, const fmpz_t n, slong prec)
acb_dirichlet_zeta_zeros(res, n, 1, prec);
}
/* Platt zeta zeros */
typedef struct
{
slong len;
arb_ptr p;
arb_struct Xa;
arb_struct Xb;
}
acb_dirichlet_platt_c_precomp_struct;
typedef acb_dirichlet_platt_c_precomp_struct acb_dirichlet_platt_c_precomp_t[1];
typedef struct
{
arb_struct c1;
arb_struct c2;
}
acb_dirichlet_platt_i_precomp_struct;
typedef acb_dirichlet_platt_i_precomp_struct acb_dirichlet_platt_i_precomp_t[1];
typedef struct
{
acb_dirichlet_platt_c_precomp_struct pre_c;
acb_dirichlet_platt_i_precomp_struct pre_i;
}
acb_dirichlet_platt_ws_precomp_struct;
typedef acb_dirichlet_platt_ws_precomp_struct acb_dirichlet_platt_ws_precomp_t[1];
/* Platt C bound */
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);
void acb_dirichlet_platt_c_bound(arb_t res,
slong sigma, const arb_t t0, const arb_t h, slong k, slong prec);
/* Platt I bound */
void acb_dirichlet_platt_i_precomp_init(acb_dirichlet_platt_i_precomp_t pre,
slong A, const arb_t H, slong sigma, slong prec);
void acb_dirichlet_platt_i_precomp_clear(acb_dirichlet_platt_i_precomp_t pre);
void acb_dirichlet_platt_i_bound_precomp(arb_t res,
const acb_dirichlet_platt_i_precomp_t pre_i,
const acb_dirichlet_platt_c_precomp_t pre_c,
const arb_t t0, slong A, const arb_t H, slong sigma, slong prec);
void acb_dirichlet_platt_i_bound(arb_t res,
const arb_t t0, slong A, const arb_t H, slong sigma, slong prec);
/* 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);
void acb_dirichlet_platt_ws_precomp_clear(acb_dirichlet_platt_ws_precomp_t pre);
void acb_dirichlet_platt_ws_interpolation_precomp(arb_t res,
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_bound_C3(arb_t res, const arb_t t0, slong A,
const arb_t H, slong Ns, slong prec);
void acb_dirichlet_platt_scaled_lambda(arb_t res, const arb_t t, slong prec);
void acb_dirichlet_platt_scaled_lambda_vec(arb_ptr res, const fmpz_t T,
slong A, slong B, slong prec);
/* Discrete Fourier Transform */
void acb_dirichlet_dft_index(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec);

253
acb_dirichlet/platt_c_bound.c Normal file
View file

@ -0,0 +1,253 @@
/*
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 "acb_dirichlet.h"
#include "arb_hypgeom.h"
static void
_arb_pow_si(arb_t res, const arb_t x, slong y, slong prec)
{
arb_t a;
arb_init(a);
arb_set_si(a, y);
arb_pow(res, x, a, prec);
arb_clear(a);
}
static void
_arb_add_d(arb_t res, const arb_t x, double y, slong prec)
{
arb_t a;
arb_init(a);
arb_set_d(a, y);
arb_add(res, x, a, prec);
arb_clear(a);
}
static void
_arb_sub_d(arb_t res, const arb_t x, double y, slong prec)
{
_arb_add_d(res, x, -y, prec);
}
/* res = 2^((6*k+5-sigma)/4) * pi^k * (sigma + 1/2)^k * h */
static void
_pre_c_Xa(arb_t res, slong sigma, const arb_t h, ulong k, slong prec)
{
arb_t pi, two, x1, x2;
arb_init(pi);
arb_init(two);
arb_init(x1);
arb_init(x2);
arb_const_pi(pi, prec);
arb_set_si(two, 2);
arb_set_si(x1, (6*k + 5 - sigma));
arb_mul_2exp_si(x1, x1, -2);
arb_pow(x1, two, x1, prec);
arb_set_si(x2, sigma);
_arb_add_d(x2, x2, 0.5, prec);
arb_mul(x2, x2, pi, prec);
arb_pow_ui(x2, x2, k, prec);
arb_mul(res, x1, x2, prec);
arb_mul(res, res, h, prec);
arb_clear(pi);
arb_clear(two);
arb_clear(x1);
arb_clear(x2);
}
/* res = 2^((6*k+7-sigma)/4) * pi^(k - 1/2) */
static void
_pre_c_Xb(arb_t res, slong sigma, ulong k, slong prec)
{
arb_t pi, two, x1, x2;
arb_init(pi);
arb_init(two);
arb_init(x1);
arb_init(x2);
arb_const_pi(pi, prec);
arb_set_si(two, 2);
arb_set_si(x1, (6*k + 7 - sigma));
arb_mul_2exp_si(x1, x1, -2);
arb_pow(x1, two, x1, prec);
arb_set_ui(x2, k);
_arb_sub_d(x2, x2, 0.5, prec);
arb_pow(x2, pi, x2, prec);
arb_mul(res, x1, x2, prec);
arb_clear(pi);
arb_clear(two);
arb_clear(x1);
arb_clear(x2);
}
/*
* res[(sigma-1)/2 - l] = binom((sigma-1)/2, l) * 2^((k+l-1)/2)
* * h^(k+l+1) * gamma((k+l+1)/2, (sigma+1/2)^2 / 2*h^2)
*/
static void
_pre_c_p(arb_ptr res, slong sigma, const arb_t h, ulong k, slong prec)
{
slong l;
slong len = (sigma - 1)/2 + 1;
arb_t two, f, f1, f2, u, base, x;
arb_init(two);
arb_init(f);
arb_init(f1);
arb_init(f2);
arb_init(u);
arb_init(base);
arb_init(x);
/* precompute 2^((k-1)/2) * h^(k+1) */
arb_set_ui(two, 2);
arb_set_si(f1, k-1);
arb_mul_2exp_si(f1, f1, -1);
arb_pow(f1, two, f1, prec);
_arb_pow_si(f2, h, k+1, prec);
arb_mul(f, f1, f2, prec);
/* precompute (sigma + 1/2)^2 / 2*h^2 */
arb_set_si(u, sigma);
_arb_add_d(u, u, 0.5, prec);
arb_div(u, u, h, prec);
arb_sqr(u, u, prec);
arb_mul_2exp_si(u, u, -1);
/* precompute powers of sqrt(2)*h */
arb_sqrt_ui(base, 2, prec);
arb_mul(base, base, h, prec);
_arb_vec_set_powers(res, base, len, prec);
for (l = 0; l < len; l++)
{
arb_mul(res + l, res + l, f, prec);
arb_bin_uiui(x, (sigma - 1)/2, l, prec);
arb_mul(res + l, res + l, x, prec);
arb_set_si(x, k + l + 1);
arb_mul_2exp_si(x, x, -1);
arb_hypgeom_gamma_upper(x, x, u, 0, prec);
arb_mul(res + l, res + l, x, prec);
}
_arb_poly_reverse(res, res, len, len);
arb_clear(two);
arb_clear(f);
arb_clear(f1);
arb_clear(f2);
arb_clear(u);
arb_clear(base);
arb_clear(x);
}
void
acb_dirichlet_platt_c_precomp_init(acb_dirichlet_platt_c_precomp_t pre,
slong sigma, const arb_t h, ulong k, slong prec)
{
if (!arb_is_positive(h))
{
flint_printf("requires positive h\n");
flint_abort();
}
else if (sigma % 2 == 0 || sigma < 3)
{
flint_printf("requires odd integer sigma >= 3 (sigma=%wd)\n", sigma);
flint_abort();
}
else
{
pre->len = (sigma - 1) / 2 + 1;
arb_init(&pre->Xa);
arb_init(&pre->Xb);
pre->p = _arb_vec_init(pre->len);
_pre_c_Xa(&pre->Xa, sigma, h, k, prec);
_pre_c_Xb(&pre->Xb, sigma, k, prec);
_pre_c_p(pre->p, sigma, h, k, prec);
}
}
void
acb_dirichlet_platt_c_precomp_clear(acb_dirichlet_platt_c_precomp_t pre)
{
arb_clear(&pre->Xa);
arb_clear(&pre->Xb);
_arb_vec_clear(pre->p, pre->len);
}
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)
{
/* requires sigma + 1/2 <= t0 */
arb_t lhs;
arb_init(lhs);
arb_set_si(lhs, sigma);
_arb_add_d(lhs, lhs, 0.5, prec);
if (!arb_le(lhs, t0))
{
arb_zero_pm_inf(res);
}
else
{
arb_t u, v;
arb_init(u);
arb_init(v);
/* u = exp((1 + sqrt(8))/(6*t0)) */
arb_sqrt_ui(u, 8, prec);
arb_add_ui(u, u, 1, prec);
arb_div_ui(u, u, 6, prec);
arb_div(u, u, t0, prec);
arb_exp(u, u, prec);
/* v = (sigma + 1/2 + t0)^((sigma - 1)/2) */
arb_add_si(v, t0, sigma, prec);
_arb_add_d(v, v, 0.5, prec);
_arb_pow_si(v, v, (sigma - 1)/2, prec);
/* res = u*(v*Xa + Xb*X(t0)) */
_arb_poly_evaluate(res, pre->p, pre->len, t0, prec);
arb_mul(res, res, &pre->Xb, prec);
arb_addmul(res, v, &pre->Xa, prec);
arb_mul(res, res, u, prec);
arb_clear(u);
arb_clear(v);
}
arb_clear(lhs);
}
void
acb_dirichlet_platt_c_bound(arb_t res,
slong sigma, const arb_t t0, const arb_t h, slong k, slong prec)
{
acb_dirichlet_platt_c_precomp_t pre;
acb_dirichlet_platt_c_precomp_init(pre, sigma, h, k, prec);
acb_dirichlet_platt_c_bound_precomp(res, pre, sigma, t0, h, k, prec);
acb_dirichlet_platt_c_precomp_clear(pre);
}

160
acb_dirichlet/platt_i_bound.c Normal file
View file

@ -0,0 +1,160 @@
/*
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 "acb_dirichlet.h"
#include "arb_hypgeom.h"
/* c1 = 4*zeta(sigma)/(2*sigma-1) * pi^((-3-2*sigma)/4) * exp(...) */
static void
_pre_i_c1(arb_t res, slong A, const arb_t H, slong sigma, slong prec)
{
arb_t pi, x1, x2, x3, x4;
arb_init(pi);
arb_init(x1);
arb_init(x2);
arb_init(x3);
arb_init(x4);
arb_const_pi(pi, prec);
/* x1 = 4*zeta(sigma) / (2*sigma - 1) */
arb_set_si(x1, sigma);
arb_zeta(x1, x1, prec);
arb_mul_2exp_si(x1, x1, 2);
arb_div_si(x1, x1, 2*sigma - 1, prec);
/* x2 = pi^((-3-2*sigma)/4) */
arb_inv(x2, pi, prec);
arb_root_ui(x2, x2, 4, prec);
arb_pow_ui(x2, x2, 3 + 2*sigma, prec);
/* x3 = (2*sigma-1)^2 / (8*H^2) */
arb_set_si(x3, 2*sigma - 1);
arb_div(x3, x3, H, prec);
arb_sqr(x3, x3, prec);
arb_mul_2exp_si(x3, x3, -3);
/* x4 = pi*A*(2*sigma-1)/2 */
arb_mul_si(x4, pi, A*(2*sigma-1), prec);
arb_mul_2exp_si(x4, x4, -1);
/* res = x1*x2*exp(x3 - x4) */
arb_sub(res, x3, x4, prec);
arb_exp(res, res, prec);
arb_mul(res, res, x1, prec);
arb_mul(res, res, x2, prec);
arb_clear(pi);
arb_clear(x1);
arb_clear(x2);
arb_clear(x3);
arb_clear(x4);
}
/* c2 = 8*pi^(1/4) */
static void
_pre_i_c2(arb_t res, slong prec)
{
arb_t pi;
arb_init(pi);
arb_const_pi(pi, prec);
arb_root_ui(res, pi, 4, prec);
arb_mul_2exp_si(res, res, 3);
arb_clear(pi);
}
void
acb_dirichlet_platt_i_precomp_init(acb_dirichlet_platt_i_precomp_t pre,
slong A, const arb_t H, slong sigma, slong prec)
{
arb_init(&pre->c1);
arb_init(&pre->c2);
_pre_i_c1(&pre->c1, A, H, sigma, prec);
_pre_i_c2(&pre->c2, prec);
}
void
acb_dirichlet_platt_i_precomp_clear(acb_dirichlet_platt_i_precomp_t pre)
{
arb_clear(&pre->c1);
arb_clear(&pre->c2);
}
void
acb_dirichlet_platt_i_bound_precomp(arb_t res,
const acb_dirichlet_platt_i_precomp_t pre_i,
const acb_dirichlet_platt_c_precomp_t pre_c,
const arb_t t0, slong A, const arb_t H,
slong sigma, slong prec)
{
arb_t pi, cbound, x1, x2, x3, x4, x5;
arb_init(pi);
arb_init(cbound);
arb_init(x1);
arb_init(x2);
arb_init(x3);
arb_init(x4);
arb_init(x5);
arb_const_pi(pi, prec);
/* x1 = (1 - 4*t0^2)/(8*H^2) */
arb_sqr(x1, t0, prec);
arb_mul_2exp_si(x1, x1, 2);
arb_sub_ui(x1, x1, 1, prec);
arb_neg(x1, x1);
arb_div(x1, x1, H, prec);
arb_div(x1, x1, H, prec);
arb_mul_2exp_si(x1, x1, -3);
/* x2 = pi*A/2 */
arb_mul_si(x2, pi, A, prec);
arb_mul_2exp_si(x2, x2, -1);
/* x3 = exp(x1 - x2) */
arb_sub(x3, x1, x2, prec);
arb_exp(x3, x3, prec);
acb_dirichlet_platt_c_bound_precomp(
cbound, pre_c, sigma, t0, H, 0, prec);
/* res = c1*cbound + c2*x3 */
arb_mul(x4, &pre_i->c1, cbound, prec);
arb_mul(x5, &pre_i->c2, x3, prec);
arb_add(res, x4, x5, prec);
arb_clear(pi);
arb_clear(cbound);
arb_clear(x1);
arb_clear(x2);
arb_clear(x3);
arb_clear(x4);
arb_clear(x5);
}
void
acb_dirichlet_platt_i_bound(arb_t res, const arb_t t0, slong A, const arb_t H,
slong sigma, slong prec)
{
acb_dirichlet_platt_c_precomp_t pre_c;
acb_dirichlet_platt_i_precomp_t pre_i;
acb_dirichlet_platt_c_precomp_init(pre_c, sigma, H, 0, prec);
acb_dirichlet_platt_i_precomp_init(pre_i, A, H, sigma, prec);
acb_dirichlet_platt_i_bound_precomp(
res, pre_i, pre_c, t0, A, H, sigma, prec);
acb_dirichlet_platt_c_precomp_clear(pre_c);
acb_dirichlet_platt_i_precomp_clear(pre_i);
}

465
acb_dirichlet/platt_ws_interpolation.c Normal file
View file

@ -0,0 +1,465 @@
/*
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 "acb_dirichlet.h"
#include "arb_hypgeom.h"
static void
_arb_inv_ui(arb_t res, ulong n, slong prec)
{
arb_set_ui(res, n);
arb_inv(res, res, prec);
}
static void
_arb_div_si_si(arb_t res, slong x, slong y, slong prec)
{
arb_set_si(res, x);
arb_div_si(res, res, y, prec);
}
static void
_arb_ui_pow_arb(arb_t res, ulong n, const arb_t x, slong prec)
{
arb_t a;
arb_init(a);
arb_set_ui(a, n);
arb_pow(res, a, x, prec);
arb_clear(a);
}
/*
* Follows the notation in https://en.wikipedia.org/wiki/Gaussian_function
* res = a*exp(-(x-b)^2 / (2*c^2))
* If 'a' is NULL a default coefficient of 1 is used.
* If 'b' is NULL a default mean of 0 is used.
*/
static void
_arb_gaussian(arb_t res, const arb_t a, const arb_t b, const arb_t c,
const arb_t x, slong prec)
{
arb_t z;
arb_init(z);
if (b == NULL)
arb_set(z, x);
else
arb_sub(z, x, b, prec);
arb_div(z, z, c, prec);
arb_sqr(z, z, prec);
arb_mul_2exp_si(z, z, -1);
arb_neg(z, z);
arb_exp(z, z, prec);
if (a != NULL)
arb_mul(res, z, a, prec);
arb_clear(z);
}
static void
_platt_lambda(arb_t res, const arb_t t, slong prec)
{
acb_t pi, s, z, s1, s2;
acb_init(pi);
acb_init(s);
acb_init(z);
acb_init(s1);
acb_init(s2);
arb_set_d(acb_realref(s), 0.5);
arb_set(acb_imagref(s), t);
acb_const_pi(pi, prec);
arb_mul_2exp_si(acb_imagref(s1), t, -1);
acb_neg(s1, s1);
acb_pow(s1, pi, s1, prec);
acb_mul_2exp_si(s2, s, -1);
acb_gamma(s2, s2, prec);
acb_zeta(z, s, prec);
acb_mul(z, z, s1, prec);
acb_mul(z, z, s2, prec);
if (!arb_contains_zero(acb_imagref(z)))
{
flint_abort();
}
arb_set(res, acb_realref(z));
acb_clear(pi);
acb_clear(s);
acb_clear(z);
acb_clear(s1);
acb_clear(s2);
}
void
acb_dirichlet_platt_scaled_lambda(arb_t res, const arb_t t, slong prec)
{
arb_t pi, lam;
arb_init(pi);
arb_init(lam);
arb_const_pi(pi, prec);
_platt_lambda(lam, t, prec);
arb_mul(res, pi, t, prec);
arb_mul_2exp_si(res, res, -2);
arb_exp(res, res, prec);
arb_mul(res, res, lam, prec);
arb_clear(pi);
arb_clear(lam);
}
void
acb_dirichlet_platt_scaled_lambda_vec(arb_ptr res,
const fmpz_t T, slong A, slong B, slong prec)
{
slong N = A*B;
if (A < 1 || B < 1 || N % 2)
{
flint_printf("requires an even number of grid points\n");
flint_abort();
}
else
{
slong i;
arb_t t;
arb_init(t);
for (i = 0; i < N; i++)
{
slong n = i - N/2;
_arb_div_si_si(t, n, A, prec);
arb_add_fmpz(t, t, T, prec);
acb_dirichlet_platt_scaled_lambda(res + i, t, prec);
}
arb_clear(t);
}
}
static void
_platt_beta(arb_t res, const arb_t t, slong prec)
{
arb_t u, v;
arb_init(u);
arb_init(v);
arb_log(u, t, prec);
arb_log(v, u, prec);
arb_div(u, v, u, prec);
_arb_inv_ui(res, 6, prec);
arb_add(res, res, u, prec);
arb_clear(u);
arb_clear(v);
}
static void
_platt_bound_C3_X(arb_t res, const arb_t t0, slong A, const arb_t H,
slong Ns, const arb_t beta, slong prec)
{
arb_t a, c, x;
arb_init(a);
arb_init(c);
arb_init(x);
/* res = gaussian(a=(t0+Ns/A)^beta, b=NULL, c=A*H, x=Ns) */
_arb_div_si_si(a, Ns, A, prec);
arb_add(a, a, t0, prec);
arb_pow(a, a, beta, prec);
arb_mul_si(c, H, A, prec);
arb_set_si(x, Ns);
_arb_gaussian(res, a, NULL, c, x, prec);
arb_clear(a);
arb_clear(c);
arb_clear(x);
}
static void
_platt_bound_C3_Y(arb_t res, const arb_t t0, slong A, const arb_t H,
slong Ns, const arb_t beta, slong prec)
{
arb_t a, b, g, g1, g2;
arb_init(a);
arb_init(b);
arb_init(g);
arb_init(g1);
arb_init(g2);
/* a = 2^((2*beta - 1)/2) */
arb_mul_2exp_si(a, beta, 1);
arb_sub_ui(a, a, 1, prec);
arb_mul_2exp_si(a, a, -1);
_arb_ui_pow_arb(a, 2, a, prec);
/* b = t0^beta */
arb_pow(b, t0, beta, prec);
/* g = incomplete_gamma(1/2, Ns^2 / (2*A^2*H^2)) */
arb_set_d(g1, 0.5);
_arb_div_si_si(g2, Ns, A, prec);
arb_div(g2, g2, H, prec);
arb_sqr(g2, g2, prec);
arb_mul_2exp_si(g2, g2, -1);
arb_hypgeom_gamma_upper(g, g1, g2, 0, prec);
/* res = a*b*A*H*g */
arb_mul_si(res, H, A, prec);
arb_mul(res, res, a, prec);
arb_mul(res, res, b, prec);
arb_mul(res, res, g, prec);
arb_clear(a);
arb_clear(b);
arb_clear(g);
arb_clear(g1);
arb_clear(g2);
}
static void
_platt_bound_C3_Z(arb_t res, const arb_t t0, slong A, const arb_t H,
const arb_t beta, slong prec)
{
arb_t a, b, g, g1, g2;
arb_init(a);
arb_init(b);
arb_init(g);
arb_init(g1);
arb_init(g2);
/* a = 2^((3*beta - 1)/2) */
arb_mul_ui(a, beta, 3, prec);
arb_sub_ui(a, a, 1, prec);
arb_mul_2exp_si(a, a, -1);
_arb_ui_pow_arb(a, 2, a, prec);
/* b = H^(beta + 1) */
arb_add_ui(b, beta, 1, prec);
arb_pow(b, H, b, prec);
/* g = incomplete_gamma((beta+1)/2, t0^2 / (2*H^2)) */
arb_add_ui(g1, beta, 1, prec);
arb_mul_2exp_si(g1, g1, -1);
arb_div(g2, t0, H, prec);
arb_sqr(g2, g2, prec);
arb_mul_2exp_si(g2, g2, -1);
arb_hypgeom_gamma_upper(g, g1, g2, 0, prec);
/* res = a*b*A*g */
arb_mul_si(res, g, A, prec);
arb_mul(res, res, a, prec);
arb_mul(res, res, b, prec);
arb_clear(a);
arb_clear(b);
arb_clear(g);
arb_clear(g1);
arb_clear(g2);
}
/*
* An improvement of Lemma C.3 is used, which is tighter by a factor of A.
* https://github.com/fredrik-johansson/arb/pull/269#issuecomment-468488347
*/
void
acb_dirichlet_platt_bound_C3(arb_t res, const arb_t t0, slong A, const arb_t H,
slong Ns, slong prec)
{
arb_t pi, ee, beta, X, Y, Z, lhs, rhs;
arb_init(pi);
arb_init(ee);
arb_init(beta);
arb_init(X);
arb_init(Y);
arb_init(Z);
arb_init(lhs);
arb_init(rhs);
/* t0 > exp(exp(1) is required for valid error bounds. */
arb_const_e(ee, prec);
arb_exp(ee, ee, prec);
if (!arb_gt(t0, ee))
{
arb_zero_pm_one(res);
goto finish;
}
/* 0 < Ns <= t0*A is required for valid error bounds. */
arb_set_si(lhs, Ns);
arb_mul_si(rhs, t0, A, prec);
if (!arb_is_positive(lhs) || !arb_le(lhs, rhs))
{
arb_zero_pm_one(res);
goto finish;
}
/* res = (X + Y + Z) * 6 / (pi * Ns) */
arb_const_pi(pi, prec);
_platt_beta(beta, t0, prec);
_platt_bound_C3_X(X, t0, A, H, Ns, beta, prec);
_platt_bound_C3_Y(Y, t0, A, H, Ns, beta, prec);
_platt_bound_C3_Z(Z, t0, A, H, beta, prec);
arb_add(res, X, Y, prec);
arb_add(res, res, Z, prec);
arb_mul_ui(res, res, 6, prec);
arb_div(res, res, pi, prec);
arb_div_si(res, res, Ns, prec);
finish:
arb_clear(pi);
arb_clear(ee);
arb_clear(beta);
arb_clear(X);
arb_clear(Y);
arb_clear(Z);
arb_clear(lhs);
arb_clear(rhs);
}
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 dt0, dt, a, s, err, total;
slong i;
slong N = A*B;
arb_init(dt0);
arb_init(dt);
arb_init(a);
arb_init(s);
arb_init(err);
arb_init(total);
arb_sub_fmpz(dt0, 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(dt, n, A, prec);
arb_sub(a, dt, dt0, prec);
arb_mul_si(a, a, A, prec);
arb_sinc_pi(a, a, prec);
arb_mul(a, a, p + i, prec);
_arb_gaussian(s, a, dt0, H, dt, prec);
arb_add(total, total, s, 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(dt0);
arb_clear(dt);
arb_clear(a);
arb_clear(s);
arb_clear(err);
arb_clear(total);
}
void
acb_dirichlet_platt_ws_precomp_init(acb_dirichlet_platt_ws_precomp_t pre,
slong A, const arb_t H, slong sigma, slong prec)
{
acb_dirichlet_platt_c_precomp_init(&pre->pre_c, sigma, H, 0, prec);
acb_dirichlet_platt_i_precomp_init(&pre->pre_i, A, H, sigma, prec);
}
void
acb_dirichlet_platt_ws_precomp_clear(acb_dirichlet_platt_ws_precomp_t pre)
{
acb_dirichlet_platt_c_precomp_clear(&pre->pre_c);
acb_dirichlet_platt_i_precomp_clear(&pre->pre_i);
}
void acb_dirichlet_platt_ws_interpolation_precomp(arb_t res,
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)
{
slong N = A*B;
if (A < 1 || B < 1 || N % 2)
{
flint_printf("requires an even number of grid points\n");
flint_abort();
}
else
{
arb_t x, dt0, dt0A, total;
arf_t lower_f;
slong n, lower_n;
arb_init(x);
arb_init(dt0);
arb_init(dt0A);
arb_init(total);
arf_init(lower_f);
arb_sub_fmpz(dt0, t0, T, prec + fmpz_clog_ui(T, 2));
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);
/*
* More than one iteration is needed only when the set of
* supporting points for interpolation is uncertain.
*/
for (n = lower_n; n == lower_n || arb_contains_si(dt0A, n); n++)
{
slong nlow = N/2 + n + 1;
slong nhigh = N/2 - n - 1;
slong Ns = FLINT_MIN(Ns_max, FLINT_MIN(nlow, nhigh));
if (Ns < 1)
{
arb_zero_pm_inf(total);
}
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)
{
arb_set(total, x);
}
else
{
arb_union(total, total, x, prec);
}
}
}
arb_set(res, total);
arb_clear(x);
arb_clear(dt0);
arb_clear(dt0A);
arb_clear(total);
arf_clear(lower_f);
}
}
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)
{
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);
acb_dirichlet_platt_ws_precomp_clear(pre);
}

101
acb_dirichlet/test/t-platt_ws_interpolation.c Normal file
View file

@ -0,0 +1,101 @@
/*
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 "acb_dirichlet.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("platt_ws_interpolation....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 5 * arb_test_multiplier(); iter++)
{
ulong A, B, iter2;
int Tbits;
fmpz_t T;
arb_ptr vec;
slong vprec;
A = 1 + n_randint(state, 10);
B = 1 + n_randint(state, 10);
if (n_randint(state, 2))
A *= 2;
else
B *= 2;
fmpz_init(T);
Tbits = 5 + n_randint(state, 20);
fmpz_set_ui(T, n_randtest_bits(state, Tbits));
vprec = 2 + n_randint(state, 300);
vec = _arb_vec_init(A*B);
acb_dirichlet_platt_scaled_lambda_vec(vec, T, A, B, vprec);
for (iter2 = 0; iter2 < 50; iter2++)
{
ulong Ns_max, sigma;
arb_t expected, observed;
arb_t t0, H;
slong prec;
arb_init(t0);
arb_init(H);
arb_init(expected);
arb_init(observed);
prec = 2 + n_randint(state, 300);
Ns_max = 1 + n_randint(state, 100);
sigma = 1 + 2*(1 + n_randint(state, 100));
arb_set_si(t0, A*B*500 - n_randint(state, A*B*1000) - 1);
arb_div_ui(t0, t0, A*1000, prec);
arb_add_fmpz(t0, t0, T, prec);
arb_set_ui(H, 1 + n_randint(state, 10000));
arb_div_ui(H, H, 1000, prec);
arb_abs(H, H);
acb_dirichlet_platt_scaled_lambda(expected, t0, prec);
acb_dirichlet_platt_ws_interpolation(observed, t0, vec,
T, A, B, Ns_max, H, sigma, prec);
if (!arb_overlaps(expected, observed))
{
flint_printf("FAIL: overlap\n\n");
flint_printf("iter = %wd\n\n", iter);
flint_printf("iter2 = %wd\n\n", iter2);
flint_printf("A = %wu\n\n", A);
flint_printf("B = %wu\n\n", B);
flint_printf("T = "); fmpz_print(T); flint_printf("\n\n");
flint_printf("vprec = %wd\n\n", vprec);
flint_printf("Ns_max = %wu\n\n", Ns_max);
flint_printf("sigma = %wu\n\n", sigma);
flint_printf("prec = %wd\n\n", prec);
flint_printf("t0 = "); arb_printd(t0, 15); flint_printf("\n\n");
flint_printf("H = "); arb_printd(H, 15); flint_printf("\n\n");
}
arb_clear(t0);
arb_clear(H);
arb_clear(expected);
arb_clear(observed);
}
fmpz_clear(T);
_arb_vec_clear(vec, A*B);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

39
doc/source/acb_dirichlet.rst
View file

@ -733,3 +733,42 @@ mpmath [Joh2018b]_ by Juan Arias de Reyna, described in [Ari2012]_.
Compute `S(g_n)` where `g_n` is the *n*-th Gram point. Requires `n \ge -1`.
Riemann zeta function zeros (Platt's method)
-------------------------------------------------------------------------------
The following functions related to the Riemann zeta function use the ideas
and formulas described by David J. Platt in [Pla2017]_.
.. function:: void acb_dirichlet_platt_scaled_lambda(arb_t res, const arb_t t, slong prec)
Compute `\Lambda(t) e^{\pi t/4}` where
.. math ::
\Lambda(t) = \pi^{-\frac{it}{2}}
\Gamma\left(\frac{\frac{1}{2}+it}{2}\right)
\zeta\left(\frac{1}{2} + it\right)
is defined in the beginning of section 3 of [Pla2017]_. As explained in
[Pla2011]_ this function has the same zeros as `\zeta(1/2 + it)` and is
real-valued by the functional equation, and the exponential factor is
designed to counteract the decay of the gamma factor as `t` increases.
.. function:: void acb_dirichlet_platt_scaled_lambda_vec(arb_ptr res, const fmpz_t T, slong A, slong B, slong prec)
Compute :func:`acb_dirichlet_platt_scaled_lambda` at `N=AB` points on a
grid, following the notation of [Pla2017]_. The first point on the grid
is `T - B/2` and the distance between grid points is `1/A`. The product
`N=AB` must be an even integer.
.. 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)
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`.
*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
interpolation. *sigma* is an odd positive integer tuning parameter
`\sigma \in 2\mathbb{Z}_{>0}+1` used in computing error bounds.

4
doc/source/credits.rst
View file

@ -253,6 +253,10 @@ Bibliography
.. [Pet1999] \K. Petras, "On the computation of the Gauss-Legendre quadrature formula with a given precision", Journal of Computational and Applied Mathematics 112 (1999), 253-267
.. [Pla2011] \D. J. Platt, "Computing degree 1 L-functions rigorously", Ph.D. Thesis, University of Bristol (2011), https://people.maths.bris.ac.uk/~madjp/thesis5.pdf
.. [Pla2017] \D. J. Platt, "Isolating some non-trivial zeros of zeta", Mathematics of Computation 86 (2017), 2449-2467, https://doi.org/10.1090/mcom/3198
.. [PP2010] \K. H. Pilehrood and T. H. Pilehrood. "Series acceleration formulas for beta values", Discrete Mathematics and Theoretical Computer Science, DMTCS, 12 (2) (2010), 223-236, https://hal.inria.fr/hal-00990465/
.. [PS1973] \M. S. Paterson and L. J. Stockmeyer, "On the number of nonscalar multiplications necessary to evaluate polynomials", SIAM J. Comput (1973)