hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-06 09:51:39 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions
arb/zeta/zeta_series_em_sum.c

293 lines
8.6 KiB
C
Raw Blame History

/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
******************************************************************************/
#include "zeta.h"
#include "fmpcb.h"
#include "fmpcb_poly.h"
#include "bernoulli.h"
/* res = src * (c + x) */
void _fmpcb_poly_mullow_cpx(fmpcb_struct * res, const fmpcb_struct * src, long len, const fmpcb_t c, long trunc, long prec)
{
long i;
if (len < trunc)
fmpcb_set(res + len, src + len - 1);
for (i = len - 1; i > 0; i--)
{
fmpcb_mul(res + i, src + i, c, prec);
fmpcb_add(res + i, res + i, src + i - 1, prec);
}
fmpcb_mul(res, src, c, prec);
}
/* series of 1/N^(c+x) */
void _fmpcb_poly_ui_invpow_cpx(fmpcb_struct * res, ulong N, const fmpcb_t c, long trunc, long prec)
{
long i;
fmpcb_t logN;
fmpcb_init(logN);
fmprb_log_ui(fmpcb_realref(logN), N, prec);
fmpcb_mul(res + 0, logN, c, prec);
fmpcb_neg(res + 0, res + 0);
fmpcb_exp(res + 0, res + 0, prec);
for (i = 1; i < trunc; i++)
{
fmpcb_mul(res + i, res + i - 1, logN, prec);
fmpcb_div_si(res + i, res + i, -i, prec);
}
fmpcb_clear(logN);
}
void _fmpcb_poly_fmpcb_invpow_cpx(fmpcb_struct * res, const fmpcb_t N, const fmpcb_t c, long trunc, long prec)
{
long i;
fmpcb_t logN;
fmpcb_init(logN);
fmpcb_log(logN, N, prec);
fmpcb_mul(res + 0, logN, c, prec);
fmpcb_neg(res + 0, res + 0);
fmpcb_exp(res + 0, res + 0, prec);
for (i = 1; i < trunc; i++)
{
fmpcb_mul(res + i, res + i - 1, logN, prec);
fmpcb_div_si(res + i, res + i, -i, prec);
}
fmpcb_clear(logN);
}
static __inline__ int
fmprb_is_int(const fmprb_t x)
{
return fmprb_is_zero(x) ||
(fmprb_is_exact(x) &&
fmpz_sgn(fmpr_expref(fmprb_midref(x))) >= 0);
}
static __inline__ int
fmpcb_is_int(const fmpcb_t z)
{
return fmprb_is_zero(fmpcb_imagref(z)) && fmprb_is_int(fmpcb_realref(z));
}
static __inline__ void
_fmpcb_vec_scalar_div_fmprb(fmpcb_struct * res, const fmpcb_struct * vec, long len, const fmprb_t c, long prec)
{
long i;
for (i = 0; i < len; i++)
{
fmprb_div(fmpcb_realref(res + i), fmpcb_realref(vec + i), c, prec);
fmprb_div(fmpcb_imagref(res + i), fmpcb_imagref(vec + i), c, prec);
}
}
void
fmpcb_zeta_series_em_sum(fmpcb_struct * z, const fmpcb_t s, const fmpcb_t a, int deflate, ulong N, ulong M, long d, long prec)
{
fmpcb_struct *t, *u, *v, *term, *sum;
fmpcb_t splus, Na, rec;
fmprb_t x;
fmpz_t c;
long i;
ulong r, n;
int aint;
bernoulli_cache_compute(2 * M + 1);
t = _fmpcb_vec_init(d + 1);
u = _fmpcb_vec_init(d);
v = _fmpcb_vec_init(d);
term = _fmpcb_vec_init(d);
sum = _fmpcb_vec_init(d);
fmpcb_init(splus);
fmpcb_init(Na);
fmpcb_init(rec);
fmprb_init(x);
fmpz_init(c);
/* sum 1/(n+a)^(s+x) */
if (fmpcb_is_one(a) && (d == 1))
{
fmpcb_struct * pows;
long j;
pows = _fmpcb_vec_init(N + 1);
fmpcb_one(pows + 1);
for (i = 2; i <= N; i++)
{
if (fmpcb_is_zero(pows + i))
{
fmprb_log_ui(fmpcb_realref(pows + i), i, prec);
fmprb_zero(fmpcb_imagref(pows + i));
fmpcb_mul(pows + i, pows + i, s, prec);
fmpcb_neg(pows + i, pows + i);
fmpcb_exp(pows + i, pows + i, prec);
}
for (j = 2; j <= i && i * j <= N; j++)
if (fmpcb_is_zero(pows + i * j))
fmpcb_mul(pows + i * j, pows + i, pows + j, prec);
}
for (i = 1; i <= N; i++)
fmpcb_add(sum, sum, pows + i, prec);
_fmpcb_vec_clear(pows, N + 1);
}
else
{
for (n = 0; n < N; n++)
{
/* printf("sum 1: %ld %ld\n", n, N); */
fmpcb_add_ui(Na, a, n, prec);
_fmpcb_poly_fmpcb_invpow_cpx(t, Na, s, d, prec);
_fmpcb_vec_add(sum, sum, t, d, prec);
}
}
/* t = 1/(N+a)^(s+x); we might need one extra term for deflation */
fmpcb_add_ui(Na, a, N, prec);
_fmpcb_poly_fmpcb_invpow_cpx(t, Na, s, d + 1, prec);
/* sum += (N+a) * 1/((s+x)-1) * t */
if (!deflate)
{
/* u = 1/(s+x) has series [1/(s-1), -1/(s-1)^2, 1/(s-1)^3, ...] */
fmpcb_sub_ui(u + 0, s, 1, prec);
fmpcb_inv(u + 0, u + 0, prec);
for (i = 1; i < d; i++)
fmpcb_mul(u + i, u + i - 1, u + 0, prec);
for (i = 1; i < d; i += 2)
fmpcb_neg(u + i, u + i);
_fmpcb_poly_mullow(v, u, d, t, d, d, prec);
_fmpcb_vec_scalar_mul(v, v, d, Na, prec);
_fmpcb_vec_add(sum, sum, v, d, prec);
}
/* sum += ((N+a)^(1-(s+x)) - 1) / ((s+x) - 1) */
else
{
/* at s = 1, this becomes (N*t - 1)/x, i.e. just remove one coeff */
if (fmpcb_is_one(s))
{
for (i = 0; i < d; i++)
fmpcb_mul(u + i, t + i + 1, Na, prec);
_fmpcb_vec_add(sum, sum, u, d, prec);
}
else
{
/* TODO: this is numerically unstable for large derivatives,
and divides by zero if s contains 1. We want a good
way to evaluate the power series ((N+a)^y - 1) / y where y has
nonzero constant term, without doing a division.
How is this best done? */
_fmpcb_vec_scalar_mul(t, t, d, Na, prec);
fmpcb_sub_ui(t + 0, t + 0, 1, prec);
fmpcb_sub_ui(u + 0, s, 1, prec);
fmpcb_inv(u + 0, u + 0, prec);
for (i = 1; i < d; i++)
fmpcb_mul(u + i, u + i - 1, u + 0, prec);
for (i = 1; i < d; i += 2)
fmpcb_neg(u + i, u + i);
_fmpcb_poly_mullow(v, u, d, t, d, d, prec);
_fmpcb_vec_add(sum, sum, v, d, prec);
_fmpcb_poly_fmpcb_invpow_cpx(t, Na, s, d, prec);
}
}
/* sum += u = 1/2 * t */
_fmpcb_vec_scalar_mul_2exp_si(u, t, d, -1L);
_fmpcb_vec_add(sum, sum, u, d, prec);
/* term = u * (s+x) / N */
_fmpcb_poly_mullow_cpx(u, u, d, s, d, prec);
_fmpcb_vec_scalar_div(term, u, d, Na, prec);
/* (N+a)^2 or 1/(N+a)^2 */
fmpcb_mul(Na, Na, Na, prec);
aint = fmpcb_is_int(Na);
if (!aint)
fmpcb_inv(Na, Na, prec);
for (r = 1; r <= M; r++)
{
/* printf("sum 2: %ld %ld\n", r, M); */
/* sum += bernoulli number * term */
fmprb_set_round_fmpz(x, fmpq_numref(bernoulli_cache + 2 * r), prec);
fmprb_div_fmpz(x, x, fmpq_denref(bernoulli_cache + 2 * r), prec);
_fmpcb_vec_scalar_mul_fmprb(u, term, d, x, prec);
_fmpcb_vec_add(sum, sum, u, d, prec);
/* multiply term by ((s+x)+2r-1)((s+x)+2r) / ((N+a)^2 * (2*r+1)*(2*r+2)) */
fmpcb_set(splus, s);
fmprb_add_ui(fmpcb_realref(splus), fmpcb_realref(splus), 2*r-1, prec);
_fmpcb_poly_mullow_cpx(term, term, d, splus, d, prec);
fmprb_add_ui(fmpcb_realref(splus), fmpcb_realref(splus), 1, prec);
_fmpcb_poly_mullow_cpx(term, term, d, splus, d, prec);
if (aint)
{
fmprb_mul_ui(x, fmpcb_realref(Na), 2*r+1, prec);
fmprb_mul_ui(x, x, 2*r+2, prec);
_fmpcb_vec_scalar_div_fmprb(term, term, d, x, prec);
}
else
{
fmpz_set_ui(c, 2*r+1);
fmpz_mul_ui(c, c, 2*r+2);
fmpcb_div_fmpz(rec, Na, c, prec);
_fmpcb_vec_scalar_mul(term, term, d, rec, prec);
}
}
_fmpcb_vec_set(z, sum, d);
_fmpcb_vec_clear(t, d + 1);
_fmpcb_vec_clear(u, d);
_fmpcb_vec_clear(v, d);
_fmpcb_vec_clear(term, d);
_fmpcb_vec_clear(sum, d);
fmpcb_clear(splus);
fmpcb_clear(Na);
fmpcb_clear(rec);
fmprb_clear(x);
fmpz_clear(c);
}
Reference in a new issue View git blame Copy permalink