mirror of
https://github.com/vale981/arb
synced 2025-03-05 17:31:38 -05:00
Merge branch 'master' of github.com:fredrik-johansson/arb
This commit is contained in:
fredrik
commit
7fa6176f4e
10 changed files with 1510 additions and 24 deletions
|
|
@ -252,7 +252,7 @@ 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,
|
||||
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,
|
||||
|
|
@ -288,6 +288,12 @@ void acb_dirichlet_platt_lemma_B2(arb_t out, slong K, const arb_t h,
|
|||
void acb_dirichlet_platt_multieval(arb_ptr out, const fmpz_t T, slong A,
|
||||
slong B, const arb_t h, slong J, slong K, slong sigma, slong prec);
|
||||
|
||||
slong _acb_dirichlet_platt_local_hardy_z_zeros(
|
||||
arb_ptr res, const fmpz_t n, slong len,
|
||||
const fmpz_t T, slong A, slong B,
|
||||
const arb_t h, slong J, slong K, slong sigma_grid,
|
||||
slong Ns_max, const arb_t H, slong sigma_interp, slong prec);
|
||||
|
||||
/* Discrete Fourier Transform */
|
||||
|
||||
void acb_dirichlet_dft_index(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec);
|
||||
|
|
|
|||
|
|
@ -1093,8 +1093,8 @@ _isolate_hardy_z_zeros(arf_interval_ptr res, const fmpz_t n, slong len)
|
|||
}
|
||||
|
||||
/* Isolate len zeros, starting from the nth zero. */
|
||||
static void
|
||||
isolate_hardy_z_zeros(arf_interval_ptr res, const fmpz_t n, slong len)
|
||||
void
|
||||
acb_dirichlet_isolate_hardy_z_zeros(arf_interval_ptr res, const fmpz_t n, slong len)
|
||||
{
|
||||
if (len <= 0)
|
||||
{
|
||||
|
|
@ -1392,7 +1392,7 @@ acb_dirichlet_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, slong prec)
|
|||
{
|
||||
slong i;
|
||||
arf_interval_ptr p = _arf_interval_vec_init(len);
|
||||
isolate_hardy_z_zeros(p, n, len);
|
||||
acb_dirichlet_isolate_hardy_z_zeros(p, n, len);
|
||||
for (i = 0; i < len; i++)
|
||||
{
|
||||
_acb_dirichlet_refine_hardy_z_zero(res + i, &p[i].a, &p[i].b, prec);
|
||||
|
|
|
|||
|
|
@ -28,7 +28,7 @@ _gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z,
|
|||
arb_init(x);
|
||||
for (i = 0; i < 5; i++)
|
||||
{
|
||||
arb_hypgeom_gamma_upper(x, s, z, regularized, prec * (1 << i));
|
||||
arb_hypgeom_gamma_upper(x, s, z, regularized, prec << i);
|
||||
if (arb_rel_accuracy_bits(x) > 1)
|
||||
{
|
||||
break;
|
||||
|
|
|
|||
|
|
@ -28,7 +28,7 @@ _gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z,
|
|||
arb_init(x);
|
||||
for (i = 0; i < 5; i++)
|
||||
{
|
||||
arb_hypgeom_gamma_upper(x, s, z, regularized, prec * (1 << i));
|
||||
arb_hypgeom_gamma_upper(x, s, z, regularized, prec << i);
|
||||
if (arb_rel_accuracy_bits(x) > 1)
|
||||
{
|
||||
break;
|
||||
|
|
|
|||
|
|
@ -29,7 +29,7 @@ _gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z,
|
|||
arb_init(x);
|
||||
for (i = 0; i < 5; i++)
|
||||
{
|
||||
arb_hypgeom_gamma_upper(x, s, z, regularized, prec * (1 << i));
|
||||
arb_hypgeom_gamma_upper(x, s, z, regularized, prec << i);
|
||||
if (arb_rel_accuracy_bits(x) > 1)
|
||||
{
|
||||
break;
|
||||
|
|
|
|||
1249
acb_dirichlet/platt_local_hardy_z_zeros.c
Normal file
1249
acb_dirichlet/platt_local_hardy_z_zeros.c
Normal file
File diff suppressed because it is too large
Load diff
|
|
@ -28,7 +28,7 @@ _gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z,
|
|||
arb_init(x);
|
||||
for (i = 0; i < 5; i++)
|
||||
{
|
||||
arb_hypgeom_gamma_upper(x, s, z, regularized, prec * (1 << i));
|
||||
arb_hypgeom_gamma_upper(x, s, z, regularized, prec << i);
|
||||
if (arb_rel_accuracy_bits(x) > 1)
|
||||
{
|
||||
break;
|
||||
|
|
@ -77,7 +77,9 @@ _arb_gaussian(arb_t res, const arb_t a, const arb_t b, const arb_t c,
|
|||
arb_mul_2exp_si(z, z, -1);
|
||||
arb_neg(z, z);
|
||||
arb_exp(z, z, prec);
|
||||
if (a != NULL)
|
||||
if (a == NULL)
|
||||
arb_set(res, z);
|
||||
else
|
||||
arb_mul(res, z, a, prec);
|
||||
arb_clear(z);
|
||||
}
|
||||
|
|
@ -339,41 +341,89 @@ _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;
|
||||
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;
|
||||
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(total);
|
||||
arb_init(pi);
|
||||
arb_init(g);
|
||||
arb_init(x);
|
||||
arb_init(c);
|
||||
|
||||
arb_const_pi(pi, prec);
|
||||
arb_sub_fmpz(dt0, t0, T, prec + fmpz_clog_ui(T, 2));
|
||||
|
||||
/* x = -N/2 - A*dt0 */
|
||||
arb_mul_si(x, dt0, A, prec);
|
||||
arb_add_si(x, x, N/2, prec);
|
||||
arb_neg(x, x);
|
||||
|
||||
/* c = sin(pi*x) / pi */
|
||||
arb_sin_pi(c, x, prec);
|
||||
arb_div(c, c, pi, prec);
|
||||
|
||||
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);
|
||||
_arb_gaussian(g, NULL, dt0, H, dt, prec);
|
||||
arb_mul(s, g, p + i, prec);
|
||||
arb_add_si(a, x, i, prec);
|
||||
arb_get_mag(m, a);
|
||||
if (mag_cmp_2exp_si(m, -1) < 0)
|
||||
{
|
||||
arb_sinc_pi(b, a, prec);
|
||||
arb_addmul(accum1, s, b, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
arb_div(b, s, a, prec);
|
||||
if (i % 2)
|
||||
{
|
||||
arb_neg(b, b);
|
||||
}
|
||||
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);
|
||||
|
||||
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(total);
|
||||
arb_clear(pi);
|
||||
arb_clear(g);
|
||||
arb_clear(x);
|
||||
arb_clear(c);
|
||||
}
|
||||
|
||||
|
||||
void
|
||||
acb_dirichlet_platt_ws_precomp_init(acb_dirichlet_platt_ws_precomp_t pre,
|
||||
slong A, const arb_t H, slong sigma, slong prec)
|
||||
|
|
@ -390,7 +440,7 @@ 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,
|
||||
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)
|
||||
{
|
||||
|
|
|
|||
87
acb_dirichlet/test/t-platt_local_hardy_z_zeros.c
Normal file
87
acb_dirichlet/test/t-platt_local_hardy_z_zeros.c
Normal file
|
|
@ -0,0 +1,87 @@
|
|||
/*
|
||||
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()
|
||||
{
|
||||
/* Check a specific combination of parameter values that is relatively fast
|
||||
* to evaluate and that has relatively tight bounds. */
|
||||
slong A, B, J, K, sigma_grid, Ns_max, sigma_interp;
|
||||
arb_t h, H;
|
||||
fmpz_t T, n;
|
||||
arb_ptr pa, pb;
|
||||
slong count, i;
|
||||
slong maxcount = 50;
|
||||
slong prec = 128;
|
||||
|
||||
flint_printf("platt_local_hardy_z_zeros....");
|
||||
fflush(stdout);
|
||||
|
||||
arb_init(h);
|
||||
arb_init(H);
|
||||
fmpz_init(T);
|
||||
fmpz_init(n);
|
||||
pa = _arb_vec_init(maxcount);
|
||||
pb = _arb_vec_init(maxcount);
|
||||
|
||||
fmpz_set_si(n, 10142);
|
||||
|
||||
/* parameters related to the location/resolution/width of the grid */
|
||||
fmpz_set_si(T, 10000);
|
||||
A = 8;
|
||||
B = 128;
|
||||
|
||||
/* tuning parameters for the evaluation of grid points */
|
||||
J = 1000;
|
||||
K = 30;
|
||||
sigma_grid = 63;
|
||||
arb_set_d(h, 4.5);
|
||||
|
||||
/* tuning parameters for interpolation on the grid */
|
||||
Ns_max = 200;
|
||||
sigma_interp = 21;
|
||||
arb_one(H);
|
||||
|
||||
count = _acb_dirichlet_platt_local_hardy_z_zeros(pa, n, maxcount,
|
||||
T, A, B, h, J, K, sigma_grid, Ns_max, H, sigma_interp, prec);
|
||||
acb_dirichlet_hardy_z_zeros(pb, n, count, prec);
|
||||
if (count != maxcount)
|
||||
{
|
||||
flint_printf("FAIL: not enough zeros were isolated\n\n");
|
||||
flint_printf("count = %wd maxcount = %wd\n\n", count, maxcount);
|
||||
flint_abort();
|
||||
}
|
||||
|
||||
for (i = 0; i < count; i++)
|
||||
{
|
||||
if (!arb_overlaps(pa+i, pb+i))
|
||||
{
|
||||
flint_printf("FAIL: overlap\n\n");
|
||||
flint_printf("observed[%wd] = ", i);
|
||||
arb_printd(pa+i, 20); flint_printf("\n\n");
|
||||
flint_printf("expected[%wd] = ", i);
|
||||
arb_printd(pb+i, 20); flint_printf("\n\n");
|
||||
flint_abort();
|
||||
}
|
||||
}
|
||||
|
||||
arb_clear(h);
|
||||
arb_clear(H);
|
||||
fmpz_clear(T);
|
||||
fmpz_clear(n);
|
||||
_arb_vec_clear(pa, maxcount);
|
||||
_arb_vec_clear(pb, maxcount);
|
||||
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
|
@ -11,6 +11,13 @@
|
|||
|
||||
#include "acb_dirichlet.h"
|
||||
|
||||
static void
|
||||
_arb_inv_si(arb_t res, slong a, slong prec)
|
||||
{
|
||||
arb_set_si(res, a);
|
||||
arb_inv(res, res, prec);
|
||||
}
|
||||
|
||||
static void
|
||||
_arb_div_si_si(arb_t res, slong a, slong b, slong prec)
|
||||
{
|
||||
|
|
@ -32,6 +39,26 @@ _arb_vec_overlaps(arb_srcptr a, arb_srcptr b, slong len)
|
|||
return 1;
|
||||
}
|
||||
|
||||
static void
|
||||
_check_containment(const char *name, const arb_t x, const char *s)
|
||||
{
|
||||
arb_t u;
|
||||
slong prec = 300;
|
||||
|
||||
arb_init(u);
|
||||
arb_set_str(u, s, prec);
|
||||
|
||||
if (!arb_contains(u, x))
|
||||
{
|
||||
flint_printf("FAIL: %s\n\n", name);
|
||||
flint_printf("observed = "); arb_printn(x, 30, 0); flint_printf("\n\n");
|
||||
flint_printf("expected = "); arb_printn(u, 30, 0); flint_printf("\n\n");
|
||||
flint_abort();
|
||||
}
|
||||
|
||||
arb_clear(u);
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
slong iter;
|
||||
|
|
@ -60,11 +87,69 @@ int main()
|
|||
|
||||
fmpz_set_si(T, 10000);
|
||||
arb_set_d(h, 4.5);
|
||||
|
||||
/* Spot-check lemma bound containment
|
||||
* in intervals calculated with PARI/GP. */
|
||||
{
|
||||
arb_t lem, xi, x, beta, t0;
|
||||
slong i = 201;
|
||||
slong k = 5;
|
||||
slong wp = 300;
|
||||
|
||||
arb_init(lem);
|
||||
arb_init(xi);
|
||||
arb_init(x);
|
||||
arb_init(t0);
|
||||
arb_init(beta);
|
||||
|
||||
_arb_inv_si(xi, B, wp);
|
||||
arb_mul_2exp_si(xi, xi, -1);
|
||||
_arb_div_si_si(x, i, B, wp);
|
||||
arb_set_fmpz(t0, T);
|
||||
acb_dirichlet_platt_beta(beta, t0, wp);
|
||||
|
||||
acb_dirichlet_platt_lemma_32(lem, h, t0, x, wp);
|
||||
_check_containment("Lemma 3.2", lem,
|
||||
"[5.3526496753240991744e-1072334 +/- 2.55e-1072354]");
|
||||
|
||||
acb_dirichlet_platt_c_bound(lem, sigma, t0, h, k, wp);
|
||||
_check_containment("Lemma A.3", lem,
|
||||
"[1.3516642396389823078e+134 +/- 2.65e+114]");
|
||||
|
||||
acb_dirichlet_platt_lemma_A5(lem, B, h, k, wp);
|
||||
_check_containment("Lemma A.5", lem,
|
||||
"[1.0075390047893384632e-30 +/- 5.57e-51]");
|
||||
|
||||
acb_dirichlet_platt_lemma_A7(lem, sigma, t0, h, k, A, wp);
|
||||
_check_containment("Lemma A.7", lem,
|
||||
"[3.0406705491484062400e-505 +/- 1.57e-525]");
|
||||
|
||||
acb_dirichlet_platt_lemma_A9(lem, sigma, t0, h, A, wp);
|
||||
_check_containment("Lemma A.9", lem,
|
||||
"[6.8953211848420326275e-536 +/- 3.52e-556]");
|
||||
|
||||
acb_dirichlet_platt_lemma_A11(lem, t0, h, B, wp);
|
||||
_check_containment("Lemma A.11", lem,
|
||||
"[3.0825745863006335768e-42 +/- 3.68e-62]");
|
||||
|
||||
acb_dirichlet_platt_lemma_B1(lem, sigma, t0, h, J, wp);
|
||||
_check_containment("Lemma B.1", lem,
|
||||
"[8.5737638613320328274e-42 +/- 7.50e-63]");
|
||||
|
||||
acb_dirichlet_platt_lemma_B2(lem, K, h, xi, wp);
|
||||
_check_containment("Lemma B.2", lem,
|
||||
"[2.0748437544358592615e-44 +/- 4.76e-64]");
|
||||
|
||||
arb_clear(lem);
|
||||
arb_clear(xi);
|
||||
arb_clear(x);
|
||||
arb_clear(t0);
|
||||
arb_clear(beta);
|
||||
}
|
||||
|
||||
/* Check a few random entries in the multieval vector. */
|
||||
vec = _arb_vec_init(N);
|
||||
|
||||
acb_dirichlet_platt_multieval(vec, T, A, B, h, J, K, sigma, prec);
|
||||
|
||||
/* Check only a few random entries in the multieval vector. */
|
||||
for (iter = 0; iter < 20; iter++)
|
||||
{
|
||||
arb_t t, r;
|
||||
|
|
@ -78,8 +163,8 @@ int main()
|
|||
if (!arb_overlaps(vec + i, r))
|
||||
{
|
||||
flint_printf("FAIL: overlap for hardcoded example\n\n");
|
||||
flint_printf("n = %wd\n\n", n);
|
||||
flint_printf("vec[i] = "); arb_printn(vec + i, 30, 0); flint_printf("\n\n");
|
||||
flint_printf("i = %wd n = %wd\n\n", i, n);
|
||||
flint_printf("vec[%wd] = ", i); arb_printn(vec + i, 30, 0); flint_printf("\n\n");
|
||||
flint_printf("r = "); arb_printn(r, 30, 0); flint_printf("\n\n");
|
||||
flint_abort();
|
||||
}
|
||||
|
|
|
|||
|
|
@ -777,3 +777,12 @@ and formulas described by David J. Platt in [Pla2017]_.
|
|||
interpolation. *sigma* is an odd positive integer tuning parameter
|
||||
`\sigma \in 2\mathbb{Z}_{>0}+1` used in computing error bounds.
|
||||
|
||||
.. function:: slong _acb_dirichlet_platt_local_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, const fmpz_t T, slong A, slong B, const arb_t h, slong J, slong K, slong sigma_grid, slong Ns_max, const arb_t H, slong sigma_interp, slong prec)
|
||||
|
||||
Sets the entries of *res* to at most *len* consecutive zeros of the
|
||||
Hardy Z-function, beginning with the *n*-th zero. The number of zeros
|
||||
isolated near *T* is returned. Requires positive *n*.
|
||||
Internally this function uses Platt's grid evaluation of the scaled
|
||||
Lambda function, and the final several parameters have the same meanings
|
||||
as in the functions :func:`acb_dirichlet_platt_multieval`
|
||||
and :func:`acb_dirichlet_platt_ws_interpolation`.
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue