mirror of
https://github.com/vale981/arb
synced 2025-03-05 17:31:38 -05:00
slight improvements to zeta code; remove old implementation
This commit is contained in:
Fredrik Johansson
parent
3f0901e71a
commit
56ea7b2e3b
8 changed files with 130 additions and 438 deletions
|
|
@ -306,67 +306,22 @@ Elementary functions
|
|||
Special functions
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
.. function:: void fmpcb_zeta_em_bound(fmpr_t bound, const fmpcb_t s, ulong N, ulong M, long prec)
|
||||
|
||||
Sets *bound* to a bound for the absolute truncation error when evaluating
|
||||
`\zeta(s)` using the Euler-Maclaurin formula and truncating the
|
||||
main series before term `N` and the tail series before term `M`.
|
||||
If `s = a + bi`, then provided `a > -2M+1`, we use Backlund's inequality
|
||||
(see [GS2003]_),
|
||||
|
||||
.. math ::
|
||||
|
||||
|R| < \left| \frac{s(s+1)\cdots(s+2M-1) B_{2M}}
|
||||
{(2M)! (a+2M-1) N^{s+2M-1}} \right|
|
||||
|
||||
Since the bound usually does not need to be as tight as possible,
|
||||
we evaluate some of the terms in the formula using order-of-magnitude
|
||||
bounds. The calculations are done using a working precision of *prec*
|
||||
bits.
|
||||
|
||||
.. function:: void fmpcb_zeta_em_choose_param(fmpr_t bound, ulong * N, ulong * M, const fmpcb_t s, long target, long prec)
|
||||
|
||||
Finds parameters *N* and *M* such that the truncation error
|
||||
for Euler-Maclaurin summation of `\zeta(s)` is smaller than
|
||||
`2^{-t}` where `t` is given by *target*, and sets *bound* to the
|
||||
corresponding error bound.
|
||||
The *prec* parameter specified the working precision used for
|
||||
the error bounding itself.
|
||||
|
||||
We use bisection on *N*, choosing `M(N)` heuristically. This strategy
|
||||
can certainly be improved.
|
||||
|
||||
.. function:: void fmpcb_zeta_em_sum(fmpcb_t z, const fmpcb_t s, ulong N, ulong M, long prec)
|
||||
|
||||
Evaluates the truncated Euler-Maclaurin sum for the Riemann zeta function
|
||||
|
||||
.. math ::
|
||||
|
||||
\zeta(s) - \epsilon = \sum_{n=1}^{N-1} \frac{1}{n^s} + \frac{N^{1-s}}{s-1}
|
||||
+ \frac{N^{-s}}{2} +
|
||||
\sum_{r=1}^{M-1} \frac{B_{2r} s (s+1) \cdots (s+2r-2)}{(2r)! N^{s+2r-1}}
|
||||
|
||||
using a working precision of *prec* bits.
|
||||
|
||||
.. function:: void fmpcb_zeta(fmpcb_t z, const fmpcb_t s, long prec)
|
||||
|
||||
Sets *z* to the value of the Riemann zeta function `\zeta(s)`.
|
||||
Uses Euler-Maclaurin summation with a working precision of *prec* bits and
|
||||
default parameters obtained from *fmpcb_zeta_em_choose_param*,
|
||||
targeting an absolute truncation error of `2^{-\operatorname{prec}}`.
|
||||
|
||||
.. function:: 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)
|
||||
|
||||
.. function:: void fmpcb_zeta_series(fmpcb_struct * z, const fmpcb_t s, const fmpcb_t a, int deflate, long d, long prec)
|
||||
|
||||
Evaluates the truncated Euler-Maclaurin sum of order `N, M` for the
|
||||
length-*d* truncated Taylor series of the Hurwitz zeta function
|
||||
`\zeta(s,a)` at `s`. With `a = 1`, this gives the usual
|
||||
Riemann zeta function.
|
||||
`\zeta(s,a)` at `s`, using a working precision of *prec* bits.
|
||||
With `a = 1`, this gives the usual Riemann zeta function.
|
||||
|
||||
If *deflate* is nonzero, `\zeta(s,a) - 1/(s-1)` is evaluated
|
||||
(which permits series expansion at `s = 1`).
|
||||
|
||||
The *fmpcb_zeta_series* function chooses default values for `N, M`
|
||||
using *fmpcb_zeta_series_em_choose_param*,
|
||||
targeting an absolute truncation error of `2^{-\operatorname{prec}}`.
|
||||
|
||||
The Euler-Maclaurin (EM) formula states that
|
||||
|
||||
.. math ::
|
||||
|
|
@ -498,3 +453,12 @@ Special functions
|
|||
I_k(A,B,C) = \frac{L_k}{(B-1)^{k+1} A^{B-1}}
|
||||
|
||||
where `L_0 = 1`, `L_k = k L_{k-1} + D^k` and `D = (B-1) (C + \log A)`.
|
||||
|
||||
.. function:: void fmpcb_zeta_series_em_choose_param(fmpr_t bound, ulong * N, ulong * M, const fmpcb_t s, const fmpcb_t a, long d, long target, long prec)
|
||||
|
||||
Chooses *N* and *M* using a default algorithm.
|
||||
|
||||
.. function:: void fmpcb_zeta(fmpcb_t z, const fmpcb_t s, long prec)
|
||||
|
||||
Sets *z* to the value of the Riemann zeta function `\zeta(s)`.
|
||||
|
||||
|
|
|
|||
13
fmpcb.h
13
fmpcb.h
|
|
@ -440,23 +440,12 @@ void fmpcb_root_exp(fmpcb_t r, const fmpcb_t a, long m, long index, long prec);
|
|||
void fmpcb_root_newton(fmpcb_t r, const fmpcb_t a, long m, long index, long prec);
|
||||
void fmpcb_root(fmpcb_t r, const fmpcb_t a, long m, long index, long prec);
|
||||
|
||||
void fmpcb_zeta_em_bound(fmpr_t bound, const fmpcb_t s, ulong N, ulong M, long prec);
|
||||
|
||||
void fmpcb_zeta_em_choose_param(fmpr_t bound, ulong * N, ulong * M, const fmpcb_t s, long target, long prec);
|
||||
|
||||
void fmpcb_zeta_em_sum(fmpcb_t z, const fmpcb_t s, ulong N, ulong M, long prec);
|
||||
|
||||
void fmpcb_zeta(fmpcb_t z, const fmpcb_t s, long 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);
|
||||
void fmpcb_zeta_series_em_choose_param(fmpr_t bound, ulong * N, ulong * M, const fmpcb_t s, const fmpcb_t a, long d, long target, long prec);
|
||||
|
||||
|
||||
void fmpcb_zeta_series_em_bound(fmpr_t bound, const fmpcb_t s, const fmpcb_t a, long N, long M, long d, long wp);
|
||||
void fmpcb_zeta_series_em_vec_bound(fmprb_struct * vec, const fmpcb_t s, const fmpcb_t a, ulong N, ulong M, long d, long wp);
|
||||
void fmpcb_zeta_series(fmpcb_struct * z, const fmpcb_t s, const fmpcb_t a, int deflate, long d, long prec);
|
||||
void fmpcb_zeta(fmpcb_t z, const fmpcb_t s, long prec);
|
||||
|
||||
|
||||
static __inline__ void
|
||||
|
|
|
|||
18
fmpcb/zeta.c
18
fmpcb/zeta.c
|
|
@ -28,16 +28,10 @@
|
|||
void
|
||||
fmpcb_zeta(fmpcb_t z, const fmpcb_t s, long prec)
|
||||
{
|
||||
ulong M, N;
|
||||
fmpr_t bound;
|
||||
|
||||
fmpr_init(bound);
|
||||
fmpcb_zeta_em_choose_param(bound, &N, &M, s, prec, FMPRB_RAD_PREC);
|
||||
|
||||
fmpcb_zeta_em_sum(z, s, N, M, prec);
|
||||
|
||||
fmprb_add_error_fmpr(fmpcb_realref(z), bound);
|
||||
fmprb_add_error_fmpr(fmpcb_imagref(z), bound);
|
||||
|
||||
fmpr_clear(bound);
|
||||
fmpcb_t a;
|
||||
fmpcb_init(a);
|
||||
fmpcb_one(a);
|
||||
fmpcb_zeta_series(z, s, a, 0, 1, prec);
|
||||
fmpcb_clear(a);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -1,141 +0,0 @@
|
|||
/*=============================================================================
|
||||
|
||||
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 "fmpcb.h"
|
||||
|
||||
/* |B_n| / n! < 4 * (2 pi)^(-n) = 2^(2-log_2(2pi)*n) */
|
||||
void
|
||||
fmpr_bernoulli_scaled_ui_ubound(fmpr_t b, ulong n, long prec)
|
||||
{
|
||||
fmpr_set_si_2exp_si(b, 1, 2 - 265*n/100);
|
||||
}
|
||||
|
||||
/* Absolute value of rising factorial (could speed up once complex gamma is available). */
|
||||
void
|
||||
fmpcb_rfac_abs_ubound(fmpr_t bound, const fmpcb_t s, ulong n, long prec)
|
||||
{
|
||||
fmpr_t term, t;
|
||||
ulong k;
|
||||
|
||||
/* M(k) = (a+k)^2 + b^2
|
||||
M(0) = a^2 + b^2
|
||||
M(k+1) = M(k) + 2*a + (2*k+1)
|
||||
*/
|
||||
fmpr_init(t);
|
||||
fmpr_init(term);
|
||||
|
||||
fmpr_one(bound);
|
||||
|
||||
/* M(0) = a^2 + b^2 */
|
||||
fmprb_get_abs_ubound_fmpr(t, fmpcb_realref(s), prec);
|
||||
fmpr_mul(term, t, t, prec, FMPR_RND_UP);
|
||||
fmprb_get_abs_ubound_fmpr(t, fmpcb_imagref(s), prec);
|
||||
fmpr_mul(t, t, t, prec, FMPR_RND_UP);
|
||||
fmpr_add(term, term, t, prec, FMPR_RND_UP);
|
||||
|
||||
/* we add t = 2*a to each term. note that this can be signed;
|
||||
we always want the most positive value */
|
||||
fmpr_add(t, fmprb_midref(fmpcb_realref(s)),
|
||||
fmprb_radref(fmpcb_realref(s)), prec, FMPR_RND_CEIL);
|
||||
fmpr_mul_2exp_si(t, t, 1);
|
||||
|
||||
for (k = 0; k < n; k++)
|
||||
{
|
||||
fmpr_mul(bound, bound, term, prec, FMPR_RND_UP);
|
||||
fmpr_add_ui(term, term, 2 * k + 1, prec, FMPR_RND_UP);
|
||||
fmpr_add(term, term, t, prec, FMPR_RND_UP);
|
||||
}
|
||||
|
||||
fmpr_sqrt(bound, bound, prec, FMPR_RND_UP);
|
||||
|
||||
fmpr_clear(t);
|
||||
fmpr_clear(term);
|
||||
}
|
||||
|
||||
void
|
||||
_fmpr_get_floor_fmpz(fmpz_t f, const fmpr_t x)
|
||||
{
|
||||
if (COEFF_IS_MPZ(*fmpr_expref(x)))
|
||||
{
|
||||
abort();
|
||||
}
|
||||
|
||||
if (fmpz_sgn(fmpr_expref(x)) >= 0)
|
||||
{
|
||||
fmpz_mul_2exp(f, fmpr_manref(x), *fmpr_expref(x));
|
||||
}
|
||||
else
|
||||
{
|
||||
fmpz_neg(f, fmpr_expref(x));
|
||||
fmpz_fdiv_q_2exp(f, fmpr_manref(x), *f);
|
||||
}
|
||||
}
|
||||
|
||||
/* the bound is rf(s,2M) * (B_{2M} / (2M)!) / ((a+2M-1) * N^(a+2M-1))
|
||||
we assume N > 0, Q > 0 */
|
||||
void
|
||||
fmpcb_zeta_em_bound(fmpr_t bound, const fmpcb_t s, ulong N, ulong M, long prec)
|
||||
{
|
||||
fmpr_t delta, t, u;
|
||||
|
||||
fmpr_init(delta);
|
||||
|
||||
/* delta = a + 2*q - 1, rounded down (we require it to be positive) */
|
||||
fmpr_sub(delta, fmprb_midref(fmpcb_realref(s)),
|
||||
fmprb_radref(fmpcb_realref(s)), prec, FMPR_RND_FLOOR);
|
||||
fmpr_add_ui(delta, delta, 2*M - 1, prec, FMPR_RND_FLOOR);
|
||||
|
||||
if (fmpr_sgn(delta) <= 0 || N < 1 || M < 1)
|
||||
{
|
||||
fmpr_clear(delta);
|
||||
fmpr_pos_inf(bound);
|
||||
return;
|
||||
}
|
||||
|
||||
fmpr_init(t);
|
||||
fmpr_init(u);
|
||||
|
||||
/* compute numerator (rounding up) */
|
||||
fmpcb_rfac_abs_ubound(t, s, 2 * M, prec);
|
||||
fmpr_bernoulli_scaled_ui_ubound(u, 2*M, prec);
|
||||
fmpr_mul(t, t, u, prec, FMPR_RND_UP);
|
||||
|
||||
/* compute denominator (rounding down) */
|
||||
fmpr_set_ui(u, N);
|
||||
{
|
||||
fmpz_t f;
|
||||
fmpz_init(f);
|
||||
_fmpr_get_floor_fmpz(f, delta);
|
||||
fmpr_pow_sloppy_fmpz(u, u, f, prec, FMPR_RND_DOWN);
|
||||
fmpz_clear(f);
|
||||
}
|
||||
|
||||
fmpr_mul(u, u, delta, prec, FMPR_RND_DOWN);
|
||||
fmpr_div(bound, t, u, prec, FMPR_RND_UP);
|
||||
|
||||
fmpr_clear(delta);
|
||||
fmpr_clear(t);
|
||||
fmpr_clear(u);
|
||||
}
|
||||
|
|
@ -1,87 +0,0 @@
|
|||
/*=============================================================================
|
||||
|
||||
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 "fmpcb.h"
|
||||
|
||||
static ulong choose_M(ulong N, long target)
|
||||
{
|
||||
return FLINT_MIN(N, target);
|
||||
}
|
||||
|
||||
void
|
||||
fmpcb_zeta_em_choose_param(fmpr_t bound, ulong * N, ulong * M, const fmpcb_t s, long target, long prec)
|
||||
{
|
||||
ulong A, B, C;
|
||||
fmpr_t Abound, Bbound, Cbound, tol;
|
||||
|
||||
fmpr_init(Abound);
|
||||
fmpr_init(Bbound);
|
||||
fmpr_init(Cbound);
|
||||
fmpr_init(tol);
|
||||
|
||||
fmpr_set_si_2exp_si(tol, 1, -target);
|
||||
|
||||
A = 1;
|
||||
B = 2;
|
||||
|
||||
fmpcb_zeta_em_bound(Bbound, s, B, choose_M(B, target), prec);
|
||||
|
||||
if (fmpr_cmp(Bbound, tol) > 0)
|
||||
{
|
||||
while (fmpr_cmp(Bbound, tol) > 0)
|
||||
{
|
||||
fmpr_set(Abound, Bbound);
|
||||
A *= 2;
|
||||
B *= 2;
|
||||
fmpcb_zeta_em_bound(Bbound, s, B, choose_M(B, target), prec);
|
||||
}
|
||||
|
||||
/* bisect (-A, B] */
|
||||
while (B > A + 4)
|
||||
{
|
||||
C = A + (B - A) / 2;
|
||||
fmpcb_zeta_em_bound(Cbound, s, C, choose_M(C, target), prec);
|
||||
if (fmpr_cmp(Cbound, tol) < 0)
|
||||
{
|
||||
B = C;
|
||||
fmpr_set(Bbound, Cbound);
|
||||
}
|
||||
else
|
||||
{
|
||||
A = C;
|
||||
fmpr_set(Abound, Cbound);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fmpr_set(bound, Bbound);
|
||||
*N = B;
|
||||
*M = choose_M(B, target);
|
||||
|
||||
fmpr_clear(Abound);
|
||||
fmpr_clear(Bbound);
|
||||
fmpr_clear(Cbound);
|
||||
fmpr_clear(tol);
|
||||
}
|
||||
|
|
@ -1,118 +0,0 @@
|
|||
/*=============================================================================
|
||||
|
||||
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 "fmpcb.h"
|
||||
#include "bernoulli.h"
|
||||
|
||||
void
|
||||
fmpcb_zeta_em_sum(fmpcb_t z, const fmpcb_t s, ulong N, ulong M, long prec)
|
||||
{
|
||||
fmpcb_t t, u, v, negs, term, sum;
|
||||
fmprb_t x;
|
||||
fmpz_t c;
|
||||
ulong r, n;
|
||||
|
||||
bernoulli_cache_compute(2 * M);
|
||||
|
||||
fmpcb_init(t);
|
||||
fmpcb_init(u);
|
||||
fmpcb_init(v);
|
||||
fmpcb_init(negs);
|
||||
fmpcb_init(term);
|
||||
fmpcb_init(sum);
|
||||
|
||||
fmprb_init(x);
|
||||
fmpz_init(c);
|
||||
|
||||
fmpcb_neg(negs, s);
|
||||
fmpcb_zero(sum);
|
||||
|
||||
/* sum 1/n^s */
|
||||
for (n = 1; n < N; n++)
|
||||
{
|
||||
fmpcb_set_ui(t, n);
|
||||
fmpcb_pow(t, t, negs, prec);
|
||||
fmpcb_add(sum, sum, t, prec);
|
||||
}
|
||||
|
||||
/* t = 1/N^s */
|
||||
fmpcb_set_ui(t, N);
|
||||
fmpcb_pow(t, t, negs, prec);
|
||||
|
||||
/* N / (s-1) * (1 / N^s) */
|
||||
fmpcb_set_ui(u, N);
|
||||
fmpcb_sub_ui(v, s, 1, prec);
|
||||
fmpcb_div(u, u, v, prec);
|
||||
fmpcb_mul(u, u, t, prec);
|
||||
fmpcb_add(sum, sum, u, prec);
|
||||
|
||||
/* 1/2 * (1 / N^s) */
|
||||
fmpcb_mul_2exp_si(u, t, -1);
|
||||
fmpcb_add(sum, sum, u, prec);
|
||||
|
||||
/* term = 1/2 * (1 / N^s) * s / N */
|
||||
fmpcb_mul(u, u, s, prec);
|
||||
fmpcb_div_ui(term, u, N, prec);
|
||||
|
||||
for (r = 1; r < M; r++)
|
||||
{
|
||||
/* 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);
|
||||
|
||||
fmprb_mul(fmpcb_realref(u), fmpcb_realref(term), x, prec);
|
||||
fmprb_mul(fmpcb_imagref(u), fmpcb_imagref(term), x, prec);
|
||||
|
||||
fmpcb_add(sum, sum, u, prec);
|
||||
|
||||
/* multiply term by (s+2r-1)(s+2r) / (N^2 * (2*r+1)*(2*r+2)) */
|
||||
fmpcb_set(u, s);
|
||||
fmprb_add_ui(fmpcb_realref(u), fmpcb_realref(s), 2*r-1, prec);
|
||||
fmpcb_mul(term, term, u, prec);
|
||||
fmprb_add_ui(fmpcb_realref(u), fmpcb_realref(u), 1, prec);
|
||||
|
||||
fmpcb_mul(term, term, u, prec);
|
||||
|
||||
fmpz_set_ui(c, N);
|
||||
fmpz_mul_ui(c, c, N);
|
||||
fmpz_mul_ui(c, c, 2*r+1);
|
||||
fmpz_mul_ui(c, c, 2*r+2);
|
||||
fmprb_div_fmpz(fmpcb_realref(term), fmpcb_realref(term), c, prec);
|
||||
fmprb_div_fmpz(fmpcb_imagref(term), fmpcb_imagref(term), c, prec);
|
||||
|
||||
}
|
||||
|
||||
fmpcb_set(z, sum);
|
||||
|
||||
fmpcb_clear(t);
|
||||
fmpcb_clear(u);
|
||||
fmpcb_clear(v);
|
||||
fmpcb_clear(negs);
|
||||
fmpcb_clear(term);
|
||||
fmpcb_clear(sum);
|
||||
|
||||
fmprb_clear(x);
|
||||
fmpz_clear(c);
|
||||
}
|
||||
|
|
@ -121,18 +121,69 @@ bound_K(fmprb_t C, const fmprb_t AN, const fmprb_t B, const fmprb_t T, long wp)
|
|||
}
|
||||
}
|
||||
|
||||
/* Absolute value of rising factorial (could speed up once complex gamma is available). */
|
||||
void
|
||||
fmpcb_rfac_abs_ubound2(fmpr_t bound, const fmpcb_t s, ulong n, long prec)
|
||||
{
|
||||
fmpr_t term, t;
|
||||
ulong k;
|
||||
|
||||
/* M(k) = (a+k)^2 + b^2
|
||||
M(0) = a^2 + b^2
|
||||
M(k+1) = M(k) + 2*a + (2*k+1)
|
||||
*/
|
||||
fmpr_init(t);
|
||||
fmpr_init(term);
|
||||
|
||||
fmpr_one(bound);
|
||||
|
||||
/* M(0) = a^2 + b^2 */
|
||||
fmprb_get_abs_ubound_fmpr(t, fmpcb_realref(s), prec);
|
||||
fmpr_mul(term, t, t, prec, FMPR_RND_UP);
|
||||
fmprb_get_abs_ubound_fmpr(t, fmpcb_imagref(s), prec);
|
||||
fmpr_mul(t, t, t, prec, FMPR_RND_UP);
|
||||
fmpr_add(term, term, t, prec, FMPR_RND_UP);
|
||||
|
||||
/* we add t = 2*a to each term. note that this can be signed;
|
||||
we always want the most positive value */
|
||||
fmpr_add(t, fmprb_midref(fmpcb_realref(s)),
|
||||
fmprb_radref(fmpcb_realref(s)), prec, FMPR_RND_CEIL);
|
||||
fmpr_mul_2exp_si(t, t, 1);
|
||||
|
||||
for (k = 0; k < n; k++)
|
||||
{
|
||||
fmpr_mul(bound, bound, term, prec, FMPR_RND_UP);
|
||||
fmpr_add_ui(term, term, 2 * k + 1, prec, FMPR_RND_UP);
|
||||
fmpr_add(term, term, t, prec, FMPR_RND_UP);
|
||||
}
|
||||
|
||||
fmpr_sqrt(bound, bound, prec, FMPR_RND_UP);
|
||||
|
||||
fmpr_clear(t);
|
||||
fmpr_clear(term);
|
||||
}
|
||||
|
||||
|
||||
void
|
||||
bound_rfac(fmprb_struct * F, const fmpcb_t s, ulong n, long len, long wp)
|
||||
{
|
||||
fmprb_struct sx[2];
|
||||
fmprb_init(sx + 0);
|
||||
fmprb_init(sx + 1);
|
||||
fmpcb_abs(sx + 0, s, wp);
|
||||
fmprb_one(sx + 1);
|
||||
_fmprb_vec_zero(F, len);
|
||||
_fmprb_poly_rfac_series_ui(F, sx, 2, n, len, wp);
|
||||
fmprb_clear(sx + 0);
|
||||
fmprb_clear(sx + 1);
|
||||
if (len == 1)
|
||||
{
|
||||
fmpcb_rfac_abs_ubound2(fmprb_midref(F + 0), s, n, wp);
|
||||
fmpr_zero(fmprb_radref(F + 0));
|
||||
}
|
||||
else
|
||||
{
|
||||
fmprb_struct sx[2];
|
||||
fmprb_init(sx + 0);
|
||||
fmprb_init(sx + 1);
|
||||
fmpcb_abs(sx + 0, s, wp);
|
||||
fmprb_one(sx + 1);
|
||||
_fmprb_vec_zero(F, len);
|
||||
_fmprb_poly_rfac_series_ui(F, sx, 2, n, len, wp);
|
||||
fmprb_clear(sx + 0);
|
||||
fmprb_clear(sx + 1);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
|
|
|
|||
|
|
@ -76,7 +76,6 @@ void _fmpcb_poly_fmpcb_invpow_cpx(fmpcb_struct * res, const fmpcb_t N, const fmp
|
|||
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++)
|
||||
|
|
@ -88,6 +87,31 @@ void _fmpcb_poly_fmpcb_invpow_cpx(fmpcb_struct * res, const fmpcb_t N, const fmp
|
|||
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)
|
||||
{
|
||||
|
|
@ -97,6 +121,7 @@ fmpcb_zeta_series_em_sum(fmpcb_struct * z, const fmpcb_t s, const fmpcb_t a, int
|
|||
fmpz_t c;
|
||||
long i;
|
||||
ulong r, n;
|
||||
int aint;
|
||||
|
||||
bernoulli_cache_compute(2 * M + 1);
|
||||
|
||||
|
|
@ -114,6 +139,8 @@ fmpcb_zeta_series_em_sum(fmpcb_struct * z, const fmpcb_t s, const fmpcb_t a, int
|
|||
/* sum 1/(n+a)^(s+x) */
|
||||
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);
|
||||
|
|
@ -178,12 +205,17 @@ fmpcb_zeta_series_em_sum(fmpcb_struct * z, const fmpcb_t s, const fmpcb_t a, int
|
|||
_fmpcb_poly_mullow_cpx(u, u, d, s, d, prec);
|
||||
_fmpcb_vec_scalar_div(term, u, d, Na, prec);
|
||||
|
||||
/* 1/(N+a)^2 */
|
||||
/* (N+a)^2 or 1/(N+a)^2 */
|
||||
fmpcb_mul(Na, Na, Na, prec);
|
||||
fmpcb_inv(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);
|
||||
|
|
@ -198,11 +230,19 @@ fmpcb_zeta_series_em_sum(fmpcb_struct * z, const fmpcb_t s, const fmpcb_t a, int
|
|||
fmprb_add_ui(fmpcb_realref(splus), fmpcb_realref(splus), 1, prec);
|
||||
_fmpcb_poly_mullow_cpx(term, term, d, splus, d, prec);
|
||||
|
||||
/* TODO: div fmpz when fmpz! */
|
||||
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);
|
||||
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);
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue