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

Merge branch 'master' of github.com:fredrik-johansson/arb

Browse source
This commit is contained in:
fredrik 2019-06-20 08:11:00 +02:00
commit 1875b0294c
15 changed files with 1742 additions and 26 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

23
acb_dirichlet.h
View file

@ -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);

34
acb_dirichlet/platt_beta.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

29
acb_dirichlet/platt_c_bound.c
View file

@ -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);
}

50
acb_dirichlet/platt_lemma_32.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

227
acb_dirichlet/platt_lemma_A11.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

124
acb_dirichlet/platt_lemma_A5.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

180
acb_dirichlet/platt_lemma_A7.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

162
acb_dirichlet/platt_lemma_A9.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

69
acb_dirichlet/platt_lemma_B1.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

62
acb_dirichlet/platt_lemma_B2.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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);
}

472
acb_dirichlet/platt_multieval.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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<N; i++)
{
n = i - N/2;
acb_set_si(t, n);
acb_div_si(t, t, A, prec);
platt_g_base(base, t, prec);
_acb_vec_set_powers(precomputed_powers, base, K, prec);
platt_g_gamma_term(gamma_term, t0, t, prec);
platt_g_exp_term(exp_term, t0, h, t, prec);
acb_mul(coeff, gamma_term, exp_term, prec);
for (k=0; k<K; k++)
{
acb_mul(table + k*N + i, coeff, precomputed_powers + k, prec);
}
}
acb_clear(t);
acb_clear(base);
acb_clear(gamma_term);
acb_clear(exp_term);
acb_clear(coeff);
_acb_vec_clear(precomputed_powers, K);
}
static void
logjsqrtpi(arb_t out, slong j, slong prec)
{
arb_const_sqrt_pi(out, prec);
arb_mul_si(out, out, j, prec);
arb_log(out, out, prec);
}
static void
_acb_add_error_arb_mag(acb_t res, const arb_t x)
{
mag_t err;
mag_init(err);
arb_get_mag(err, x);
acb_add_error_mag(res, err);
mag_clear(err);
}
static void
_acb_vec_scalar_add_error_mag(acb_ptr res, slong len, const mag_t err)
{
slong i;
for (i = 0; i < len; i++)
{
acb_add_error_mag(res + i, err);
}
}
static void
_acb_vec_scalar_add_error_arb_mag(acb_ptr res, slong len, const arb_t x)
{
mag_t err;
mag_init(err);
arb_get_mag(err, x);
_acb_vec_scalar_add_error_mag(res, len, err);
mag_clear(err);
}
static void
get_smk_index(slong *out, slong B, slong j, slong prec)
{
arb_t pi, x;
fmpz_t z;
arb_init(pi);
arb_init(x);
fmpz_init(z);
while (1)
{
arb_const_pi(pi, prec);
logjsqrtpi(x, j, prec);
arb_div(x, x, pi, prec);
arb_mul_2exp_si(x, x, -1);
arb_mul_si(x, x, B, prec);
_arb_add_d(x, x, 0.5, prec);
arb_floor(x, x, prec);
if (arb_get_unique_fmpz(z, x))
{
*out = fmpz_get_si(z);
break;
}
else
{
prec *= 2;
}
}
arb_clear(pi);
arb_clear(x);
fmpz_clear(z);
}
static void
platt_smk(acb_ptr table,
const arb_t t0, slong A, slong B, slong J, slong K, slong prec)
{
slong j, k, m;
slong N = A * B;
acb_ptr row;
arb_ptr diff_powers;
arb_t rpi, rsqrtj, um, a, base;
acb_t z;
arb_init(rpi);
arb_init(rsqrtj);
arb_init(um);
arb_init(a);
arb_init(base);
acb_init(z);
diff_powers = _arb_vec_init(K);
arb_const_pi(rpi, prec);
arb_inv(rpi, rpi, prec);
for (j = 1; j <= J; j++)
{
logjsqrtpi(a, j, prec);
arb_mul(a, a, rpi, prec);
arb_rsqrt_ui(rsqrtj, (ulong) j, prec);
acb_set_arb(z, t0);
acb_mul_arb(z, z, a, prec);
acb_neg(z, z);
acb_exp_pi_i(z, z, prec);
acb_mul_arb(z, z, rsqrtj, prec);
get_smk_index(&m, B, j, prec);
arb_set_si(um, m);
arb_div_si(um, um, B, prec);
arb_mul_2exp_si(base, a, -1);
arb_sub(base, base, um, prec);
_arb_vec_set_powers(diff_powers, base, K, prec);
for (k = 0; k < K; k++)
{
row = table + N*k;
acb_addmul_arb(row + m, z, diff_powers + k, prec);
}
}
arb_clear(rpi);
arb_clear(rsqrtj);
arb_clear(um);
arb_clear(a);
arb_clear(base);
acb_clear(z);
_arb_vec_clear(diff_powers, K);
}
static void
do_convolutions(acb_ptr out_table,
const acb_ptr table, const acb_ptr S_table,
slong N, slong K, slong prec)
{
slong i, k;
acb_ptr padded_table_row, padded_S_table_row, padded_out_table;
acb_ptr fp, gp;
acb_dft_pre_t pre;
padded_table_row = _acb_vec_init(N*2);
padded_S_table_row = _acb_vec_init(N*2);
padded_out_table = _acb_vec_init(N*2);
fp = _acb_vec_init(N*2);
gp = _acb_vec_init(N*2);
acb_dft_precomp_init(pre, N*2, prec);
for (k = 0; k < K; k++)
{
_acb_vec_zero(padded_table_row, N*2);
_acb_vec_zero(padded_S_table_row, N*2);
_acb_vec_zero(padded_out_table, N*2);
_acb_vec_set(padded_table_row, table + k*N, N);
_acb_vec_set(padded_S_table_row, S_table + k*N, N);
for (i = 1; i < N; i++)
{
acb_swap(padded_S_table_row + i, padded_S_table_row + N*2 - i);
}
acb_dft_precomp(fp, padded_S_table_row, pre, prec);
acb_dft_precomp(gp, padded_table_row, pre, prec);
_acb_vec_kronecker_mul(gp, gp, fp, N*2, prec);
acb_dft_inverse_precomp(padded_out_table, gp, pre, prec);
for (i = 0; i <= N/2; i++)
{
acb_add(out_table + i,
out_table + i, padded_out_table + i, prec);
}
}
_acb_vec_clear(padded_table_row, N*2);
_acb_vec_clear(padded_S_table_row, N*2);
_acb_vec_clear(padded_out_table, N*2);
_acb_vec_clear(fp, N*2);
_acb_vec_clear(gp, N*2);
acb_dft_precomp_clear(pre);
}
static void
remove_gaussian_window(arb_ptr out, slong A, slong B, const arb_t h, slong prec)
{
slong i, n;
slong N = A*B;
arb_t t, x;
arb_init(t);
arb_init(x);
for (i = 0; i < N; i++)
{
n = i - N/2;
arb_set_si(t, n);
arb_div_si(t, t, A, prec);
arb_div(x, t, h, prec);
arb_sqr(x, x, prec);
arb_mul_2exp_si(x, x, -1);
arb_exp(x, x, prec);
arb_mul(out + i, out + i, x, prec);
}
arb_clear(t);
arb_clear(x);
}
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 N = A*B;
slong i, k;
acb_ptr table, S_table, out_a, out_b;
acb_ptr row;
arb_t t0, t, x, k_factorial, err, ratio, c, xi;
acb_t z;
acb_dft_pre_t pre_N;
arb_init(t0);
arb_init(t);
arb_init(x);
arb_init(k_factorial);
arb_init(err);
arb_init(ratio);
arb_init(c);
arb_init(xi);
acb_init(z);
table = _acb_vec_init(K*N);
S_table = _acb_vec_init(K*N);
out_a = _acb_vec_init(N);
out_b = _acb_vec_init(N);
acb_dft_precomp_init(pre_N, N, prec);
arb_set_fmpz(t0, T);
_arb_inv_si(xi, B, prec);
arb_mul_2exp_si(xi, xi, -1);
platt_smk(S_table, t0, A, B, J, K, prec);
platt_g_table(table, A, B, t0, h, K, prec);
for (k = 0; k < K; k++)
{
acb_dirichlet_platt_lemma_A5(err, B, h, k, prec);
_acb_vec_scalar_add_error_arb_mag(table + N*k, N, err);
}
for (k = 0; k < K; k++)
{
row = table + N*k;
for (i = 0; i < N/2; i++)
{
acb_swap(row + i, row + i + N/2);
}
acb_dft_precomp(row, row, pre_N, prec);
}
_acb_vec_scalar_div_ui(table, table, N*K, (ulong) A, prec);
for (k = 0; k < K; k++)
{
acb_dirichlet_platt_lemma_A7(err, sigma, t0, h, k, A, prec);
_acb_vec_scalar_add_error_arb_mag(table + N*k, N, err);
}
arb_one(k_factorial);
for (k = 2; k < K; k++)
{
row = table + N*k;
arb_mul_ui(k_factorial, k_factorial, (ulong) k, prec);
_acb_vec_scalar_div_arb(row, row, N, k_factorial, prec);
}
do_convolutions(out_a, table, S_table, N, K, prec);
for (i = 0; i < N/2 + 1; i++)
{
arb_set_si(x, i);
arb_div_si(x, x, B, prec);
acb_dirichlet_platt_lemma_32(err, h, t0, x, prec);
_acb_add_error_arb_mag(out_a + i, err);
}
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_mul_2exp_si(c, c, 1);
arb_sub_ui(c, c, 1, prec);
acb_dirichlet_platt_lemma_B2(err, K, h, xi, prec);
arb_mul(err, err, c, prec);
_acb_vec_scalar_add_error_arb_mag(out_a, N/2 + 1, err);
for (i = 1; i < N/2; i++)
{
acb_conj(out_a + N - i, out_a + i);
}
acb_dirichlet_platt_lemma_A9(err, sigma, t0, h, A, prec);
_acb_vec_scalar_add_error_arb_mag(out_a, N, err);
acb_dft_inverse_precomp(out_b, out_a, pre_N, prec);
_acb_vec_scalar_mul_ui(out_b, out_b, N, (ulong) A, prec);
for (i = 0; i < N/2; i++)
{
acb_swap(out_b + i, out_b + i + N/2);
}
acb_dirichlet_platt_lemma_A11(err, t0, h, B, prec);
_acb_vec_scalar_add_error_arb_mag(out_b, N, err);
for (i = 0; i < N; i++)
{
arb_swap(out + i, acb_realref(out_b + i));
}
remove_gaussian_window(out, A, B, h, prec);
arb_clear(t0);
arb_clear(t);
arb_clear(x);
arb_clear(k_factorial);
arb_clear(err);
arb_clear(ratio);
arb_clear(c);
arb_clear(xi);
acb_clear(z);
_acb_vec_clear(table, K*N);
_acb_vec_clear(S_table, K*N);
_acb_vec_clear(out_a, N);
_acb_vec_clear(out_b, N);
acb_dft_precomp_clear(pre_N);
}

52
acb_dirichlet/platt_ws_interpolation.c
View file

@ -12,11 +12,31 @@
#include "acb_dirichlet.h"
#include "arb_hypgeom.h"
/* Increase precision adaptively. */
static void
_arb_inv_ui(arb_t res, ulong n, slong prec)
_gamma_upper_workaround(arb_t res, const arb_t s, const arb_t z,
int regularized, slong prec)
{
arb_set_ui(res, n);
arb_inv(res, res, 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
@ -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);

126
acb_dirichlet/test/t-platt_beta.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

151
acb_dirichlet/test/t-platt_multieval.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

7
doc/source/acb_dirichlet.rst
View file

@ -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)