mirror of
https://github.com/vale981/arb
synced 2025-03-04 17:01:40 -05:00
add acb_dirichlet_zeta_jet + first derivative using Riemann-Siegel; also fix a bug in Riemann-Siegel code
This commit is contained in:
Fredrik Johansson
parent
0fbf524a22
commit
537fae5110
7 changed files with 476 additions and 96 deletions
|
|
@ -42,6 +42,9 @@ void acb_dirichlet_zeta_rs_r(acb_t res, const acb_t s, slong K, slong prec);
|
|||
void acb_dirichlet_zeta_rs(acb_t res, const acb_t s, slong K, slong prec);
|
||||
void acb_dirichlet_zeta(acb_t res, const acb_t s, slong prec);
|
||||
|
||||
void acb_dirichlet_zeta_jet_rs(acb_ptr res, const acb_t s, slong len, slong prec);
|
||||
void acb_dirichlet_zeta_jet(acb_t res, const acb_t s, int deflate, slong len, slong prec);
|
||||
|
||||
void acb_dirichlet_hurwitz(acb_t res, const acb_t s, const acb_t a, slong prec);
|
||||
|
||||
typedef struct
|
||||
|
|
|
|||
|
|
@ -36,16 +36,9 @@ acb_dirichlet_l_jet(acb_ptr res, const acb_t s,
|
|||
if (G == NULL || G->q == 1)
|
||||
{
|
||||
if (len == 1 && !deflate)
|
||||
{
|
||||
acb_dirichlet_zeta(res, s, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_init(a);
|
||||
acb_one(a);
|
||||
_acb_poly_zeta_cpx_reflect(res, s, a, deflate, len, prec);
|
||||
acb_clear(a);
|
||||
}
|
||||
acb_dirichlet_zeta_jet(res, s, deflate, len, prec);
|
||||
return;
|
||||
}
|
||||
|
||||
|
|
|
|||
113
acb_dirichlet/test/t-zeta_jet_rs.c
Normal file
113
acb_dirichlet/test/t-zeta_jet_rs.c
Normal file
|
|
@ -0,0 +1,113 @@
|
|||
/*
|
||||
Copyright (C) 2017 Fredrik Johansson
|
||||
|
||||
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"
|
||||
|
||||
int main()
|
||||
{
|
||||
slong iter;
|
||||
flint_rand_t state;
|
||||
|
||||
flint_printf("zeta_jet_rs....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
for (iter = 0; iter < 200 * arb_test_multiplier(); iter++)
|
||||
{
|
||||
acb_t s1, s2, a;
|
||||
acb_ptr r1, r2;
|
||||
slong prec1, prec2, len;
|
||||
|
||||
acb_init(a);
|
||||
acb_init(s1);
|
||||
acb_init(s2);
|
||||
r1 = _acb_vec_init(2);
|
||||
r2 = _acb_vec_init(2);
|
||||
|
||||
len = 1 + n_randint(state, 2);
|
||||
|
||||
if (n_randint(state, 2))
|
||||
arb_set_d(acb_realref(s1), 0.5);
|
||||
else
|
||||
arb_randtest(acb_realref(s1), state, 2 + n_randint(state, 100), 2);
|
||||
|
||||
arb_randtest(acb_imagref(s1), state, 2 + n_randint(state, 100), 2);
|
||||
arb_add_ui(acb_imagref(s1), acb_imagref(s1), n_randtest(state) % 1000, 100);
|
||||
acb_set(s2, s1);
|
||||
|
||||
acb_randtest(r1, state, 2 + n_randint(state, 100), 3);
|
||||
acb_randtest(r2, state, 2 + n_randint(state, 100), 3);
|
||||
|
||||
if (len > 1)
|
||||
{
|
||||
acb_randtest(r1 + 1, state, 2 + n_randint(state, 100), 3);
|
||||
acb_randtest(r2 + 1, state, 2 + n_randint(state, 100), 3);
|
||||
}
|
||||
|
||||
prec1 = 2 + n_randint(state, 150);
|
||||
prec2 = 2 + n_randint(state, 150);
|
||||
|
||||
if (n_randint(state, 4) == 0)
|
||||
{
|
||||
mag_zero(arb_radref(acb_realref(s2)));
|
||||
mag_zero(arb_radref(acb_imagref(s2)));
|
||||
}
|
||||
|
||||
if (n_randint(state, 4) == 0)
|
||||
{
|
||||
if (n_randint(state, 2))
|
||||
arb_get_ubound_arf(arb_midref(acb_realref(s2)), acb_realref(s2), ARF_PREC_EXACT);
|
||||
else
|
||||
arb_get_lbound_arf(arb_midref(acb_realref(s2)), acb_realref(s2), ARF_PREC_EXACT);
|
||||
mag_zero(arb_radref(acb_realref(s2)));
|
||||
}
|
||||
|
||||
if (n_randint(state, 4) == 0)
|
||||
{
|
||||
if (n_randint(state, 2))
|
||||
arb_get_ubound_arf(arb_midref(acb_imagref(s2)), acb_imagref(s2), ARF_PREC_EXACT);
|
||||
else
|
||||
arb_get_lbound_arf(arb_midref(acb_imagref(s2)), acb_imagref(s2), ARF_PREC_EXACT);
|
||||
mag_zero(arb_radref(acb_imagref(s2)));
|
||||
}
|
||||
|
||||
acb_dirichlet_zeta_jet_rs(r1, s1, len, prec1);
|
||||
acb_one(a);
|
||||
|
||||
_acb_poly_zeta_cpx_series(r2, s2, a, 0, len, prec2);
|
||||
|
||||
if (!acb_overlaps(r1, r2) || !acb_overlaps(r1 + 1, r2 + 1))
|
||||
{
|
||||
flint_printf("FAIL: overlap\n\n");
|
||||
flint_printf("iter = %wd, len = %wd\n", iter, len);
|
||||
flint_printf("s1 = "); acb_printn(s1, 50, 0); flint_printf("\n\n");
|
||||
flint_printf("s2 = "); acb_printn(s2, 50, 0); flint_printf("\n\n");
|
||||
flint_printf("r1[0] = "); acb_printn(r1, 50, 0); flint_printf("\n\n");
|
||||
flint_printf("r2[0] = "); acb_printn(r2, 50, 0); flint_printf("\n\n");
|
||||
flint_printf("r1[1] = "); acb_printn(r1 + 1, 50, 0); flint_printf("\n\n");
|
||||
flint_printf("r2[1] = "); acb_printn(r2 + 1, 50, 0); flint_printf("\n\n");
|
||||
flint_abort();
|
||||
}
|
||||
|
||||
acb_clear(a);
|
||||
acb_clear(s1);
|
||||
acb_clear(s2);
|
||||
_acb_vec_clear(r1, 2);
|
||||
_acb_vec_clear(r2, 2);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
||||
143
acb_dirichlet/zeta_jet.c
Normal file
143
acb_dirichlet/zeta_jet.c
Normal file
|
|
@ -0,0 +1,143 @@
|
|||
/*
|
||||
Copyright (C) 2017 Fredrik Johansson
|
||||
|
||||
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_zeta_jet(acb_t t, const acb_t h, int deflate, slong len, slong prec)
|
||||
{
|
||||
acb_t a;
|
||||
acb_init(a);
|
||||
acb_one(a);
|
||||
|
||||
/* use reflection formula */
|
||||
if (arf_sgn(arb_midref(acb_realref(h))) < 0)
|
||||
{
|
||||
/* zeta(s) = (2*pi)**s * sin(pi*s/2) / pi * gamma(1-s) * zeta(1-s) */
|
||||
acb_t pi, hcopy;
|
||||
acb_ptr f, s1, s2, s3, s4, u;
|
||||
slong i;
|
||||
|
||||
acb_init(pi);
|
||||
acb_init(hcopy);
|
||||
f = _acb_vec_init(2);
|
||||
s1 = _acb_vec_init(len);
|
||||
s2 = _acb_vec_init(len);
|
||||
s3 = _acb_vec_init(len);
|
||||
s4 = _acb_vec_init(len);
|
||||
u = _acb_vec_init(len);
|
||||
acb_set(hcopy, h);
|
||||
|
||||
acb_const_pi(pi, prec);
|
||||
|
||||
/* s1 = (2*pi)**s */
|
||||
acb_mul_2exp_si(pi, pi, 1);
|
||||
_acb_poly_pow_cpx(s1, pi, h, len, prec);
|
||||
acb_mul_2exp_si(pi, pi, -1);
|
||||
|
||||
/* s2 = sin(pi*s/2) / pi */
|
||||
acb_set(f, h);
|
||||
acb_one(f + 1);
|
||||
acb_mul_2exp_si(f, f, -1);
|
||||
acb_mul_2exp_si(f + 1, f + 1, -1);
|
||||
_acb_poly_sin_pi_series(s2, f, 2, len, prec);
|
||||
_acb_vec_scalar_div(s2, s2, len, pi, prec);
|
||||
|
||||
/* s3 = gamma(1-s) */
|
||||
acb_sub_ui(f, hcopy, 1, prec);
|
||||
acb_neg(f, f);
|
||||
acb_set_si(f + 1, -1);
|
||||
_acb_poly_gamma_series(s3, f, 2, len, prec);
|
||||
|
||||
/* s4 = zeta(1-s) */
|
||||
acb_sub_ui(f, hcopy, 1, prec);
|
||||
acb_neg(f, f);
|
||||
_acb_poly_zeta_cpx_series(s4, f, a, 0, len, prec);
|
||||
for (i = 1; i < len; i += 2)
|
||||
acb_neg(s4 + i, s4 + i);
|
||||
|
||||
_acb_poly_mullow(u, s1, len, s2, len, len, prec);
|
||||
_acb_poly_mullow(s1, s3, len, s4, len, len, prec);
|
||||
_acb_poly_mullow(t, u, len, s1, len, len, prec);
|
||||
|
||||
/* add 1/(1-(s+t)) = 1/(1-s) + t/(1-s)^2 + ... */
|
||||
if (deflate)
|
||||
{
|
||||
acb_sub_ui(u, hcopy, 1, prec);
|
||||
acb_neg(u, u);
|
||||
acb_inv(u, u, prec);
|
||||
for (i = 1; i < len; i++)
|
||||
acb_mul(u + i, u + i - 1, u, prec);
|
||||
_acb_vec_add(t, t, u, len, prec);
|
||||
}
|
||||
|
||||
acb_clear(pi);
|
||||
acb_clear(hcopy);
|
||||
_acb_vec_clear(f, 2);
|
||||
_acb_vec_clear(s1, len);
|
||||
_acb_vec_clear(s2, len);
|
||||
_acb_vec_clear(s3, len);
|
||||
_acb_vec_clear(s4, len);
|
||||
_acb_vec_clear(u, len);
|
||||
}
|
||||
else
|
||||
{
|
||||
_acb_poly_zeta_cpx_series(t, h, a, deflate, len, prec);
|
||||
}
|
||||
|
||||
acb_clear(a);
|
||||
}
|
||||
|
||||
/* todo: should adjust precision to input accuracy */
|
||||
void
|
||||
acb_dirichlet_zeta_jet(acb_t res, const acb_t s, int deflate, slong len, slong prec)
|
||||
{
|
||||
double cutoff;
|
||||
|
||||
if (len == 1 && deflate == 0)
|
||||
{
|
||||
acb_zeta(res, s, prec);
|
||||
return;
|
||||
}
|
||||
|
||||
if (deflate == 0 && (arb_contains_zero(acb_imagref(s))
|
||||
&& arb_contains_si(acb_realref(s), 1)))
|
||||
{
|
||||
_acb_vec_indeterminate(res, len);
|
||||
return;
|
||||
}
|
||||
|
||||
if (len > 2 || deflate != 0)
|
||||
{
|
||||
_acb_dirichlet_zeta_jet(res, s, deflate, len, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
cutoff = 24.0 * prec * sqrt(prec);
|
||||
|
||||
if (arb_is_exact(acb_realref(s)) &&
|
||||
arf_cmp_2exp_si(arb_midref(acb_realref(s)), -1) == 0)
|
||||
cutoff *= 2.5;
|
||||
else
|
||||
cutoff *= 4.0;
|
||||
|
||||
if (arf_cmpabs_d(arb_midref(acb_imagref(s)), cutoff) >= 0 &&
|
||||
arf_cmpabs_d(arb_midref(acb_realref(s)), 10 + prec * 0.1) <= 0)
|
||||
{
|
||||
acb_dirichlet_zeta_jet_rs(res, s, len, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
_acb_dirichlet_zeta_jet(res, s, deflate, len, prec);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
199
acb_dirichlet/zeta_jet_rs.c
Normal file
199
acb_dirichlet/zeta_jet_rs.c
Normal file
|
|
@ -0,0 +1,199 @@
|
|||
/*
|
||||
Copyright (C) 2017 Fredrik Johansson
|
||||
|
||||
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"
|
||||
|
||||
/*
|
||||
f(s) = (f(s+hi) + f(s-hi))/2 + err0
|
||||
f'(s) = (f(s+hi) - f(s-hi))/(2hi) + err1
|
||||
|
||||
|err0| <= |f(s +/- R)| (1/R)^2 h^2 / (1 - h/R)
|
||||
|err1| <= |f(s +/- R)| (1/R)^3 h^2 / (1 - h/R)
|
||||
|
||||
assuming h/R < 1
|
||||
*/
|
||||
static void
|
||||
acb_dirichlet_zeta_jet_rs_mid(acb_ptr res, const acb_t s, slong prec)
|
||||
{
|
||||
acb_t t, u;
|
||||
arb_t hh;
|
||||
mag_t h, R, err0, err1, tmp;
|
||||
slong hexp;
|
||||
|
||||
acb_init(t);
|
||||
acb_init(u);
|
||||
arb_init(hh);
|
||||
mag_init(h);
|
||||
mag_init(R);
|
||||
mag_init(err0);
|
||||
mag_init(err1);
|
||||
mag_init(tmp);
|
||||
|
||||
hexp = -(prec / 2) - 10;
|
||||
mag_set_ui_2exp_si(h, 1, hexp);
|
||||
mag_set_ui_2exp_si(R, 1, -5);
|
||||
|
||||
acb_set(t, s);
|
||||
mag_add(arb_radref(acb_realref(t)), arb_radref(acb_realref(t)), R);
|
||||
mag_add(arb_radref(acb_imagref(t)), arb_radref(acb_imagref(t)), R);
|
||||
/* tmp = |f(s +/- R)| */
|
||||
acb_dirichlet_zeta_bound(tmp, t);
|
||||
|
||||
/* err0 = (1/R)^2 h^2 / (1 - h/R) */
|
||||
mag_div(err0, h, R);
|
||||
mag_geom_series(err0, err0, 2);
|
||||
/* err0 *= tmp */
|
||||
mag_mul(err0, err0, tmp);
|
||||
/* err1 = err0 / R */
|
||||
mag_div(err1, err0, R);
|
||||
|
||||
/* hh = h as an arb_t */
|
||||
arb_one(hh);
|
||||
arb_mul_2exp_si(hh, hh, hexp);
|
||||
|
||||
acb_set(t, s);
|
||||
acb_set(u, s);
|
||||
arb_add(acb_imagref(t), acb_imagref(t), hh, 10 * prec);
|
||||
arb_sub(acb_imagref(u), acb_imagref(u), hh, 10 * prec);
|
||||
|
||||
/* zeta(mid(s)+ih) */
|
||||
acb_dirichlet_zeta_rs(res, t, 0, 1.5 * prec + 10);
|
||||
/* zeta(mid(s)-ih) */
|
||||
acb_dirichlet_zeta_rs(res + 1, u, 0, 1.5 * prec + 10);
|
||||
|
||||
acb_sub(t, res, res + 1, prec);
|
||||
acb_add(res, res, res + 1, prec);
|
||||
acb_swap(res + 1, t);
|
||||
|
||||
acb_mul_2exp_si(res, res, -1);
|
||||
acb_mul_2exp_si(res + 1, res + 1, -1);
|
||||
acb_mul_2exp_si(res + 1, res + 1, -hexp);
|
||||
acb_div_onei(res + 1, res + 1);
|
||||
|
||||
acb_add_error_mag(res, err0);
|
||||
acb_add_error_mag(res + 1, err1);
|
||||
|
||||
acb_clear(t);
|
||||
acb_clear(u);
|
||||
arb_clear(hh);
|
||||
mag_clear(h);
|
||||
mag_clear(R);
|
||||
mag_clear(err0);
|
||||
mag_clear(err1);
|
||||
mag_clear(tmp);
|
||||
}
|
||||
|
||||
void
|
||||
acb_dirichlet_zeta_jet_rs(acb_ptr res, const acb_t s, slong len, slong prec)
|
||||
{
|
||||
if (len > 2)
|
||||
{
|
||||
flint_printf("acb_dirichlet_zeta_jet_rs: len > 2 not implemented\n");
|
||||
flint_abort();
|
||||
}
|
||||
|
||||
if (len <= 0)
|
||||
return;
|
||||
|
||||
if (len == 1)
|
||||
{
|
||||
acb_dirichlet_zeta_rs(res, s, 0, prec);
|
||||
}
|
||||
else if (acb_is_exact(s))
|
||||
{
|
||||
acb_dirichlet_zeta_jet_rs_mid(res, s, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
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
|
||||
|
||||
f(s) = f(m) + err0
|
||||
f'(s) = f'(m) + err1
|
||||
|
||||
|err1| <= |f''(s)| r
|
||||
|err0| <= min(|f'(s)| r, |f'(m)| r + 0.5 |f''(s)| r^2)
|
||||
= r min(|f'(s)|, |f'(m)| + 0.5 |f''(s)| r)
|
||||
|
||||
|f'(s)| <= |f(s +/- R)| / R
|
||||
|f''(s)| <= 2 |f(s +/- R)| / R^2
|
||||
*/
|
||||
|
||||
acb_init(t);
|
||||
acb_init(u);
|
||||
mag_init(r);
|
||||
mag_init(R);
|
||||
mag_init(err0);
|
||||
mag_init(err1);
|
||||
mag_init(der1);
|
||||
mag_init(der2);
|
||||
mag_init(M);
|
||||
|
||||
/* r = rad(s) */
|
||||
mag_hypot(r, arb_radref(acb_realref(s)), arb_radref(acb_imagref(s)));
|
||||
|
||||
/* R = 1/8 */
|
||||
mag_set_ui_2exp_si(R, 1, -3);
|
||||
|
||||
/* t = s +/- R */
|
||||
acb_set(t, s);
|
||||
mag_add(arb_radref(acb_realref(t)), arb_radref(acb_realref(t)), R);
|
||||
mag_add(arb_radref(acb_imagref(t)), arb_radref(acb_imagref(t)), R);
|
||||
/* M = |f(s +/- R)| */
|
||||
acb_dirichlet_zeta_bound(M, t);
|
||||
/* der1 = |f'(s)| */
|
||||
mag_div(der1, M, R);
|
||||
/* der2 = |f''(s)| */
|
||||
mag_div(der2, der1, R);
|
||||
mag_mul_2exp_si(der2, der2, 1);
|
||||
|
||||
/* f(m), f'(m) */
|
||||
acb_set(t, s);
|
||||
mag_zero(arb_radref(acb_realref(t)));
|
||||
mag_zero(arb_radref(acb_imagref(t)));
|
||||
acb_dirichlet_zeta_jet_rs_mid(res, t, prec);
|
||||
|
||||
/* err1 = |f''(s)| r */
|
||||
mag_mul(err1, der2, r);
|
||||
|
||||
/* err0 = |f'(m)| + 0.5 |f''(s)| r */
|
||||
acb_get_mag(M, res + 1);
|
||||
mag_mul_2exp_si(err0, err1, -1);
|
||||
mag_add(err0, err0, M);
|
||||
/* err0 = min(err0, |f'(s)| */
|
||||
mag_min(err0, err0, der1);
|
||||
/* err0 = err0 * r */
|
||||
mag_mul(err0, err0, r);
|
||||
|
||||
acb_add_error_mag(res, err0);
|
||||
acb_add_error_mag(res + 1, err1);
|
||||
|
||||
acb_clear(t);
|
||||
acb_clear(u);
|
||||
mag_clear(r);
|
||||
mag_clear(R);
|
||||
mag_clear(err0);
|
||||
mag_clear(err1);
|
||||
mag_clear(der1);
|
||||
mag_clear(der2);
|
||||
mag_clear(M);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -92,6 +92,12 @@ acb_dirichlet_zeta_rs(acb_t res, const acb_t s, slong K, slong prec)
|
|||
{
|
||||
acb_t t;
|
||||
mag_t rad, R, err;
|
||||
slong acc;
|
||||
|
||||
acc = acb_rel_accuracy_bits(s);
|
||||
acc = FLINT_MAX(acc, 0);
|
||||
acc = FLINT_MIN(acc, prec);
|
||||
prec = FLINT_MIN(prec, acc + 20);
|
||||
|
||||
acb_init(t);
|
||||
mag_init(rad);
|
||||
|
|
@ -101,18 +107,16 @@ acb_dirichlet_zeta_rs(acb_t res, const acb_t s, slong K, slong prec)
|
|||
/* rad = rad(s) */
|
||||
mag_hypot(rad, arb_radref(acb_realref(s)), arb_radref(acb_imagref(s)));
|
||||
|
||||
/* R >= rad(s) */
|
||||
mag_set_d(R, 0.125);
|
||||
mag_max(R, R, rad);
|
||||
|
||||
/* mid(t) = mid(s); rad(t) = R */
|
||||
/* s +/- R */
|
||||
mag_set_ui_2exp_si(R, 1, -3);
|
||||
acb_set(t, s);
|
||||
mag_zero(arb_radref(acb_realref(t)));
|
||||
mag_zero(arb_radref(acb_imagref(t)));
|
||||
mag_add(arb_radref(acb_realref(t)), arb_radref(acb_realref(t)), R);
|
||||
mag_add(arb_radref(acb_imagref(t)), arb_radref(acb_imagref(t)), R);
|
||||
|
||||
/* |zeta'(s)| <= |zeta(t)| / R */
|
||||
/* |zeta'(s)| <= |zeta(s +/- R)| / R */
|
||||
acb_dirichlet_zeta_bound(err, t);
|
||||
mag_div(err, err, R);
|
||||
|
||||
/* error <= |zeta'(s)| * rad(s) */
|
||||
mag_mul(err, err, rad);
|
||||
|
||||
|
|
|
|||
|
|
@ -10,6 +10,7 @@
|
|||
*/
|
||||
|
||||
#include "acb_poly.h"
|
||||
#include "acb_dirichlet.h"
|
||||
|
||||
void
|
||||
_acb_poly_zeta_cpx_series(acb_ptr z, const acb_t s, const acb_t a, int deflate, slong d, slong prec)
|
||||
|
|
@ -77,85 +78,6 @@ _acb_poly_zeta_cpx_series(acb_ptr z, const acb_t s, const acb_t a, int deflate,
|
|||
_arb_vec_clear(vb, d);
|
||||
}
|
||||
|
||||
void
|
||||
_acb_poly_zeta_cpx_reflect(acb_ptr t, const acb_t h, const acb_t a, int deflate, slong len, slong prec)
|
||||
{
|
||||
/* use reflection formula */
|
||||
if (arf_sgn(arb_midref(acb_realref(h))) < 0 && acb_is_one(a))
|
||||
{
|
||||
/* zeta(s) = (2*pi)**s * sin(pi*s/2) / pi * gamma(1-s) * zeta(1-s) */
|
||||
acb_t pi, hcopy;
|
||||
acb_ptr f, s1, s2, s3, s4, u;
|
||||
slong i;
|
||||
|
||||
acb_init(pi);
|
||||
acb_init(hcopy);
|
||||
f = _acb_vec_init(2);
|
||||
s1 = _acb_vec_init(len);
|
||||
s2 = _acb_vec_init(len);
|
||||
s3 = _acb_vec_init(len);
|
||||
s4 = _acb_vec_init(len);
|
||||
u = _acb_vec_init(len);
|
||||
acb_set(hcopy, h);
|
||||
|
||||
acb_const_pi(pi, prec);
|
||||
|
||||
/* s1 = (2*pi)**s */
|
||||
acb_mul_2exp_si(pi, pi, 1);
|
||||
_acb_poly_pow_cpx(s1, pi, h, len, prec);
|
||||
acb_mul_2exp_si(pi, pi, -1);
|
||||
|
||||
/* s2 = sin(pi*s/2) / pi */
|
||||
acb_set(f, h);
|
||||
acb_one(f + 1);
|
||||
acb_mul_2exp_si(f, f, -1);
|
||||
acb_mul_2exp_si(f + 1, f + 1, -1);
|
||||
_acb_poly_sin_pi_series(s2, f, 2, len, prec);
|
||||
_acb_vec_scalar_div(s2, s2, len, pi, prec);
|
||||
|
||||
/* s3 = gamma(1-s) */
|
||||
acb_sub_ui(f, hcopy, 1, prec);
|
||||
acb_neg(f, f);
|
||||
acb_set_si(f + 1, -1);
|
||||
_acb_poly_gamma_series(s3, f, 2, len, prec);
|
||||
|
||||
/* s4 = zeta(1-s) */
|
||||
acb_sub_ui(f, hcopy, 1, prec);
|
||||
acb_neg(f, f);
|
||||
_acb_poly_zeta_cpx_series(s4, f, a, 0, len, prec);
|
||||
for (i = 1; i < len; i += 2)
|
||||
acb_neg(s4 + i, s4 + i);
|
||||
|
||||
_acb_poly_mullow(u, s1, len, s2, len, len, prec);
|
||||
_acb_poly_mullow(s1, s3, len, s4, len, len, prec);
|
||||
_acb_poly_mullow(t, u, len, s1, len, len, prec);
|
||||
|
||||
/* add 1/(1-(s+t)) = 1/(1-s) + t/(1-s)^2 + ... */
|
||||
if (deflate)
|
||||
{
|
||||
acb_sub_ui(u, hcopy, 1, prec);
|
||||
acb_neg(u, u);
|
||||
acb_inv(u, u, prec);
|
||||
for (i = 1; i < len; i++)
|
||||
acb_mul(u + i, u + i - 1, u, prec);
|
||||
_acb_vec_add(t, t, u, len, prec);
|
||||
}
|
||||
|
||||
acb_clear(pi);
|
||||
acb_clear(hcopy);
|
||||
_acb_vec_clear(f, 2);
|
||||
_acb_vec_clear(s1, len);
|
||||
_acb_vec_clear(s2, len);
|
||||
_acb_vec_clear(s3, len);
|
||||
_acb_vec_clear(s4, len);
|
||||
_acb_vec_clear(u, len);
|
||||
}
|
||||
else
|
||||
{
|
||||
_acb_poly_zeta_cpx_series(t, h, a, deflate, len, prec);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
_acb_poly_zeta_series(acb_ptr res, acb_srcptr h, slong hlen, const acb_t a, int deflate, slong len, slong prec)
|
||||
{
|
||||
|
|
@ -165,7 +87,10 @@ _acb_poly_zeta_series(acb_ptr res, acb_srcptr h, slong hlen, const acb_t a, int
|
|||
t = _acb_vec_init(len);
|
||||
u = _acb_vec_init(len);
|
||||
|
||||
_acb_poly_zeta_cpx_reflect(t, h, a, deflate, len, prec);
|
||||
if (acb_is_one(a))
|
||||
acb_dirichlet_zeta_jet(t, h, deflate, len, prec);
|
||||
else
|
||||
_acb_poly_zeta_cpx_series(t, h, a, deflate, len, prec);
|
||||
|
||||
/* compose with nonconstant part */
|
||||
acb_zero(u);
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue