This commit is contained in:
Fredrik Johansson 2022-05-23 23:38:01 +02:00
commit fc657c38fa
11 changed files with 377 additions and 176 deletions

View file

@ -287,23 +287,25 @@ void acb_dirichlet_platt_lemma_A9(arb_t out, slong sigma, const arb_t t0,
void acb_dirichlet_platt_lemma_A11(arb_t out, const arb_t t0, const arb_t h,
slong B, slong prec);
void acb_dirichlet_platt_lemma_B1(arb_t out, slong sigma, const arb_t t0,
const arb_t h, slong J, slong prec);
const arb_t h, const fmpz_t J, slong prec);
void acb_dirichlet_platt_lemma_B2(arb_t out, slong K, const arb_t h,
const arb_t xi, slong prec);
/* Platt DFT grid evaluation of scaled Lambda */
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 A, slong B, const arb_t h, const fmpz_t J, slong K,
slong sigma, slong prec);
void acb_dirichlet_platt_multieval_threaded(arb_ptr out, const fmpz_t T,
slong A, slong B, const arb_t h, slong J, slong K, slong sigma, slong prec);
slong A, slong B, const arb_t h, const fmpz_t J, slong K,
slong sigma, slong prec);
/* Platt Hardy Z zeros */
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,
const arb_t h, const fmpz_t J, slong K, slong sigma_grid,
slong Ns_max, const arb_t H, slong sigma_interp, slong prec);
slong acb_dirichlet_platt_local_hardy_z_zeros(
arb_ptr res, const fmpz_t n, slong len, slong prec);

View file

@ -11,22 +11,11 @@
#include "acb_dirichlet.h"
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);
}
void
acb_dirichlet_platt_lemma_B1(arb_t out,
slong sigma, const arb_t t0, const arb_t h, slong J, slong prec)
slong sigma, const arb_t t0, const arb_t h, const fmpz_t J, slong prec)
{
arb_t pi, C;
arb_t x1, x2, x3;
arb_t pi, C, x1, x2, x3, Ja;
if (sigma % 2 == 0 || sigma < 3)
{
@ -39,9 +28,11 @@ acb_dirichlet_platt_lemma_B1(arb_t out,
arb_init(x1);
arb_init(x2);
arb_init(x3);
arb_init(Ja);
arb_const_pi(pi, prec);
acb_dirichlet_platt_c_bound(C, sigma, t0, h, 0, prec);
arb_set_fmpz(Ja, J);
arb_set_si(x1, 2*sigma - 1);
arb_div(x1, x1, h, prec);
@ -54,7 +45,7 @@ acb_dirichlet_platt_lemma_B1(arb_t out,
arb_pow(x2, pi, x2, prec);
arb_set_si(x3, 1 - sigma);
_arb_ui_pow_arb(x3, (ulong) J, x3, prec);
arb_pow(x3, Ja, x3, prec);
arb_div_si(x3, x3, sigma - 1, prec);
arb_mul(out, x1, x2, prec);
@ -66,4 +57,5 @@ acb_dirichlet_platt_lemma_B1(arb_t out,
arb_clear(x1);
arb_clear(x2);
arb_clear(x3);
arb_clear(Ja);
}

View file

@ -55,7 +55,7 @@ typedef const platt_ctx_struct * platt_ctx_srcptr;
static void
platt_ctx_init(platt_ctx_t ctx,
const fmpz_t T, slong A, slong B,
const arb_t h, slong J, slong K, slong sigma_grid,
const arb_t h, const fmpz_t J, slong K, slong sigma_grid,
slong Ns_max, const arb_t H, slong sigma_interp, slong prec)
{
fmpz_init(&ctx->T);
@ -1097,7 +1097,7 @@ slong
_acb_dirichlet_platt_isolate_local_hardy_z_zeros(
arf_interval_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,
const arb_t h, const fmpz_t J, slong K, slong sigma_grid,
slong Ns_max, const arb_t H, slong sigma_interp, slong prec)
{
slong zeros_count;
@ -1296,7 +1296,7 @@ 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,
const arb_t h, const fmpz_t J, slong K, slong sigma_grid,
slong Ns_max, const arb_t H, slong sigma_interp, slong prec)
{
slong zeros_count, i;
@ -1345,9 +1345,9 @@ static platt_ctx_ptr
_create_heuristic_context(const fmpz_t n, slong prec)
{
platt_ctx_ptr p = NULL;
slong J, K, A, B, Ns_max, sigma_grid, sigma_interp;
slong K, A, B, Ns_max, sigma_grid, sigma_interp;
slong kbits;
fmpz_t T, k;
fmpz_t J, T, k;
arb_t g, h, H, logT;
double dlogJ, dK, dgrid, dh, dH, dinterp;
double x, x2, x3, x4;
@ -1436,7 +1436,7 @@ _create_heuristic_context(const fmpz_t n, slong prec)
arb_set_d(h, dh);
arb_set_d(H, dH);
J = (slong) exp(dlogJ);
fmpz_set_si(J, (slong) exp(dlogJ));
K = (slong) dK;
sigma_grid = ((slong) (dgrid/2))*2 + 1;
sigma_interp = ((slong) (dinterp/2))*2 + 1;

View file

@ -167,10 +167,10 @@ platt_g_table(acb_ptr table, slong A, slong B,
}
static void
logjsqrtpi(arb_t out, slong j, slong prec)
logjsqrtpi(arb_t out, const fmpz_t j, slong prec)
{
arb_const_sqrt_pi(out, prec);
arb_mul_si(out, out, j, prec);
arb_mul_fmpz(out, out, j, prec);
arb_log(out, out, prec);
}
@ -209,15 +209,17 @@ _acb_vec_scalar_add_error_arb_mag(acb_ptr res, slong len, const arb_t x)
* log(j*sqrt(pi))/(2*pi) >= m/B - 1/(2*B)
*/
void
get_smk_points(slong * res, slong A, slong B)
get_smk_points(fmpz * res, slong A, slong B)
{
slong m, N, prec;
arb_t x, u, v;
fmpz_t z;
arb_init(x);
arb_init(u); /* pi / B */
arb_init(v); /* 1 / sqrt(pi) */
fmpz_init(z);
N = A*B;
prec = 4;
arb_indeterminate(u);
@ -233,7 +235,7 @@ get_smk_points(slong * res, slong A, slong B)
arb_ceil(x, x, prec);
if (arb_get_unique_fmpz(z, x))
{
res[m] = fmpz_get_si(z);
fmpz_set(res + m, z);
break;
}
else
@ -253,7 +255,7 @@ get_smk_points(slong * res, slong A, slong B)
}
slong
platt_get_smk_index(slong B, slong j, slong prec)
platt_get_smk_index(slong B, const fmpz_t j, slong prec)
{
slong m;
arb_t pi, x;
@ -356,11 +358,12 @@ smk_block_accumulate(smk_block_t p, acb_ptr res, slong prec)
void
_platt_smk(acb_ptr table, acb_ptr startvec, acb_ptr stopvec,
const slong * smk_points, const arb_t t0, slong A, slong B,
slong jstart, slong jstop, slong mstart, slong mstop,
const fmpz * smk_points, const arb_t t0, slong A, slong B,
const fmpz_t jstart, const fmpz_t jstop, slong mstart, slong mstop,
slong K, slong prec)
{
slong j, k, m;
fmpz_t j, j_plus_one;
slong k, m;
slong N = A * B;
smk_block_t block;
acb_ptr accum;
@ -368,6 +371,8 @@ _platt_smk(acb_ptr table, acb_ptr startvec, acb_ptr stopvec,
arb_t rpi, logsqrtpi, rsqrtj, um, a, base;
acb_t z;
fmpz_init(j);
fmpz_init(j_plus_one);
arb_init(rpi);
arb_init(logsqrtpi);
arb_init(rsqrtj);
@ -387,13 +392,15 @@ _platt_smk(acb_ptr table, acb_ptr startvec, acb_ptr stopvec,
m = platt_get_smk_index(B, jstart, prec);
_arb_div_si_si(um, m, B, prec);
for (j = jstart; j <= jstop; j++)
for (fmpz_set(j, jstart); fmpz_cmp(j, jstop) <= 0; fmpz_add_ui(j, j, 1))
{
arb_log_ui(a, (ulong) j, prec);
arb_log_fmpz(a, j, prec);
arb_add(a, a, logsqrtpi, prec);
arb_mul(a, a, rpi, prec);
arb_rsqrt_ui(rsqrtj, (ulong) j, prec);
/* todo arb_rsqrt_fmpz */
arb_sqrt_fmpz(rsqrtj, j, prec);
arb_inv(rsqrtj, rsqrtj, prec);
acb_set_arb(z, t0);
acb_mul_arb(z, z, a, prec);
@ -401,7 +408,7 @@ _platt_smk(acb_ptr table, acb_ptr startvec, acb_ptr stopvec,
acb_exp_pi_i(z, z, prec);
acb_mul_arb(z, z, rsqrtj, prec);
while (m < N - 1 && smk_points[m + 1] <= j)
while (m < N - 1 && fmpz_cmp(smk_points + m + 1, j) <= 0)
{
m += 1;
_arb_div_si_si(um, m, B, prec);
@ -421,8 +428,12 @@ _platt_smk(acb_ptr table, acb_ptr startvec, acb_ptr stopvec,
smk_block_increment(block, z, diff_powers);
{
int j_stops = j == jstop;
int m_increases = m < N - 1 && smk_points[m + 1] <= j + 1;
int j_stops, m_increases;
fmpz_add_ui(j_plus_one, j, 1);
j_stops = fmpz_equal(j, jstop);
m_increases = (m < N - 1 &&
fmpz_cmp(smk_points + m + 1, j_plus_one) <= 0);
if (j_stops || m_increases || smk_block_is_full(block))
{
smk_block_accumulate(block, accum, prec);
@ -448,6 +459,8 @@ _platt_smk(acb_ptr table, acb_ptr startvec, acb_ptr stopvec,
}
}
fmpz_clear(j);
fmpz_clear(j_plus_one);
arb_clear(rpi);
arb_clear(logsqrtpi);
arb_clear(rsqrtj);
@ -539,7 +552,7 @@ remove_gaussian_window(arb_ptr out, slong A, slong B, const arb_t h, slong prec)
void
_acb_dirichlet_platt_multieval(arb_ptr out, acb_srcptr S_table,
const arb_t t0, slong A, slong B, const arb_t h, slong J,
const arb_t t0, slong A, slong B, const arb_t h, const fmpz_t J,
slong K, slong sigma, slong prec)
{
slong N = A*B;
@ -612,7 +625,7 @@ _acb_dirichlet_platt_multieval(arb_ptr out, acb_srcptr S_table,
acb_dirichlet_platt_lemma_B1(err, sigma, t0, h, J, prec);
_acb_vec_scalar_add_error_arb_mag(out_a, N/2 + 1, err);
arb_sqrt_ui(c, (ulong) J, prec);
arb_sqrt_fmpz(c, J, prec);
arb_mul_2exp_si(c, c, 1);
arb_sub_ui(c, c, 1, prec);
acb_dirichlet_platt_lemma_B2(err, K, h, xi, prec);
@ -660,7 +673,7 @@ _acb_dirichlet_platt_multieval(arb_ptr out, acb_srcptr S_table,
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)
const arb_t h, const fmpz_t J, slong K, slong sigma, slong prec)
{
if (flint_get_num_threads() > 1)
{
@ -672,23 +685,29 @@ acb_dirichlet_platt_multieval(arb_ptr out, const fmpz_t T, slong A, slong B,
slong N = A*B;
acb_ptr S;
arb_t t0;
slong * smk_points;
fmpz_t one;
fmpz * smk_points;
smk_points = flint_malloc(N * sizeof(slong));
smk_points = _fmpz_vec_init(N);
get_smk_points(smk_points, A, B);
fmpz_init(one);
fmpz_one(one);
arb_init(t0);
S = _acb_vec_init(K*N);
arb_set_fmpz(t0, T);
_platt_smk(S, NULL, NULL, smk_points, t0, A, B, 1, J, 0, N-1, K, prec);
_platt_smk(S, NULL, NULL, smk_points,
t0, A, B, one, J, 0, N-1, K, prec);
_acb_dirichlet_platt_multieval(out, S, t0, A, B, h, J, K, sigma, prec);
arb_clear(t0);
fmpz_clear(one);
_acb_vec_clear(S, K*N);
flint_free(smk_points);
_fmpz_vec_clear(smk_points, N);
}
}

View file

@ -13,16 +13,16 @@
#include "acb_dirichlet.h"
#include "pthread.h"
slong platt_get_smk_index(slong B, slong j, slong prec);
slong platt_get_smk_index(slong B, const fmpz_t j, slong prec);
void get_smk_points(slong * res, slong A, slong B);
void _platt_smk(acb_ptr table, acb_ptr startvec, acb_ptr stopvec,
const slong * smk_points, const arb_t t0, slong A, slong B,
slong jstart, slong jstop, slong mstart, slong mstop,
const fmpz_t jstart, const fmpz_t jstop, slong mstart, slong mstop,
slong K, slong prec);
void _acb_dirichlet_platt_multieval(arb_ptr out, acb_srcptr S_table,
const arb_t t0, slong A, slong B, const arb_t h, slong J,
const arb_t t0, slong A, slong B, const arb_t h, const fmpz_t J,
slong K, slong sigma, slong prec);
typedef struct
@ -35,8 +35,8 @@ typedef struct
slong A;
slong B;
slong K;
slong jstart;
slong jstop;
fmpz_t jstart;
fmpz_t jstop;
slong mstart;
slong mstop;
slong prec;
@ -56,21 +56,31 @@ _platt_smk_thread(void * arg_ptr)
void
acb_dirichlet_platt_multieval_threaded(arb_ptr out, const fmpz_t T, slong A,
slong B, const arb_t h, slong J, slong K, slong sigma, slong prec)
slong B, const arb_t h, const fmpz_t J, slong K,
slong sigma, slong prec)
{
slong i, num_threads, N, threadtasks;
slong * smk_points;
slong i, num_threads, N;
fmpz * smk_points;
pthread_t * threads;
platt_smk_arg_t * args;
acb_ptr S;
arb_t t0;
fmpz_t threadtasks;
num_threads = flint_get_num_threads();
if (num_threads < 1)
{
flint_printf("no threads available\n");
flint_abort();
}
N = A*B;
num_threads = flint_get_num_threads();
fmpz_init(threadtasks);
threads = flint_malloc(sizeof(pthread_t) * num_threads);
args = flint_malloc(sizeof(platt_smk_arg_t) * num_threads);
threadtasks = (J+num_threads-1)/num_threads;
smk_points = flint_malloc(N * sizeof(slong));
fmpz_add_si(threadtasks, J, num_threads - 1);
fmpz_tdiv_q_ui(threadtasks, threadtasks, (ulong) num_threads);
smk_points = _fmpz_vec_init(N);
arb_init(t0);
get_smk_points(smk_points, A, B);
@ -88,12 +98,15 @@ acb_dirichlet_platt_multieval_threaded(arb_ptr out, const fmpz_t T, slong A,
args[i].B = B;
args[i].K = K;
args[i].prec = prec;
args[i].jstart = i*threadtasks + 1;
args[i].jstop = (i+1)*threadtasks;
fmpz_init(args[i].jstart);
fmpz_init(args[i].jstop);
fmpz_mul_si(args[i].jstart, threadtasks, i);
fmpz_add_ui(args[i].jstart, args[i].jstart, 1);
fmpz_mul_si(args[i].jstop, threadtasks, i + 1);
args[i].mstart = platt_get_smk_index(B, args[i].jstart, prec);
args[i].mstop = platt_get_smk_index(B, args[i].jstop, prec);
}
args[num_threads-1].jstop = J;
fmpz_set(args[num_threads-1].jstop, J);
args[num_threads-1].mstop = platt_get_smk_index(B, J, prec);
for (i = 0; i < num_threads; i++)
@ -119,13 +132,15 @@ acb_dirichlet_platt_multieval_threaded(arb_ptr out, const fmpz_t T, slong A,
}
_acb_vec_clear(args[i].startvec, K);
_acb_vec_clear(args[i].stopvec, K);
fmpz_clear(args[i].jstart);
fmpz_clear(args[i].jstop);
}
_acb_dirichlet_platt_multieval(out, S, t0, A, B, h, J, K, sigma, prec);
arb_clear(t0);
_acb_vec_clear(S, K*N);
flint_free(smk_points);
_fmpz_vec_clear(smk_points, N);
flint_free(args);
flint_free(threads);

View file

@ -15,9 +15,9 @@ 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;
slong A, B, K, sigma_grid, Ns_max, sigma_interp;
arb_t h, H;
fmpz_t T, n;
fmpz_t J, T, n;
arb_ptr pa, pb;
slong count, i;
slong maxcount = 50;
@ -28,6 +28,7 @@ int main()
arb_init(h);
arb_init(H);
fmpz_init(J);
fmpz_init(T);
fmpz_init(n);
pa = _arb_vec_init(maxcount);
@ -41,7 +42,7 @@ int main()
B = 128;
/* tuning parameters for the evaluation of grid points */
J = 1000;
fmpz_set_si(J, 1000);
K = 30;
sigma_grid = 63;
arb_set_d(h, 4.5);
@ -76,6 +77,7 @@ int main()
arb_clear(h);
arb_clear(H);
fmpz_clear(J);
fmpz_clear(T);
fmpz_clear(n);
_arb_vec_clear(pa, maxcount);

View file

@ -74,17 +74,18 @@ int main()
slong A = 8;
slong B = 128;
slong N = A*B;
slong J = 1000;
slong K = 30;
slong sigma = 63;
slong prec = 128;
fmpz_t T;
fmpz_t J, T;
arb_t h;
arb_ptr vec;
arb_init(h);
fmpz_init(J);
fmpz_init(T);
fmpz_set_si(J, 1000);
fmpz_set_si(T, 10000);
arb_set_d(h, 4.5);
@ -172,6 +173,7 @@ int main()
arb_clear(r);
}
fmpz_clear(J);
fmpz_clear(T);
arb_clear(h);
_arb_vec_clear(vec, N);
@ -180,16 +182,15 @@ int main()
for (iter = 0; iter < 10 * arb_test_multiplier(); iter++)
{
slong prec;
ulong A, B, N, J, K;
ulong A, B, N, K;
slong sigma, Tbits;
fmpz_t T;
fmpz_t J, T;
arb_t h;
arb_ptr v1, v2;
/* better but slower limits are in parentheses below */
prec = 2 + n_randint(state, 300);
sigma = 1 + 2*(1 + n_randint(state, 100)); /* (200) */
J = 1 + n_randint(state, 100); /* (10000) */
K = 1 + n_randint(state, 20); /* (50) */
A = 1 + n_randint(state, 10);
B = 1 + n_randint(state, 10); /* (500) */
@ -199,7 +200,9 @@ int main()
B *= 2;
N = A*B;
fmpz_init(J);
fmpz_init(T);
fmpz_set_si(J, 1 + n_randint(state, 100)); /* (10000) */
Tbits = 5 + n_randint(state, 15);
fmpz_set_ui(T, n_randtest_bits(state, Tbits));
@ -217,13 +220,15 @@ int main()
flint_printf("FAIL: overlap\n\n");
flint_printf("iter = %wd prec = %wd\n\n", iter, prec);
flint_printf("sigma = %wd\n\n", sigma);
flint_printf("A = %wu B = %wu J = %wu K = %wu\n\n", A, B, J, K);
flint_printf("A = %wu B = %wu K = %wu\n\n", A, B, K);
flint_printf("J = "); fmpz_print(J); flint_printf("\n\n");
flint_printf("T = "); fmpz_print(T); flint_printf("\n\n");
flint_printf("h = "); arb_printn(h, 30, 0); flint_printf("\n\n");
flint_abort();
}
arb_clear(h);
fmpz_clear(J);
fmpz_clear(T);
_arb_vec_clear(v1, N);
_arb_vec_clear(v2, N);

View file

@ -47,17 +47,19 @@ int main()
slong A = 8;
slong B = 128;
slong N = A*B;
slong J = 1000;
slong K = 30;
slong sigma = 63;
slong prec = 128;
fmpz_t J;
fmpz_t T;
arb_t h;
arb_ptr vec;
arb_init(h);
fmpz_init(J);
fmpz_init(T);
fmpz_set_si(J, 1000);
fmpz_set_si(T, 10000);
arb_set_d(h, 4.5);
@ -88,6 +90,7 @@ int main()
arb_clear(r);
}
fmpz_clear(J);
fmpz_clear(T);
arb_clear(h);
_arb_vec_clear(vec, N);
@ -96,16 +99,15 @@ int main()
for (iter = 0; iter < 10 * arb_test_multiplier(); iter++)
{
slong prec;
ulong A, B, N, J, K;
ulong A, B, N, K;
slong sigma, Tbits;
fmpz_t T;
fmpz_t J, T;
arb_t h;
arb_ptr v1, v2;
/* better but slower limits are in parentheses below */
prec = 2 + n_randint(state, 300);
sigma = 1 + 2*(1 + n_randint(state, 100)); /* (200) */
J = 1 + n_randint(state, 100); /* (10000) */
K = 1 + n_randint(state, 20); /* (50) */
A = 1 + n_randint(state, 10);
B = 1 + n_randint(state, 10); /* (500) */
@ -115,7 +117,9 @@ int main()
B *= 2;
N = A*B;
fmpz_init(J);
fmpz_init(T);
fmpz_set_si(J, 1 + n_randint(state, 100)); /* (10000) */
Tbits = 5 + n_randint(state, 15);
fmpz_set_ui(T, n_randtest_bits(state, Tbits));
@ -135,13 +139,15 @@ int main()
flint_printf("FAIL: overlap\n\n");
flint_printf("iter = %wd prec = %wd\n\n", iter, prec);
flint_printf("sigma = %wd\n\n", sigma);
flint_printf("A = %wu B = %wu J = %wu K = %wu\n\n", A, B, J, K);
flint_printf("A = %wu B = %wu K = %wu\n\n", A, B, K);
flint_printf("J = "); fmpz_print(J); flint_printf("\n\n");
flint_printf("T = "); fmpz_print(T); flint_printf("\n\n");
flint_printf("h = "); arb_printn(h, 30, 0); flint_printf("\n\n");
flint_abort();
}
arb_clear(h);
fmpz_clear(J);
fmpz_clear(T);
_arb_vec_clear(v1, N);
_arb_vec_clear(v2, N);

View file

@ -791,9 +791,9 @@ and formulas described by David J. Platt in [Pla2017]_.
.. function:: void acb_dirichlet_platt_scaled_lambda_vec(arb_ptr res, const fmpz_t T, slong A, slong B, slong prec)
.. function:: void acb_dirichlet_platt_multieval(arb_ptr res, const fmpz_t T, slong A, slong B, const arb_t h, slong J, slong K, slong sigma, slong prec)
.. function:: void acb_dirichlet_platt_multieval(arb_ptr res, const fmpz_t T, slong A, slong B, const arb_t h, const fmpz_t J, slong K, slong sigma, slong prec)
.. function:: void acb_dirichlet_platt_multieval_threaded(arb_ptr res, const fmpz_t T, slong A, slong B, const arb_t h, slong J, slong K, slong sigma, slong prec)
.. function:: void acb_dirichlet_platt_multieval_threaded(arb_ptr res, const fmpz_t T, slong A, slong B, const arb_t h, const fmpz_t J, slong K, slong sigma, 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
@ -818,7 +818,7 @@ 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)
.. 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, const fmpz_t J, slong K, slong sigma_grid, slong Ns_max, const arb_t H, slong sigma_interp, slong prec)
.. function:: slong acb_dirichlet_platt_local_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, slong prec)
.. function:: slong acb_dirichlet_platt_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, slong prec)

View file

@ -526,6 +526,24 @@ squarefree factorization::
-1.0000
cpu/wall(s): 0 0.001
zeta_zeros.c
-------------------------------------------------------------------------------
This program finds the imaginary parts of consecutive nontrivial zeros
of the Riemann zeta function by calling either
:func:`acb_dirichlet_hardy_z_zeros` or
:func:`acb_dirichlet_platt_local_hardy_z_zeros` depending on the height
of the zeros and the number of zeros requested.
The program takes the following arguments::
zeta_zeros [-n n] [-count n] [-prec n] [-threads n] [-platt] [-noplatt] [-v] [-verbose] [-h] [-help]
> build/examples/zeta_zeros -n 1048449114 -count 2
1048449114 [388858886.0022851217767970582 +/- 7.46e-20]
1048449115 [388858886.0023936897027167201 +/- 7.59e-20]
cpu/wall(s): 0.255 0.255
virt/peak/res/peak(MB): 26.77 26.77 7.88 7.88
complex_plot.c
-------------------------------------------------------------------------------

View file

@ -1,137 +1,279 @@
/* This file is public domain. Author: Fredrik Johansson. */
/* This file is public domain. Author: D.H.J. Polymath. */
#include <string.h>
#include "acb_dirichlet.h"
#include "flint/profiler.h"
void print_zeros(arb_srcptr p, const fmpz_t n, slong len, slong digits)
{
slong i;
fmpz_t k;
fmpz_init_set(k, n);
for (i = 0; i < len; i++)
{
fmpz_print(k);
flint_printf("\t");
arb_printn(p+i, digits, 0);
flint_printf("\n");
fmpz_add_ui(k, k, 1);
}
fmpz_clear(k);
}
void print_help()
{
flint_printf("zeta_zeros [-n n] [-count n] [-prec n] [-threads n] "
"[-platt] [-noplatt] [-v] [-verbose] [-h] [-help]\n\n");
flint_printf("Reports the imaginary parts of consecutive nontrivial zeros "
"of the Riemann zeta function.\n");
}
void requires_value(int argc, char *argv[], slong i)
{
if (i == argc-1)
{
flint_printf("the argument %s requires a value\n", argv[i]);
flint_abort();
}
}
void invalid_value(char *argv[], slong i)
{
flint_printf("invalid value for the argument %s: %s\n", argv[i], argv[i+1]);
flint_abort();
}
int main(int argc, char *argv[])
{
arb_ptr res;
slong i, prec, num, digits, num_threads, found, alloc;
int platt, suppress, stream;
fmpz_t start;
if (argc < 2)
{
flint_printf("usage: zeta_zeros [-start n] [-num n] [-threads t] [-digits d] [-suppress] [-stream] [-platt]\n");
flint_printf("compute nontrivial zeros of the Riemann zeta function\n");
return 1;
}
fmpz_init(start);
fmpz_one(start);
digits = 30;
num_threads = 1;
suppress = 0;
num = 1;
platt = 0;
stream = 0;
const slong max_buffersize = 30000;
int verbose = 0;
int platt = 0;
int noplatt = 0;
slong i, buffersize, prec;
fmpz_t requested, count, nstart, n;
arb_ptr p;
for (i = 1; i < argc; i++)
{
if (!strcmp(argv[i], "-suppress"))
if (!strcmp(argv[i], "-h") || !strcmp(argv[i], "-help"))
{
suppress = 1;
print_help();
return 0;
}
if (!strcmp(argv[i], "-stream"))
}
fmpz_init(requested);
fmpz_init(count);
fmpz_init(nstart);
fmpz_init(n);
fmpz_one(nstart);
fmpz_set_si(requested, -1);
buffersize = max_buffersize;
prec = -1;
for (i = 1; i < argc; i++)
{
if (!strcmp(argv[i], "-noplatt"))
{
stream = 1;
}
else if (!strcmp(argv[i], "-num"))
{
num = atol(argv[i+1]);
i++;
}
else if (!strcmp(argv[i], "-threads"))
{
num_threads = atol(argv[i+1]);
i++;
}
else if (!strcmp(argv[i], "-start"))
{
fmpz_set_str(start, argv[i+1], 10);
i++;
}
else if (!strcmp(argv[i], "-digits"))
{
digits = atol(argv[i+1]);
i++;
noplatt = 1;
}
else if (!strcmp(argv[i], "-platt"))
{
platt = 1;
}
else if (!strcmp(argv[i], "-v") || !strcmp(argv[i], "-verbose"))
{
verbose = 1;
}
else if (!strcmp(argv[i], "-threads"))
{
slong threads;
requires_value(argc, argv, i);
threads = atol(argv[i+1]);
if (threads < 1)
{
invalid_value(argv, i);
}
flint_set_num_threads(threads);
}
else if (!strcmp(argv[i], "-n"))
{
requires_value(argc, argv, i);
if (fmpz_set_str(nstart, argv[i+1], 10) || fmpz_sgn(nstart) < 1)
{
invalid_value(argv, i);
}
}
else if (!strcmp(argv[i], "-count"))
{
requires_value(argc, argv, i);
if (fmpz_set_str(requested, argv[i+1], 10) ||
fmpz_sgn(requested) < 1)
{
invalid_value(argv, i);
}
if (fmpz_cmp_si(requested, buffersize) < 0)
{
buffersize = fmpz_get_si(requested);
}
}
else if (!strcmp(argv[i], "-prec"))
{
requires_value(argc, argv, i);
prec = atol(argv[i+1]);
if (prec < 2)
{
invalid_value(argv, i);
}
}
}
flint_set_num_threads(num_threads);
if (platt && noplatt)
{
flint_printf("conflicting arguments platt and noplatt\n");
flint_abort();
}
prec = digits * 3.3219280948873623 + 3;
if (platt && fmpz_cmp_si(nstart, 10000) < 0)
{
flint_printf("this implementation of the platt algorithm "
"is not valid below the 10000th zero\n");
flint_abort();
}
flint_printf("Computing %wd zeros starting with zero #", num);
fmpz_print(start);
flint_printf(", with prec = %wd bits...\n", prec);
fflush(stdout);
/* Above n~1e15 the large height method is better.
* Above n~1e11 the large height method is better for more than 100 zeros.
* Don't worry about crossing the threshold, just use the method
* that is better at the beginning of the run.
*/
if (!noplatt && !platt)
{
fmpz_t threshold;
fmpz_init(threshold);
fmpz_set_si(threshold, 10);
fmpz_pow_ui(threshold, threshold, 15);
if (fmpz_sgn(requested) < 0 || fmpz_cmp_si(requested, 100) > 0)
{
fmpz_set_si(threshold, 10);
fmpz_pow_ui(threshold, threshold, 11);
}
if (fmpz_cmp(nstart, threshold) < 0)
{
noplatt = 1;
}
else
{
platt = 1;
}
fmpz_clear(threshold);
}
if (prec == -1)
{
prec = 64 + fmpz_clog_ui(nstart, 2);
if (platt) prec *= 2;
}
if (verbose)
{
flint_printf("n: "); fmpz_print(nstart); flint_printf("\n");
flint_printf("count: "); fmpz_print(requested); flint_printf("\n");
flint_printf("threads: %wd\n", flint_get_num_threads());
flint_printf("prec: %wd\n", prec);
if (platt)
{
flint_printf("method: platt (good for large heights, "
"many consecutive zeros, and lower precision; "
"interprets prec as a working precision)\n");
}
else
{
flint_printf("method: noplatt (good for small heights, "
"few consecutive zeros, and greater precision; "
"interprets prec as a goal precision)\n");
}
}
p = _arb_vec_init(buffersize);
fmpz_set(n, nstart);
fmpz_zero(count);
TIMEIT_ONCE_START
/* This is the low height method. */
if (noplatt)
{
fmpz_t iter;
fmpz_init(iter);
while (fmpz_sgn(requested) < 0 || fmpz_cmp(count, requested) < 0)
{
slong num = buffersize;
if (fmpz_sgn(requested) >= 0)
{
fmpz_t remaining;
fmpz_init(remaining);
fmpz_sub(remaining, requested, count);
if (fmpz_cmp_si(remaining, num) < 0)
{
num = fmpz_get_si(remaining);
}
fmpz_clear(remaining);
}
if (fmpz_cmp_si(iter, 30) < 0)
{
num = FLINT_MIN(1 << fmpz_get_si(iter), num);
}
acb_dirichlet_hardy_z_zeros(p, n, num, prec);
print_zeros(p, n, num, prec*0.3 + 2);
fmpz_add_si(n, n, num);
fmpz_add_si(count, count, num);
fmpz_add_si(iter, iter, 1);
}
fmpz_clear(iter);
}
/* This is the large height method. */
if (platt)
{
alloc = num;
res = _arb_vec_init(alloc);
found = acb_dirichlet_platt_hardy_z_zeros(res, start, num, prec);
}
else if (stream)
{
fmpz_t n;
fmpz_init(n);
alloc = 1;
res = _arb_vec_init(alloc);
found = 0;
for (i = 0; i < num; i++)
while (fmpz_sgn(requested) < 0 || fmpz_cmp(count, requested) < 0)
{
fmpz_add_ui(n, start, i);
acb_dirichlet_hardy_z_zero(res, n, prec);
if (!suppress)
slong found;
slong num = buffersize;
if (fmpz_sgn(requested) >= 0)
{
arb_printn(res, digits, ARB_STR_NO_RADIUS);
flint_printf("\n");
fmpz_t remaining;
fmpz_init(remaining);
fmpz_sub(remaining, requested, count);
if (fmpz_cmp_si(remaining, num) < 0)
{
num = fmpz_get_si(remaining);
}
fmpz_clear(remaining);
}
found++;
found = acb_dirichlet_platt_local_hardy_z_zeros(p, n, num, prec);
if (!found)
{
flint_printf("Failed to find some zeros.\n");
flint_printf("Maybe prec is not high enough or something "
"is wrong with the internal tuning parameters.");
flint_abort();
}
print_zeros(p, n, found, prec*0.3 + 2);
fmpz_add_si(n, n, found);
fmpz_add_si(count, count, found);
}
fmpz_clear(n);
}
else
{
alloc = num;
res = _arb_vec_init(alloc);
acb_dirichlet_hardy_z_zeros(res, start, num, prec);
found = num;
}
TIMEIT_ONCE_STOP
SHOW_MEMORY_USAGE
if (found < num)
flint_printf("Found %wd zeros\n", found);
_arb_vec_clear(p, buffersize);
fmpz_clear(nstart);
fmpz_clear(n);
fmpz_clear(requested);
fmpz_clear(count);
if (!suppress && !stream)
{
for (i = 0; i < found; i++)
{
arb_printn(res + i, digits, ARB_STR_NO_RADIUS);
flint_printf("\n");
}
}
flint_printf("\n");
_arb_vec_clear(res, alloc);
fmpz_clear(start);
flint_cleanup_master();
flint_cleanup();
return 0;
}