diff --git a/acb_dirichlet.h b/acb_dirichlet.h index b2fd4020..ee52d55c 100644 --- a/acb_dirichlet.h +++ b/acb_dirichlet.h @@ -265,6 +265,29 @@ 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); +/* Platt lemma bounds of errors in the DFT grid evaluation of scaled Lambda */ + +void acb_dirichlet_platt_beta(arb_t res, const arb_t t, slong prec); +void acb_dirichlet_platt_lemma_32(arb_t out, const arb_t h, const arb_t t0, + const arb_t x, slong prec); +void acb_dirichlet_platt_lemma_A5(arb_t out, slong B, const arb_t h, slong k, + slong prec); +void acb_dirichlet_platt_lemma_A7(arb_t out, slong sigma, const arb_t t0, + const arb_t h, slong k, slong A, slong prec); +void acb_dirichlet_platt_lemma_A9(arb_t out, slong sigma, const arb_t t0, + const arb_t h, slong A, slong prec); +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); +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); + /* Discrete Fourier Transform */ void acb_dirichlet_dft_index(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec); diff --git a/acb_dirichlet/platt_beta.c b/acb_dirichlet/platt_beta.c new file mode 100644 index 00000000..6894f16e --- /dev/null +++ b/acb_dirichlet/platt_beta.c @@ -0,0 +1,34 @@ +/* + 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 . +*/ + +#include "acb_dirichlet.h" + +static void +_arb_inv_ui(arb_t res, ulong n, slong prec) +{ + arb_set_ui(res, n); + arb_inv(res, res, prec); +} + +void +acb_dirichlet_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); +} diff --git a/acb_dirichlet/platt_c_bound.c b/acb_dirichlet/platt_c_bound.c index bfea8874..107aca1b 100644 --- a/acb_dirichlet/platt_c_bound.c +++ b/acb_dirichlet/platt_c_bound.c @@ -12,6 +12,33 @@ #include "acb_dirichlet.h" #include "arb_hypgeom.h" +/* Increase precision adaptively. */ +static void +_gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z, + int regularized, slong prec) +{ + if (!arb_is_finite(s) || !arb_is_finite(z)) + { + arb_indeterminate(res); + } + else + { + arb_t x; + slong i; + arb_init(x); + for (i = 0; i < 5; i++) + { + arb_hypgeom_gamma_upper(x, s, z, regularized, prec * (1 << i)); + if (arb_rel_accuracy_bits(x) > 1) + { + break; + } + } + arb_swap(res, x); + arb_clear(x); + } +} + static void _arb_pow_si(arb_t res, const arb_t x, slong y, slong prec) { @@ -148,7 +175,7 @@ _pre_c_p(arb_ptr res, slong sigma, const arb_t h, ulong k, slong prec) arb_set_si(x, k + l + 1); arb_mul_2exp_si(x, x, -1); - arb_hypgeom_gamma_upper(x, x, u, 0, prec); + _gamma_upper_workaround(x, x, u, 0, prec); arb_mul(res + l, res + l, x, prec); } diff --git a/acb_dirichlet/platt_lemma_32.c b/acb_dirichlet/platt_lemma_32.c new file mode 100644 index 00000000..8578d679 --- /dev/null +++ b/acb_dirichlet/platt_lemma_32.c @@ -0,0 +1,50 @@ +/* + 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 . +*/ + +#include "acb_dirichlet.h" + +/* out = 2*pi^(5/4)*exp(1/(8*h^2) - (t0^2)/(2*h^2) - pi*x) */ +void +acb_dirichlet_platt_lemma_32(arb_t out, const arb_t h, const arb_t t0, + const arb_t x, slong prec) +{ + arb_t pi, one_fourth; + arb_t x1, x2; + + arb_init(pi); + arb_init(one_fourth); + arb_init(x1); + arb_init(x2); + + arb_const_pi(pi, prec); + arb_set_d(one_fourth, 0.25); + + arb_set_d(x1, 1.25); + arb_pow(x1, pi, x1, prec); + arb_mul_2exp_si(x1, x1, 1); + + arb_sqr(x2, t0, prec); + arb_sub(x2, x2, one_fourth, prec); + arb_div(x2, x2, h, prec); + arb_div(x2, x2, h, prec); + arb_mul_2exp_si(x2, x2, -1); + + arb_mul(out, pi, x, prec); + arb_add(out, out, x2, prec); + arb_neg(out, out); + arb_exp(out, out, prec); + arb_mul(out, out, x1, prec); + + arb_clear(pi); + arb_clear(one_fourth); + arb_clear(x1); + arb_clear(x2); +} diff --git a/acb_dirichlet/platt_lemma_A11.c b/acb_dirichlet/platt_lemma_A11.c new file mode 100644 index 00000000..18eaaa4b --- /dev/null +++ b/acb_dirichlet/platt_lemma_A11.c @@ -0,0 +1,227 @@ +/* + 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 . +*/ + +#include "acb_dirichlet.h" +#include "arb_hypgeom.h" + +/* Increase precision adaptively. */ +static void +_gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z, + int regularized, slong prec) +{ + if (!arb_is_finite(s) || !arb_is_finite(z)) + { + arb_indeterminate(res); + } + else + { + arb_t x; + slong i; + arb_init(x); + for (i = 0; i < 5; i++) + { + arb_hypgeom_gamma_upper(x, s, z, regularized, prec * (1 << i)); + if (arb_rel_accuracy_bits(x) > 1) + { + break; + } + } + arb_swap(res, x); + arb_clear(x); + } +} + +static void +_platt_lemma_A11_X(arb_t out, + const arb_t t0, const arb_t h, slong B, const arb_t beta, slong prec) +{ + arb_t x1, x2; + + arb_init(x1); + arb_init(x2); + + arb_set_si(x1, B); + arb_mul_2exp_si(x1, x1, -1); + arb_add(x1, x1, t0, prec); + arb_pow(x1, x1, beta, prec); + + arb_set_si(x2, B); + arb_div(x2, x2, h, prec); + arb_sqr(x2, x2, prec); + arb_mul_2exp_si(x2, x2, -3); + arb_neg(x2, x2); + arb_exp(x2, x2, prec); + + arb_mul(out, x1, x2, prec); + + arb_clear(x1); + arb_clear(x2); +} + +static void +_platt_lemma_A11_Y(arb_t out, + const arb_t t0, const arb_t h, slong B, const arb_t beta, slong prec) +{ + arb_t x1, x2, x3, x4, x5; + + arb_init(x1); + arb_init(x2); + arb_init(x3); + arb_init(x4); + arb_init(x5); + + arb_rsqrt_ui(x1, 2, prec); + + arb_pow(x2, t0, beta, prec); + + arb_one(x3); + arb_mul_2exp_si(x3, x3, -1); + + arb_set_si(x4, B); + arb_div(x4, x4, h, prec); + arb_sqr(x4, x4, prec); + arb_mul_2exp_si(x4, x4, -3); + + _gamma_upper_workaround(x5, x3, x4, 0, prec); + + arb_mul(out, x1, x2, prec); + arb_mul(out, out, x5, prec); + + arb_clear(x1); + arb_clear(x2); + arb_clear(x3); + arb_clear(x4); + arb_clear(x5); +} + +static void +_platt_lemma_A11_Z(arb_t out, + const arb_t t0, const arb_t h, const arb_t beta, slong prec) +{ + arb_t two; + arb_t x1, x2, x3, x4, x5; + + arb_init(two); + arb_init(x1); + arb_init(x2); + arb_init(x3); + arb_init(x4); + arb_init(x5); + + arb_set_ui(two, 2); + + arb_sub_ui(x1, beta, 1, prec); + arb_mul_2exp_si(x1, x1, -1); + arb_pow(x1, two, x1, prec); + + arb_pow(x2, h, beta, prec); + + arb_add_ui(x3, beta, 1, prec); + arb_mul_2exp_si(x3, x3, -1); + + arb_div(x4, t0, h, prec); + arb_sqr(x4, x4, prec); + arb_mul_2exp_si(x4, x4, -1); + + _gamma_upper_workaround(x5, x3, x4, 0, prec); + + arb_mul(out, x1, x2, prec); + arb_mul(out, out, x5, prec); + + arb_clear(two); + arb_clear(x1); + arb_clear(x2); + arb_clear(x3); + arb_clear(x4); + arb_clear(x5); +} + +static int +_platt_lemma_A11_constraint(const arb_t t0, const arb_t h, slong B, + const arb_t beta, slong prec) +{ + int result; + + arb_t a, b, expe; + arb_init(a); + arb_init(b); + arb_init(expe); + + /* expe = exp(e) */ + arb_const_e(expe, prec); + arb_exp(expe, expe, prec); + + /* a = beta*h^2 / t0 */ + arb_sqr(a, h, prec); + arb_mul(a, a, beta, prec); + arb_div(a, a, t0, prec); + + /* b = B/2 */ + arb_set_si(b, B); + arb_mul_2exp_si(b, b, -1); + + result = arb_le(a, b) && arb_le(b, t0) && arb_gt(t0, expe); + + arb_clear(a); + arb_clear(b); + arb_clear(expe); + + return result; +} + +void +acb_dirichlet_platt_lemma_A11(arb_t out, const arb_t t0, const arb_t h, + slong B, slong prec) +{ + arb_t beta; + arb_init(beta); + acb_dirichlet_platt_beta(beta, t0, prec); + + if (_platt_lemma_A11_constraint(t0, h, B, beta, prec)) + { + arb_t X, Y, Z; + arb_t x1, x2; + + arb_init(X); + arb_init(Y); + arb_init(Z); + arb_init(x1); + arb_init(x2); + + _platt_lemma_A11_X(X, t0, h, B, beta, prec); + _platt_lemma_A11_Y(Y, t0, h, B, beta, prec); + _platt_lemma_A11_Z(Z, t0, h, beta, prec); + + arb_set_ui(x1, 2); + arb_pow(x1, x1, beta, prec); + arb_mul(x1, x1, h, prec); + arb_div_si(x1, x1, B, prec); + + arb_add(x2, Y, Z, prec); + arb_mul(x2, x2, x1, prec); + arb_add(x2, x2, X, prec); + arb_mul_ui(x2, x2, 6, prec); + + arb_set(out, x2); + + arb_clear(X); + arb_clear(Y); + arb_clear(Z); + arb_clear(x1); + arb_clear(x2); + } + else + { + arb_zero_pm_inf(out); + } + + arb_clear(beta); +} diff --git a/acb_dirichlet/platt_lemma_A5.c b/acb_dirichlet/platt_lemma_A5.c new file mode 100644 index 00000000..81f515b9 --- /dev/null +++ b/acb_dirichlet/platt_lemma_A5.c @@ -0,0 +1,124 @@ +/* + 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 . +*/ + +#include "acb_dirichlet.h" +#include "arb_hypgeom.h" + + +/* Increase precision adaptively. */ +static void +_gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z, + int regularized, slong prec) +{ + if (!arb_is_finite(s) || !arb_is_finite(z)) + { + arb_indeterminate(res); + } + else + { + arb_t x; + slong i; + arb_init(x); + for (i = 0; i < 5; i++) + { + arb_hypgeom_gamma_upper(x, s, z, regularized, prec * (1 << i)); + if (arb_rel_accuracy_bits(x) > 1) + { + break; + } + } + arb_swap(res, x); + arb_clear(x); + } +} + +/* Lemma A5 requires B > h*sqrt(k) */ +static int +_platt_lemma_A5_constraint(slong B, const arb_t h, slong k, slong prec) +{ + int result; + arb_t lhs, rhs; + arb_init(rhs); + arb_init(lhs); + arb_set_si(lhs, B); + arb_sqrt_ui(rhs, (ulong) k, prec); + arb_mul(rhs, rhs, h, prec); + result = arb_gt(lhs, rhs); + arb_clear(rhs); + arb_clear(lhs); + return result; +} + + +void +acb_dirichlet_platt_lemma_A5(arb_t out, slong B, const arb_t h, + slong k, slong prec) +{ + arb_t a, b; + arb_t x1, x2, x3, x4, x5, x6; + + if (!_platt_lemma_A5_constraint(B, h, k, prec)) + { + arb_zero_pm_inf(out); + return; + } + + arb_init(a); + arb_init(b); + arb_init(x1); + arb_init(x2); + arb_init(x3); + arb_init(x4); + arb_init(x5); + arb_init(x6); + + arb_const_pi(x1, prec); + arb_mul_si(x1, x1, B, prec); + arb_pow_ui(x1, x1, (ulong) k, prec); + arb_mul_2exp_si(x1, x1, 3); + + arb_set_si(a, B); + arb_div(a, a, h, prec); + arb_sqr(a, a, prec); + arb_mul_2exp_si(a, a, -3); + + arb_neg(x2, a); + arb_exp(x2, x2, prec); + + arb_set_si(x3, 3*k - 1); + arb_mul_2exp_si(x3, x3, -1); + + arb_set_ui(x4, 2); + arb_pow(x4, x4, x3, prec); + + arb_set_si(b, k+1); + + arb_div_si(x5, h, B, prec); + arb_pow_ui(x5, x5, (ulong) (k+1), prec); + + arb_mul_2exp_si(x6, b, -1); + + _gamma_upper_workaround(x6, x6, a, 0, prec); + + arb_mul(out, x4, x5, prec); + arb_mul(out, out, x6, prec); + arb_add(out, out, x2, prec); + arb_mul(out, out, x1, prec); + + arb_clear(a); + arb_clear(b); + arb_clear(x1); + arb_clear(x2); + arb_clear(x3); + arb_clear(x4); + arb_clear(x5); + arb_clear(x6); +} diff --git a/acb_dirichlet/platt_lemma_A7.c b/acb_dirichlet/platt_lemma_A7.c new file mode 100644 index 00000000..26c24d79 --- /dev/null +++ b/acb_dirichlet/platt_lemma_A7.c @@ -0,0 +1,180 @@ +/* + 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 . +*/ + +#include "acb_dirichlet.h" + +static void +_platt_lemma_A7_S(arb_t out, slong sigma, + const arb_t t0, const arb_t h, slong k, slong A, slong prec) +{ + slong l; + arb_t total, summand; + arb_t pi, half; + arb_t a; + arb_t l_factorial, kd2, t02; + arb_t x1, x2, x3, x4, x5; + + arb_init(total); + arb_init(summand); + arb_init(pi); + arb_init(half); + arb_init(a); + arb_init(l_factorial); + arb_init(kd2); + arb_init(t02); + arb_init(x1); + arb_init(x2); + arb_init(x3); + arb_init(x4); + arb_init(x5); + + arb_one(half); + arb_mul_2exp_si(half, half, -1); + arb_const_pi(pi, prec); + arb_one(l_factorial); + arb_set_si(kd2, k); + arb_mul_2exp_si(kd2, kd2, -1); + arb_sqr(t02, t0, prec); + + for (l=0; l<=(sigma-1)/2; l++) + { + if (l > 1) + { + arb_mul_si(l_factorial, l_factorial, l, prec); + } + + arb_mul_si(a, pi, 4*l+1, prec); + arb_mul_si(a, a, A, prec); + + arb_inv(x1, a, prec); + arb_add_ui(x1, x1, 1, prec); + + arb_add_si(x2, half, 2*l, prec); + arb_sqr(x2, x2, prec); + arb_add(x2, x2, t02, prec); + arb_pow(x2, x2, kd2, prec); + arb_div(x2, x2, l_factorial, prec); + + arb_set_si(x3, 4*l + 1); + arb_div(x3, x3, h, prec); + arb_sqr(x3, x3, prec); + arb_mul_2exp_si(x3, x3, -3); + + arb_mul_2exp_si(x4, a, -1); + + arb_sub(x5, x3, x4, prec); + arb_exp(x5, x5, prec); + + arb_mul(summand, x1, x2, prec); + arb_mul(summand, summand, x5, prec); + + arb_add(total, total, summand, prec); + } + + arb_set(out, total); + + arb_clear(total); + arb_clear(summand); + arb_clear(pi); + arb_clear(half); + arb_clear(a); + arb_clear(l_factorial); + arb_clear(kd2); + arb_clear(t02); + arb_clear(x1); + arb_clear(x2); + arb_clear(x3); + arb_clear(x4); + arb_clear(x5); +} + +void +acb_dirichlet_platt_lemma_A7(arb_t out, slong sigma, + const arb_t t0, const arb_t h, slong k, slong A, slong prec) +{ + arb_t S, C; + arb_t pi, a; + arb_t x1, x2; + arb_t y1, y2, y3, y4; + arb_t z1, z2; + + if (sigma % 2 == 0 || sigma < 3) + { + arb_zero_pm_inf(out); + return; + } + + arb_init(S); + arb_init(C); + arb_init(pi); + arb_init(a); + arb_init(x1); + arb_init(x2); + arb_init(y1); + arb_init(y2); + arb_init(y3); + arb_init(y4); + arb_init(z1); + arb_init(z2); + + arb_const_pi(pi, prec); + + arb_pow_ui(x1, pi, (ulong) (k+1), prec); + arb_mul_2exp_si(x1, x1, k+3); + + arb_div(x2, t0, h, prec); + arb_sqr(x2, x2, prec); + arb_mul_2exp_si(x2, x2, -1); + arb_neg(x2, x2); + arb_exp(x2, x2, prec); + + _platt_lemma_A7_S(S, sigma, t0, h, k, A, prec); + + arb_mul(z1, x1, x2, prec); + arb_mul(z1, z1, S, prec); + + arb_mul_si(a, pi, 2*sigma-1, prec); + arb_mul_si(a, a, A, prec); + + arb_inv(y1, a, prec); + arb_add_ui(y1, y1, 1, prec); + + arb_set_si(y2, 2*sigma + 1); + arb_div(y2, y2, h, prec); + arb_sqr(y2, y2, prec); + arb_mul_2exp_si(y2, y2, -3); + + arb_mul_2exp_si(y3, a, -1); + + arb_sub(y4, y2, y3, prec); + arb_exp(y4, y4, prec); + + acb_dirichlet_platt_c_bound(C, sigma, t0, h, k, prec); + + arb_mul(z2, y1, y4, prec); + arb_mul(z2, z2, C, prec); + arb_mul_2exp_si(z2, z2, 1); + + arb_add(out, z1, z2, prec); + + arb_clear(S); + arb_clear(C); + arb_clear(pi); + arb_clear(a); + arb_clear(x1); + arb_clear(x2); + arb_clear(y1); + arb_clear(y2); + arb_clear(y3); + arb_clear(y4); + arb_clear(z1); + arb_clear(z2); +} diff --git a/acb_dirichlet/platt_lemma_A9.c b/acb_dirichlet/platt_lemma_A9.c new file mode 100644 index 00000000..6400b6f2 --- /dev/null +++ b/acb_dirichlet/platt_lemma_A9.c @@ -0,0 +1,162 @@ +/* + 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 . +*/ + +#include "acb_dirichlet.h" + +static void +_arb_pow_d(arb_t z, const arb_t x, double d, slong prec) +{ + arb_t u; + arb_init(u); + arb_set_d(u, d); + arb_pow(z, x, u, prec); + arb_clear(u); +} + +static void +_platt_lemma_A9_a(arb_t out, slong sigma, + const arb_t t0, const arb_t h, slong A, slong prec) +{ + arb_t a, pi, C; + arb_t y1, y2, y3, y4; + arb_t z1, z2, z3; + + arb_init(a); + arb_init(pi); + arb_init(C); + arb_init(y1); + arb_init(y2); + arb_init(y3); + arb_init(y4); + arb_init(z1); + arb_init(z2); + arb_init(z3); + + arb_const_pi(pi, prec); + + arb_mul_si(a, pi, 2*sigma-1, prec); + arb_mul_si(a, a, A, prec); + + arb_inv(y1, a, prec); + arb_add_ui(y1, y1, 1, prec); + + arb_set_si(y2, 2*sigma - 1); + arb_div(y2, y2, h, prec); + arb_sqr(y2, y2, prec); + arb_mul_2exp_si(y2, y2, -3); + + arb_mul_2exp_si(y3, a, -1); + + arb_sub(y4, y2, y3, prec); + arb_exp(y4, y4, prec); + + acb_dirichlet_platt_c_bound(C, sigma, t0, h, 0, prec); + + arb_zeta_ui(z1, (ulong) sigma, prec); + arb_mul_2exp_si(z1, z1, 1); + + arb_set_si(z2, 1-2*sigma); + arb_mul_2exp_si(z2, z2, -2); + arb_pow(z2, pi, z2, prec); + + arb_sub(z3, y2, y3, prec); + arb_exp(z3, z3, prec); + + arb_mul(out, z1, z2, prec); + arb_mul(out, out, z3, prec); + arb_mul(out, out, C, prec); + arb_mul(out, out, y1, prec); + + arb_clear(a); + arb_clear(pi); + arb_clear(C); + arb_clear(y1); + arb_clear(y2); + arb_clear(y3); + arb_clear(y4); + arb_clear(z1); + arb_clear(z2); + arb_clear(z3); +} + + +static void +_platt_lemma_A9_b(arb_t out, + const arb_t t0, const arb_t h, slong A, slong prec) +{ + arb_t pi; + arb_t x1, x2, x3, x4, x5; + + arb_init(pi); + arb_init(x1); + arb_init(x2); + arb_init(x3); + arb_init(x4); + arb_init(x5); + + arb_const_pi(pi, prec); + + _arb_pow_d(x1, pi, 1.25, prec); + arb_mul_2exp_si(x1, x1, 2); + + arb_sqr(x2, t0, prec); + arb_mul_2exp_si(x2, x2, 2); + arb_sub_ui(x2, x2, 1, prec); + arb_neg(x2, x2); + arb_div(x2, x2, h, prec); + arb_div(x2, x2, h, prec); + arb_mul_2exp_si(x2, x2, -3); + + arb_mul_si(x3, pi, A, prec); + arb_mul_2exp_si(x3, x3, -1); + + arb_mul_si(x4, pi, A, prec); + arb_inv(x4, x4, prec); + arb_add_ui(x4, x4, 1, prec); + + arb_sub(x5, x2, x3, prec); + arb_exp(x5, x5, prec); + + arb_mul(out, x1, x4, prec); + arb_mul(out, out, x5, prec); + + arb_clear(pi); + arb_clear(x1); + arb_clear(x2); + arb_clear(x3); + arb_clear(x4); + arb_clear(x5); +} + + +void +acb_dirichlet_platt_lemma_A9(arb_t out, slong sigma, + const arb_t t0, const arb_t h, slong A, slong prec) +{ + arb_t a, b; + + if (sigma % 2 == 0 || sigma < 3) + { + arb_zero_pm_inf(out); + return; + } + + arb_init(a); + arb_init(b); + + _platt_lemma_A9_a(a, sigma, t0, h, A, prec); + _platt_lemma_A9_b(b, t0, h, A, prec); + + arb_add(out, a, b, prec); + + arb_clear(a); + arb_clear(b); +} diff --git a/acb_dirichlet/platt_lemma_B1.c b/acb_dirichlet/platt_lemma_B1.c new file mode 100644 index 00000000..61f4b7d5 --- /dev/null +++ b/acb_dirichlet/platt_lemma_B1.c @@ -0,0 +1,69 @@ +/* + 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 . +*/ + +#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) +{ + arb_t pi, C; + arb_t x1, x2, x3; + + if (sigma % 2 == 0 || sigma < 3) + { + arb_zero_pm_inf(out); + return; + } + + arb_init(pi); + arb_init(C); + arb_init(x1); + arb_init(x2); + arb_init(x3); + + arb_const_pi(pi, prec); + acb_dirichlet_platt_c_bound(C, sigma, t0, h, 0, prec); + + arb_set_si(x1, 2*sigma - 1); + arb_div(x1, x1, h, prec); + arb_sqr(x1, x1, prec); + arb_mul_2exp_si(x1, x1, -3); + arb_exp(x1, x1, prec); + + arb_set_si(x2, 1 - 2*sigma); + arb_mul_2exp_si(x2, x2, -2); + arb_pow(x2, pi, x2, prec); + + arb_set_si(x3, 1 - sigma); + _arb_ui_pow_arb(x3, (ulong) J, x3, prec); + arb_div_si(x3, x3, sigma - 1, prec); + + arb_mul(out, x1, x2, prec); + arb_mul(out, out, x3, prec); + arb_mul(out, out, C, prec); + + arb_clear(pi); + arb_clear(C); + arb_clear(x1); + arb_clear(x2); + arb_clear(x3); +} diff --git a/acb_dirichlet/platt_lemma_B2.c b/acb_dirichlet/platt_lemma_B2.c new file mode 100644 index 00000000..672bb5e2 --- /dev/null +++ b/acb_dirichlet/platt_lemma_B2.c @@ -0,0 +1,62 @@ +/* + 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 . +*/ + +#include "acb_dirichlet.h" + +void +acb_dirichlet_platt_lemma_B2(arb_t out, slong K, const arb_t h, + const arb_t xi, slong prec) +{ + arb_t two, half, pi; + arb_t x1, x2, x3, x4, x5; + + arb_init(two); + arb_init(half); + arb_init(pi); + arb_init(x1); + arb_init(x2); + arb_init(x3); + arb_init(x4); + arb_init(x5); + + arb_set_ui(two, 2); + arb_mul_2exp_si(half, two, -2); + arb_const_pi(pi, prec); + + arb_set_si(x1, K + 5); + arb_mul_2exp_si(x1, x1, -1); + arb_pow(x1, two, x1, prec); + + arb_add_si(x2, half, K, prec); + arb_pow(x2, pi, x2, prec); + + arb_pow_ui(x3, h, (ulong) (K + 1), prec); + + arb_pow_ui(x4, xi, (ulong) K, prec); + + arb_set_si(x5, K + 2); + arb_mul_2exp_si(x5, x5, -1); + arb_rgamma(x5, x5, prec); + + arb_mul(out, x1, x2, prec); + arb_mul(out, out, x3, prec); + arb_mul(out, out, x4, prec); + arb_mul(out, out, x5, prec); + + arb_clear(two); + arb_clear(half); + arb_clear(pi); + arb_clear(x1); + arb_clear(x2); + arb_clear(x3); + arb_clear(x4); + arb_clear(x5); +} diff --git a/acb_dirichlet/platt_multieval.c b/acb_dirichlet/platt_multieval.c new file mode 100644 index 00000000..6d5e9c65 --- /dev/null +++ b/acb_dirichlet/platt_multieval.c @@ -0,0 +1,472 @@ +/* + 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 . +*/ + +#include "acb_dirichlet.h" +#include "arb_hypgeom.h" +#include "acb_calc.h" +#include "acb_dft.h" + + +static void +_arb_add_d(arb_t z, const arb_t x, double d, slong prec) +{ + arb_t u; + arb_init(u); + arb_set_d(u, d); + arb_add(z, x, u, prec); + arb_clear(u); +} + +static void +_arb_inv_si(arb_t z, slong n, slong prec) +{ + arb_set_si(z, n); + arb_inv(z, z, prec); +} + +static void +platt_g_gamma_term(acb_t out, const arb_t t0, const acb_t t, slong prec) +{ + acb_t z; + acb_init(z); + acb_add_arb(z, t, t0, prec); + acb_mul_onei(z, z); + acb_mul_2exp_si(z, z, 1); + acb_add_ui(z, z, 1, prec); + acb_mul_2exp_si(z, z, -2); + acb_gamma(out, z, prec); + acb_clear(z); +} + +static void +platt_g_exp_term(acb_t out, + const arb_t t0, const arb_t h, const acb_t t, slong prec) +{ + arb_t pi; + acb_t z1, z2; + arb_init(pi); + acb_init(z1); + acb_init(z2); + arb_const_pi(pi, prec); + acb_add_arb(z1, t, t0, prec); + acb_mul_arb(z1, z1, pi, prec); + acb_mul_2exp_si(z1, z1, -2); + acb_div_arb(z2, t, h, prec); + acb_sqr(z2, z2, prec); + acb_mul_2exp_si(z2, z2, -1); + acb_sub(out, z1, z2, prec); + acb_exp(out, out, prec); + arb_clear(pi); + acb_clear(z1); + acb_clear(z2); +} + +static void +platt_g_base(acb_t out, const acb_t t, slong prec) +{ + arb_t pi; + arb_init(pi); + arb_const_pi(pi, prec); + acb_mul_arb(out, t, pi, prec); + acb_mul_onei(out, out); + acb_mul_2exp_si(out, out, 1); + acb_neg(out, out); + arb_clear(pi); +} + + +static void +platt_g_table(acb_ptr table, slong A, slong B, + const arb_t t0, const arb_t h, slong K, slong prec) +{ + slong N = A*B; + slong i, n, k; + acb_t t, base; + acb_t gamma_term, exp_term, coeff; + acb_ptr precomputed_powers; + + acb_init(t); + acb_init(base); + acb_init(gamma_term); + acb_init(exp_term); + acb_init(coeff); + + precomputed_powers = _acb_vec_init(K); + + for (i=0; i 1) + { + break; + } + } + arb_swap(res, x); + arb_clear(x); + } } static void @@ -145,22 +165,6 @@ acb_dirichlet_platt_scaled_lambda_vec(arb_ptr res, } } - -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) @@ -211,7 +215,7 @@ _platt_bound_C3_Y(arb_t res, const arb_t t0, slong A, const arb_t H, 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); + _gamma_upper_workaround(g, g1, g2, 0, prec); /* res = a*b*A*H*g */ arb_mul_si(res, H, A, prec); @@ -254,7 +258,7 @@ _platt_bound_C3_Z(arb_t res, const arb_t t0, slong A, const arb_t H, 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); + _gamma_upper_workaround(g, g1, g2, 0, prec); /* res = a*b*A*g */ arb_mul_si(res, g, A, prec); @@ -293,7 +297,7 @@ acb_dirichlet_platt_bound_C3(arb_t res, const arb_t t0, slong A, const arb_t H, arb_exp(ee, ee, prec); if (!arb_gt(t0, ee)) { - arb_zero_pm_one(res); + arb_zero_pm_inf(res); goto finish; } @@ -302,13 +306,13 @@ acb_dirichlet_platt_bound_C3(arb_t res, const arb_t t0, slong A, const arb_t H, arb_mul_si(rhs, t0, A, prec); if (!arb_is_positive(lhs) || !arb_le(lhs, rhs)) { - arb_zero_pm_one(res); + arb_zero_pm_inf(res); goto finish; } /* res = (X + Y + Z) * 6 / (pi * Ns) */ arb_const_pi(pi, prec); - _platt_beta(beta, t0, prec); + acb_dirichlet_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); diff --git a/acb_dirichlet/test/t-platt_beta.c b/acb_dirichlet/test/t-platt_beta.c new file mode 100644 index 00000000..8ca8e697 --- /dev/null +++ b/acb_dirichlet/test/t-platt_beta.c @@ -0,0 +1,126 @@ +/* + 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 . +*/ + +#include "acb_dirichlet.h" + +static void +_arb_div_ui_ui(arb_t res, ulong a, ulong b, slong prec) +{ + arb_set_ui(res, a); + arb_div_ui(res, res, b, prec); +} + +static int +_arb_lt_d(const arb_t a, double d) +{ + int result; + arb_t x; + arb_init(x); + arb_set_d(x, d); + result = arb_lt(a, x); + arb_clear(x); + return result; +} + +int main() +{ + slong iter; + flint_rand_t state; + arb_t x, t, t0, expe; + + flint_printf("platt_beta...."); + fflush(stdout); + flint_randinit(state); + + arb_init(x); + arb_init(t); + arb_init(t0); + arb_init(expe); + + for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++) + { + arb_t beta, a, b, c; + acb_t z; + slong prec, mbits; + + prec = 2 + n_randint(state, 300); + mbits = 2 + n_randint(state, 20); + arb_randtest(t, state, prec, mbits); + arb_randtest(t0, state, prec, mbits); + arb_abs(t, t); + arb_abs(t0, t0); + arb_add(x, t, t0, prec); + arb_const_e(expe, prec); + arb_exp(expe, expe, prec); + + if (!arb_is_nonnegative(t) || !arb_gt(t0, expe) || !_arb_lt_d(x, 1e8)) + { + continue; + } + + arb_init(beta); + arb_init(a); + arb_init(b); + arb_init(c); + acb_init(z); + + acb_dirichlet_platt_beta(beta, t0, prec); + + /* Lemma A.10 in "Isolating some non-trivial zeros of zeta" */ + arb_pow(a, x, beta, prec); + arb_mul_ui(a, a, 3, prec); + acb_dirichlet_platt_scaled_lambda(c, t, prec); + arb_abs(c, c); + if (arb_gt(c, a)) + { + flint_printf("FAIL: Lemma A.10 |f(t)|\n\n"); + flint_printf("iter = %wd\n\n", iter); + flint_printf("prec = %wd\n\n", prec); + flint_printf("t = "); arb_printd(t, 15); flint_printf("\n\n"); + flint_printf("t0 = "); arb_printd(t0, 15); flint_printf("\n\n"); + flint_abort(); + } + + /* In Lemma A.10 proof in "Isolating some non-trivial zeros of zeta" */ + arb_pow(a, x, beta, prec); + _arb_div_ui_ui(b, 732, 1000, prec); + arb_mul(a, a, b, prec); + acb_set_d(z, 0.5); + arb_set(acb_imagref(z), x); + acb_zeta(z, z, prec); + acb_abs(c, z, prec); + if (arb_gt(c, a)) + { + flint_printf("FAIL: Lemma A.10 |zeta(1/2 + i(t + t0))|\n\n"); + flint_printf("iter = %wd\n\n", iter); + flint_printf("prec = %wd\n\n", prec); + flint_printf("t = "); arb_printd(t, 15); flint_printf("\n\n"); + flint_printf("t0 = "); arb_printd(t0, 15); flint_printf("\n\n"); + flint_abort(); + } + + arb_clear(beta); + arb_clear(a); + arb_clear(b); + arb_clear(c); + acb_clear(z); + } + + arb_clear(x); + arb_clear(t); + arb_clear(t0); + arb_clear(expe); + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/acb_dirichlet/test/t-platt_multieval.c b/acb_dirichlet/test/t-platt_multieval.c new file mode 100644 index 00000000..9ed40ee8 --- /dev/null +++ b/acb_dirichlet/test/t-platt_multieval.c @@ -0,0 +1,151 @@ +/* + 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 . +*/ + +#include "acb_dirichlet.h" + +static void +_arb_div_si_si(arb_t res, slong a, slong b, slong prec) +{ + arb_set_si(res, a); + arb_div_si(res, res, b, prec); +} + +static int +_arb_vec_overlaps(arb_srcptr a, arb_srcptr b, slong len) +{ + slong i; + for (i = 0; i < len; i++) + { + if (!arb_overlaps(a + i, b + i)) + { + return 0; + } + } + return 1; +} + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("platt_multieval...."); + fflush(stdout); + flint_randinit(state); + + /* Check a specific combination of parameter values that is relatively fast + * to evaluate and that has relatively tight bounds. */ + { + 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; + arb_t h; + arb_ptr vec; + + arb_init(h); + fmpz_init(T); + + fmpz_set_si(T, 10000); + arb_set_d(h, 4.5); + 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; + slong i = n_randint(state, N); + slong n = i - N/2; + arb_init(t); + arb_init(r); + _arb_div_si_si(t, n, A, prec); + arb_add_fmpz(t, t, T, prec); + acb_dirichlet_platt_scaled_lambda(r, t, prec); + 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("r = "); arb_printn(r, 30, 0); flint_printf("\n\n"); + flint_abort(); + } + arb_clear(t); + arb_clear(r); + } + + fmpz_clear(T); + arb_clear(h); + _arb_vec_clear(vec, N); + } + + for (iter = 0; iter < 10 * arb_test_multiplier(); iter++) + { + slong prec; + ulong A, B, N, J, K; + slong sigma, Tbits; + fmpz_t 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) */ + if (n_randint(state, 2)) + A *= 2; + else + B *= 2; + N = A*B; + + fmpz_init(T); + Tbits = 5 + n_randint(state, 15); + fmpz_set_ui(T, n_randtest_bits(state, Tbits)); + + arb_init(h); + arb_set_si(h, 1 + n_randint(state, 20000)); + arb_div_si(h, h, 1000, prec); + + v1 = _arb_vec_init(N); + v2 = _arb_vec_init(N); + acb_dirichlet_platt_scaled_lambda_vec(v1, T, A, B, prec); + acb_dirichlet_platt_multieval(v2, T, A, B, h, J, K, sigma, prec); + + if (!_arb_vec_overlaps(v1, v2, N)) + { + 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("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(T); + _arb_vec_clear(v1, N); + _arb_vec_clear(v2, N); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/doc/source/acb_dirichlet.rst b/doc/source/acb_dirichlet.rst index a3259ee2..0f2f55c1 100644 --- a/doc/source/acb_dirichlet.rst +++ b/doc/source/acb_dirichlet.rst @@ -756,10 +756,15 @@ 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) + 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. + `N=AB` must be an even integer. The multieval variant evaluates the + function at all points on the grid simultaneously using forward and inverse + discrete Fourier transforms, and it requires the four additional tuning + parameters *h*, *J*, *K*, and *sigma*. .. function:: void acb_dirichlet_platt_ws_interpolation(arb_t res, const arb_t t0, arb_srcptr p, const fmpz_t T, slong A, slong B, slong Ns_max, const arb_t H, slong sigma, slong prec)