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

some tweaks to zeta and dirichlet evaluation

Browse source
This commit is contained in:
Fredrik Johansson 2017-11-21 17:38:38 +01:00
parent 4d34deaec0
commit 75d8f23521
3 changed files with 37 additions and 8 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

35
acb_dirichlet/l_euler_product.c
View file

@ -18,11 +18,11 @@ acb_dirichlet_l_euler_product(acb_t res, const acb_t s,
const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
{
arf_t left;
slong wp, powprec, left_s;
slong wp, powprec, left_s, acc;
ulong val, p, p_limit;
double p_needed_approx, powmag, logp, errmag;
int is_real;
acb_t t, u, v, c;
int is_real, is_int;
acb_t t, u, v, c, negs;
acb_dirichlet_roots_t roots;
mag_t err;
@ -60,6 +60,13 @@ acb_dirichlet_l_euler_product(acb_t res, const acb_t s,
left_s = arf_get_si(left, ARF_RND_FLOOR);
/* Adjust precision based on possible accuracy. */
acc = acb_rel_accuracy_bits(s);
acc = FLINT_MAX(acc, 0);
acc = FLINT_MIN(acc, prec);
acc += left_s;
prec = FLINT_MIN(prec, acc + 10);
/* Heuristic. */
wp = prec + FLINT_BIT_COUNT(prec) + (prec / left_s) + 4;
@ -82,7 +89,11 @@ acb_dirichlet_l_euler_product(acb_t res, const acb_t s,
acb_init(u);
acb_init(v);
acb_init(c);
acb_init(negs);
is_int = acb_is_int(s);
acb_neg(negs, s);
acb_one(v);
for (p = 2; p < p_limit; p = n_nextprime(p, 1))
@ -106,9 +117,20 @@ acb_dirichlet_l_euler_product(acb_t res, const acb_t s,
{
acb_dirichlet_root(c, roots, val, powprec);
acb_set_ui(t, p);
acb_pow(t, t, s, powprec);
acb_set_round(u, v, powprec);
acb_div(t, u, t, powprec);
if (is_int)
{
acb_pow(t, t, s, powprec);
acb_set_round(u, v, powprec);
acb_div(t, u, t, powprec);
}
else
{
acb_pow(t, t, negs, powprec);
acb_set_round(u, v, powprec);
acb_mul(t, u, t, powprec);
}
acb_mul(t, t, c, powprec);
acb_sub(v, v, t, wp);
}
@ -129,6 +151,7 @@ acb_dirichlet_l_euler_product(acb_t res, const acb_t s,
acb_clear(u);
acb_clear(v);
acb_clear(c);
acb_clear(negs);
arf_clear(left);
}

3
acb_dirichlet/zeta_jet_rs.c
View file

@ -115,12 +115,15 @@ acb_dirichlet_zeta_jet_rs(acb_ptr res, const acb_t s, slong len, slong prec)
{
acb_t t, u;
mag_t r, R, err0, err1, der1, der2, M;
/*
slong acc;
acc = acb_rel_accuracy_bits(s);
acc = FLINT_MAX(acc, 0);
acc = FLINT_MIN(acc, prec);
prec = FLINT_MIN(prec, acc + 20);
*/
/*
assume s = m +/- r

7
acb_dirichlet/zeta_rs_r.c
View file

@ -20,7 +20,7 @@ acb_dirichlet_zeta_rs_r(acb_t res, const acb_t s, slong K, slong prec)
acb_t U, S, u, v;
fmpz_t N;
mag_t err;
slong j, k, wp;
slong j, k, wp, K_limit;
/* determinate K automatically */
if (K <= 0)
@ -41,7 +41,10 @@ acb_dirichlet_zeta_rs_r(acb_t res, const acb_t s, slong K, slong prec)
best_log2err = 1e300;
/* todo: also break if too slow rate of decay? */
for (K = 1; K < 10 + prec * 0.25; K++)
K_limit = 10 + prec * 0.25;
K_limit += pow(t, 0.2); /* possibly useful for off-strip evaluation */
for (K = 1; K < K_limit; K++)
{
if (sigma < 0 && K + sigma < 3)
continue;