mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
Merge pull request #279 from p15-git-acc/platt_multieval
platt multieval
This commit is contained in:
Fredrik Johansson
GitHub
commit
f44f23effa
15 changed files with 1742 additions and 26 deletions
|
|
@ -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
34
acb_dirichlet/platt_beta.c
Normal 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);
|
||||
}
|
||||
|
|
@ -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
50
acb_dirichlet/platt_lemma_32.c
Normal 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
227
acb_dirichlet/platt_lemma_A11.c
Normal 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
124
acb_dirichlet/platt_lemma_A5.c
Normal 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
180
acb_dirichlet/platt_lemma_A7.c
Normal 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
162
acb_dirichlet/platt_lemma_A9.c
Normal 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
69
acb_dirichlet/platt_lemma_B1.c
Normal 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
62
acb_dirichlet/platt_lemma_B2.c
Normal 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
472
acb_dirichlet/platt_multieval.c
Normal 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);
|
||||
}
|
||||
|
|
@ -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
126
acb_dirichlet/test/t-platt_beta.c
Normal 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
151
acb_dirichlet/test/t-platt_multieval.c
Normal 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;
|
||||
}
|
||||
|
|
@ -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)
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue