mirror of
https://github.com/vale981/arb
synced 2025-03-06 09:51:39 -05:00
366 lines
14 KiB
ReStructuredText
366 lines
14 KiB
ReStructuredText
**zeta.h** -- support for the zeta function
|
|
===============================================================================
|
|
|
|
Integer zeta values
|
|
-------------------------------------------------------------------------------
|
|
|
|
.. function:: void fmprb_const_zeta3_bsplit(fmprb_t x, long prec)
|
|
|
|
Sets *x* to Apery's constant `\zeta(3)`, computed by applying binary
|
|
splitting to a hypergeometric series.
|
|
|
|
.. function:: void fmprb_zeta_ui_asymp(fmprb_t z, ulong s, long prec)
|
|
|
|
Assuming `s \ge 2`, approximates `\zeta(s)` by `1 + 2^{-s}` along with
|
|
a correct error bound. We use the following bounds: for `s > b`,
|
|
`\zeta(s) - 1 < 2^{-b}`, and generally,
|
|
`\zeta(s) - (1 + 2^{-s}) < 2^{2-\lfloor 3 s/2 \rfloor}`.
|
|
|
|
.. function:: void fmprb_zeta_ui_euler_product(fmprb_t z, ulong s, long prec)
|
|
|
|
Computes `\zeta(s)` using the Euler product. This is fast only if *s*
|
|
is large compared to the precision.
|
|
|
|
Writing `P(a,b) = \prod_{a \le p \le b} (1 - p^{-s})`, we have
|
|
`1/\zeta(s) = P(a,M) P(M+1,\infty)`.
|
|
|
|
To bound the error caused by truncating
|
|
the product at `M`, we write `P(M+1,\infty) = 1 - \epsilon(s,M)`.
|
|
Since `0 < P(a,M) \le 1`, the absolute error for `\zeta(s)` is
|
|
bounded by `\epsilon(s,M)`.
|
|
|
|
According to the analysis in [Fil1992]_, it holds for all `s \ge 6` and `M \ge 1`
|
|
that `1/P(M+1,\infty) - 1 \le f(s,M) \equiv 2 M^{1-s} / (s/2 - 1)`.
|
|
Thus, we have `1/(1-\epsilon(s,M)) - 1 \le f(s,M)`, and expanding
|
|
the geometric series allows us to conclude that
|
|
`\epsilon(M) \le f(s,M)`.
|
|
|
|
.. function:: void fmprb_zeta_ui_bernoulli(fmprb_t x, ulong n, long prec)
|
|
|
|
Computes `\zeta(n)` for even *n* via the corresponding Bernoulli number.
|
|
|
|
.. function:: void fmprb_zeta_ui_vec_borwein(fmprb_struct * z, ulong start, long num, ulong step, long prec)
|
|
|
|
Evaluates `\zeta(s)` at `\mathrm{num}` consecutive integers *s* beginning
|
|
with *start* and proceeding in increments of *step*.
|
|
Uses Borwein's formula ([Bor2000]_, [GS2003]_),
|
|
implemented to support fast multi-evaluation
|
|
(but also works well for a single *s*).
|
|
|
|
Requires `\mathrm{start} \ge 2`. For efficiency, the largest *s*
|
|
should be at most about as
|
|
large as *prec*. Arguments approaching *LONG_MAX* will cause
|
|
overflows.
|
|
One should therefore only use this function for *s* up to about *prec*, and
|
|
then switch to the Euler product.
|
|
|
|
The algorithm for single *s* is basically identical to the one used in MPFR
|
|
(see [MPFR2012]_ for a detailed description).
|
|
In particular, we evaluate the sum backwards to avoid storing more than one
|
|
`d_k` coefficient, and use integer arithmetic throughout since it
|
|
is convenient and the terms turn out to be slightly larger than
|
|
`2^\mathrm{prec}`.
|
|
The only numerical error in the main loop comes from the division by `k^s`,
|
|
which adds less than 1 unit of error per term.
|
|
For fast multi-evaluation, we repeatedly divide by `k^{\mathrm{step}}`.
|
|
Each division reduces the input error and adds at most 1 unit of
|
|
additional rounding error, so by induction, the error per term
|
|
is always smaller than 2 units.
|
|
|
|
.. function:: void fmprb_zeta_ui_bsplit(fmprb_t x, ulong s, long prec)
|
|
|
|
Computes `\zeta(s)` for arbitrary `s \ge 2` using a binary splitting
|
|
implementation of Borwein's algorithm. This has quasilinear complexity
|
|
with respect to the precision (assuming that `s` is fixed).
|
|
We have
|
|
|
|
.. math ::
|
|
|
|
\zeta(s) = \frac{1}{d_n (1-2^{1-s})}
|
|
\sum_{k=0}^{n-1} \frac{(-1)^k(d_n-d_k)}{(k+1)^s} + \gamma_n(s)
|
|
|
|
where
|
|
|
|
.. math ::
|
|
|
|
d_k = n \sum_{i=0}^k \frac{(n+i-1)! 4^i}{(n-i)! (2i)!}.
|
|
|
|
On the domain of interest, `|\gamma_n(s)| \le 3 / (3 + \sqrt 8)^n`.
|
|
|
|
We write the summation as a system of first-order recurrences for
|
|
`(s_k, d_k, t_k)` where `t_k = d_k - d_{k-1}`. This system is
|
|
described by the matrix equation
|
|
|
|
.. math ::
|
|
|
|
\begin{pmatrix} s_{k+1} \\ d_{k+2} \\ t_{k+3} \end{pmatrix}
|
|
=
|
|
\begin{pmatrix}
|
|
1 & (-1)^k (k+1)^{-s} & 0 \\
|
|
0 & 1 & 1 \\
|
|
0 & 0 & u(k)
|
|
\end{pmatrix}
|
|
\begin{pmatrix} s_k \\ d_{k+1} \\ t_{k+2} \end{pmatrix}.
|
|
|
|
We derive the binary splitting scheme by considering a product
|
|
of an arbitrary pair in the chain `M_{n-1} M_{n-2} \cdots M_1 M_0`.
|
|
This gives
|
|
|
|
.. math ::
|
|
|
|
\begin{pmatrix}
|
|
1 & A_L & B_L \\
|
|
0 & 1 & C_L \\
|
|
0 & 0 & D_L
|
|
\end{pmatrix}
|
|
\begin{pmatrix}
|
|
1 & A_R & B_R \\
|
|
0 & 1 & C_R \\
|
|
0 & 0 & D_R
|
|
\end{pmatrix} =
|
|
\begin{pmatrix}
|
|
1 & A_L+A_R & B_R+A_L C_R+B_L D_R \\
|
|
0 & 1 & C_R+C_L D_R \\
|
|
0 & 0 & D_L D_R
|
|
\end{pmatrix}.
|
|
|
|
The next step is to clear denominators. Instead of putting the
|
|
whole matrix on a common denominator, we optimize by putting `C, D` on a
|
|
denominator `Q_1` (the product of denominators of `u`) and `A, B` on
|
|
a common denominator `Q_3 = Q_1 Q_2` (where `Q_2` is the product of
|
|
`(k+1)^s` factors). This gives a small efficiency improvement. Thus,
|
|
we have
|
|
|
|
.. math ::
|
|
|
|
\begin{pmatrix}
|
|
1 & \dfrac{A_L}{Q_{3L}} & \dfrac{B_L}{Q_{3L}} \\[3ex]
|
|
0 & 1 & \dfrac{C_L}{Q_{1L}} \\[3ex]
|
|
0 & 0 & \dfrac{D_L}{Q_{1L}}
|
|
\end{pmatrix}
|
|
\begin{pmatrix}
|
|
1 & \dfrac{A_R}{Q_{3R}} & \dfrac{B_R}{Q_{3R}} \\[3ex]
|
|
0 & 1 & \dfrac{C_R}{Q_{1R}} \\[3ex]
|
|
0 & 0 & \dfrac{D_R}{Q_{1R}}
|
|
\end{pmatrix} =
|
|
\begin{pmatrix}
|
|
1 & \dfrac{Q_{3L} A_R + A_L Q_{3R}}{Q_{3L} Q_{3R}} & \dfrac{Q_{3L} B_R + A_L C_R Q_{2R} + B_L D_R Q_{2R}}{Q_{3L} Q_{3R}} \\[3ex]
|
|
0 & 1 & \dfrac{Q_{1L} C_R + C_L D_R}{Q_{1L} Q_{1R}} \\[3ex]
|
|
0 & 0 & \dfrac{D_L D_R}{Q_{1L} Q_{1R}}
|
|
\end{pmatrix}.
|
|
|
|
In the final matrix, we note that
|
|
`A / Q_3 = \sum_k (-1)^k (k+1)^{-s}`, and `C / Q_1 = d_n`.
|
|
Thus `(1 / d_n) \sum_k (-1)^k (k+1)^{-s} (d_n - d_k)` is given by
|
|
`A/Q_3 - (B/Q_3) / (C/Q_1) = (A C - B Q_1) / (Q_3 C)`.
|
|
|
|
.. function:: void fmprb_zeta_ui(fmprb_t x, ulong s, long prec)
|
|
|
|
Computes `\zeta(s)` for nonnegative integer `s \ne 1`, automatically
|
|
choosing an appropriate algorithm.
|
|
|
|
.. function:: void fmprb_zeta_ui_vec(fmprb_struct * x, ulong start, long num, long prec)
|
|
|
|
.. function:: void fmprb_zeta_ui_vec_even(fmprb_struct * x, ulong start, long num, long prec)
|
|
|
|
.. function:: void fmprb_zeta_ui_vec_odd(fmprb_struct * x, ulong start, long num, long prec)
|
|
|
|
Computes `\zeta(s)` at num consecutive integers (respectively num
|
|
even or num odd integers) beginning with `s = \mathrm{start} \ge 2`,
|
|
automatically choosing an appropriate algorithm.
|
|
|
|
|
|
Related constants
|
|
-------------------------------------------------------------------------------
|
|
|
|
.. function:: void fmprb_const_khinchin(fmprb_t res, long prec)
|
|
|
|
Sets *res* to Khinchin's constant `K_0`, computed as
|
|
|
|
.. math ::
|
|
|
|
\log K_0 = \frac{1}{\log 2} \left[
|
|
\sum_{k=2}^{N-1} \log \left(\frac{k-1}{k} \right) \log \left(\frac{k+1}{k} \right)
|
|
+ \sum_{n=1}^\infty
|
|
\frac {\zeta (2n,N)}{n} \sum_{k=1}^{2n-1} \frac{(-1)^{k+1}}{k}
|
|
\right]
|
|
|
|
where `N \ge 2` is a free parameter that can be used for tuning [BBC1997]_.
|
|
If the infinite series is truncated after `n = M`, the remainder
|
|
is smaller in absolute value than
|
|
|
|
.. math ::
|
|
|
|
\sum_{n=M+1}^{\infty} \zeta(2n, N) =
|
|
\sum_{n=M+1}^{\infty} \sum_{k=0}^{\infty} (k+N)^{-2n} \le
|
|
\sum_{n=M+1}^{\infty} \left( N^{-2n} + \int_0^{\infty} (t+N)^{-2n} dt \right)
|
|
|
|
= \sum_{n=M+1}^{\infty} \frac{1}{N^{2n}} \left(1 + \frac{N}{2n-1}\right)
|
|
\le \sum_{n=M+1}^{\infty} \frac{N+1}{N^{2n}} = \frac{1}{N^{2M} (N-1)}
|
|
\le \frac{1}{N^{2M}}.
|
|
|
|
Thus, for an error of at most `2^{-p}` in the series,
|
|
it is sufficient to choose `M \ge p / (2 \log_2 N)`.
|
|
|
|
|
|
Euler-Maclaurin summation
|
|
-------------------------------------------------------------------------------
|
|
|
|
.. 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`, 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 ::
|
|
|
|
\sum_{k=N}^U f(k) = \int_N^U f(t) dt + \frac{1}{2} \left(f(N) + f(U)\right)
|
|
|
|
+ \sum_{k=1}^{M} \frac{B_{2k}}{(2k)!} \left( f^{(2k-1)}(U) - f^{(2k-1)}(N) \right)
|
|
|
|
- \int_N^U \frac{\tilde B_{2M}(t)}{(2M)!} f^{(2M)}(t) dt
|
|
|
|
where `f` is a sufficiently differentiable function (for example,
|
|
analytic), `B_n` is a Bernoulli number, and
|
|
`\tilde B_n(t) = B_n(t-\lfloor t\rfloor)` is a periodic Bernoulli
|
|
polynomial. If `f` decreases sufficiently rapidly, the formula
|
|
remains valid after letting `U \to \infty`.
|
|
|
|
To evaluate the Hurwitz zeta function, we set `f(k) = (a + k)^{-s}`,
|
|
giving `\zeta(s,a) = \sum_{k=0}^{N-1} f(k) + \sum_{k=N}^{\infty} f(k)`,
|
|
where EM summation is applied to the right sum.
|
|
By choosing `M` and `N` large enough, and taking the standard
|
|
logarithm branch cut, the EM formula gives an analytic
|
|
continuation of `\zeta(s,a)` to all `a, s \in \mathbb{C}`
|
|
(except for poles at `s = 1` and
|
|
`\mathrm{Re}(s) > 0, a = 0, -1, -2, \ldots`). In order to
|
|
evaluate derivatives with respect to `s` of `\zeta(s,a)`, we
|
|
substitute `s \to s + x \in \mathbb{C}[[x]]`.
|
|
|
|
We choose `N` such that `\Re(a+N) > 0`. Then the first integral is
|
|
well-defined for `s` with `\Re(s) > 1` and has the closed form
|
|
|
|
.. math ::
|
|
|
|
\int_N^{\infty} f(t) dt = \int_N^{\infty}
|
|
(a + t)^{-s}dt = \frac{(a+N)^{1-s}}{s-1},
|
|
|
|
providing analytic continuation of this term with respect to `s`.
|
|
Removing the singularity from this term also conveniently allows us
|
|
to evaluate derivatives of `\zeta(s,a) - 1/(s-1)` at `s = 1`.
|
|
|
|
The derivatives of `f(k)` are given by
|
|
|
|
.. math ::
|
|
|
|
f^{(r)}(k) = \frac{(-1)^r (s)_{r}}{(a+k)^{s+r}}
|
|
|
|
where `(s)_{r} = s (s+1) \cdots (s+r-1)` denotes a rising factorial.
|
|
Thus, the remainder integral becomes
|
|
|
|
.. math ::
|
|
|
|
R(s) = \int_N^{\infty} \frac{\tilde B_{2M}(t)}{(2M)!} \frac{(s)_{2M}}{(a+t)^{s+2M}} dt,
|
|
|
|
valid when `\Re(a+N) > 0` and `\Re(s+2M-1) > 0`. We will use the
|
|
stronger condition `\Re(a+N) > 1`.
|
|
|
|
If `F = \sum_k f_k x^k \in \mathbb{C}[[x]]`, define
|
|
`|F| = \sum_k |f_k| x^k` and `|F| \le |G|` if
|
|
`\forall_k : |f_k| \le |g_k|`. It is easy to check that
|
|
`|F + G| \le |F| + |G|` and `|FG| \le |F||G|`. With this notation,
|
|
|
|
.. math ::
|
|
|
|
|R(s+x)| = \left|\int_N^{\infty} \frac{\tilde B_{2M}(t)}{(2M)!}
|
|
\frac{(s+x)_{2M}}{(a+t)^{s+x+2M}} dt\right|
|
|
|
|
\le \int_N^{\infty} \left| \frac{\tilde B_{2M}(t)}{(2M)!}
|
|
\frac{(s+x)_{2M}}{(a+t)^{s+x+2M}} \right| dt
|
|
|
|
\le \frac{4 \left| (s+x)_{2M} \right|}{(2 \pi)^{2M}}
|
|
\int_N^{\infty} \left| \frac{1}{(a+t)^{s+x+2M}} \right| dt
|
|
|
|
where the fact that `|\tilde B_{2M}(x)| < 4 (2M)! / (2\pi)^{2M}` has
|
|
been invoked. Thus it remains to bound the coefficients `R_k` satisfying
|
|
|
|
.. math ::
|
|
|
|
\int_N^{\infty} \left| \frac{1}{(a+t)^{s+x+2M}} \right| dt =
|
|
\sum_k R_k x^k, \quad
|
|
R_k = \int_N^{\infty} \frac{1}{k!}
|
|
\left| \frac{\log(a+t)^k}{(a+t)^{s+2M}} \right| dt.
|
|
|
|
Writing `a = \alpha + \beta i`, where by assumption
|
|
`\alpha + t \ge \alpha + N \ge 1`, we have
|
|
|
|
.. math ::
|
|
|
|
|\log(\alpha + \beta i + t)|
|
|
= \left|\log(\alpha + t) +
|
|
\log\left( 1 + \frac{\beta i}{\alpha + t}\right) \right|
|
|
\le \log(\alpha + t) + \left|\log\left(1 + \frac{\beta i}{\alpha+t}\right)\right|
|
|
|
|
= \log(\alpha+t) + \left|\frac{1}{2}\log\left(1+\frac{\beta^2}{(\alpha+t)^2}\right)
|
|
+ i\tan^{-1}\left(\frac{\beta}{\alpha + t}\right)\right| \le \log(\alpha + t) + C
|
|
|
|
where
|
|
|
|
.. math ::
|
|
|
|
C = \frac{1}{2}\log\left(1+\frac{\beta^2}{(\alpha+N)^2}\right) +
|
|
\tan^{-1}\left(\frac{|\beta|}{\alpha+N}\right) \le \frac{\beta^2}{2 (\alpha+N)^2}
|
|
+ \frac{|\beta|}{(\alpha+N)}.
|
|
|
|
Also writing `s = \sigma + \tau i`, where by assumption `\sigma + 2M > 1`,
|
|
we have
|
|
|
|
.. math ::
|
|
|
|
\frac{1}{|(\alpha+\beta i+t)^{\sigma+\tau i + 2M}|}
|
|
= \frac{e^{\tau \operatorname{arg}(\alpha+\beta i+t)}}{|\alpha+\beta i+t|^{\sigma+2M}}
|
|
\le \frac{K}{(\alpha + t)^{\sigma+2M}}
|
|
|
|
where `K = \exp(\max(0, \tau \tan^{-1}(\beta / (\alpha + N))))`. Finally,
|
|
|
|
.. math ::
|
|
|
|
R_k \le \frac{K}{k!} \, I_k(N+\alpha, \sigma + 2M, C)
|
|
|
|
with the `K` and `C` defined above, where `I_k(A,B,C)` denotes the
|
|
sequence of integrals
|
|
|
|
.. math ::
|
|
|
|
I_k(A,B,C) \equiv \int_A^{\infty} t^{-B} (C + \log t)^k dt
|
|
|
|
which can be evaluated as
|
|
|
|
.. math ::
|
|
|
|
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)`.
|
|
|
|
|
|
|