diff --git a/acb_dirichlet.h b/acb_dirichlet.h index ce1862ce..d79eb277 100644 --- a/acb_dirichlet.h +++ b/acb_dirichlet.h @@ -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); diff --git a/acb_dirichlet/platt_lemma_B1.c b/acb_dirichlet/platt_lemma_B1.c index 61f4b7d5..3373687c 100644 --- a/acb_dirichlet/platt_lemma_B1.c +++ b/acb_dirichlet/platt_lemma_B1.c @@ -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); } diff --git a/acb_dirichlet/platt_local_hardy_z_zeros.c b/acb_dirichlet/platt_local_hardy_z_zeros.c index c8bf6fde..f6930163 100644 --- a/acb_dirichlet/platt_local_hardy_z_zeros.c +++ b/acb_dirichlet/platt_local_hardy_z_zeros.c @@ -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; diff --git a/acb_dirichlet/platt_multieval.c b/acb_dirichlet/platt_multieval.c index 18b4628d..e6735055 100644 --- a/acb_dirichlet/platt_multieval.c +++ b/acb_dirichlet/platt_multieval.c @@ -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); } } diff --git a/acb_dirichlet/platt_multieval_threaded.c b/acb_dirichlet/platt_multieval_threaded.c index 38017739..4ae43d79 100644 --- a/acb_dirichlet/platt_multieval_threaded.c +++ b/acb_dirichlet/platt_multieval_threaded.c @@ -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); diff --git a/acb_dirichlet/test/t-platt_local_hardy_z_zeros.c b/acb_dirichlet/test/t-platt_local_hardy_z_zeros.c index 26c34eb0..67c6ee26 100644 --- a/acb_dirichlet/test/t-platt_local_hardy_z_zeros.c +++ b/acb_dirichlet/test/t-platt_local_hardy_z_zeros.c @@ -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); diff --git a/acb_dirichlet/test/t-platt_multieval.c b/acb_dirichlet/test/t-platt_multieval.c index 4277c32d..c8ff39ba 100644 --- a/acb_dirichlet/test/t-platt_multieval.c +++ b/acb_dirichlet/test/t-platt_multieval.c @@ -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); diff --git a/acb_dirichlet/test/t-platt_multieval_threaded.c b/acb_dirichlet/test/t-platt_multieval_threaded.c index f36a80d9..224310af 100644 --- a/acb_dirichlet/test/t-platt_multieval_threaded.c +++ b/acb_dirichlet/test/t-platt_multieval_threaded.c @@ -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); diff --git a/doc/source/acb_dirichlet.rst b/doc/source/acb_dirichlet.rst index e2ed5fb7..ecd5c60c 100644 --- a/doc/source/acb_dirichlet.rst +++ b/doc/source/acb_dirichlet.rst @@ -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) diff --git a/doc/source/examples.rst b/doc/source/examples.rst index 62d8d612..36035af9 100644 --- a/doc/source/examples.rst +++ b/doc/source/examples.rst @@ -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 ------------------------------------------------------------------------------- diff --git a/examples/zeta_zeros.c b/examples/zeta_zeros.c index aaed43ac..fafdeca4 100644 --- a/examples/zeta_zeros.c +++ b/examples/zeta_zeros.c @@ -1,137 +1,279 @@ -/* This file is public domain. Author: Fredrik Johansson. */ +/* This file is public domain. Author: D.H.J. Polymath. */ #include #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; } -