mirror of
https://github.com/vale981/arb
synced 2025-03-05 17:31:38 -05:00
128 lines
5.2 KiB
ReStructuredText
128 lines
5.2 KiB
ReStructuredText
.. _bernoulli:
|
|
|
|
**bernoulli.h** -- support for Bernoulli numbers
|
|
===============================================================================
|
|
|
|
This module provides helper functions for exact or approximate calculation
|
|
of the Bernoulli numbers, which are defined by the exponential
|
|
generating function
|
|
|
|
.. math ::
|
|
|
|
\frac{x}{e^x-1} = \sum_{n=0}^{\infty} B_n \frac{x^n}{n!}.
|
|
|
|
Efficient algorithms are implemented for both multi-evaluation
|
|
and calculation of isolated Bernoulli numbers.
|
|
A global (or thread-local) cache is also provided,
|
|
to support fast repeated evaluation of various special functions
|
|
that depend on the Bernoulli numbers (including the gamma function
|
|
and the Riemann zeta function).
|
|
|
|
|
|
Generation of Bernoulli numbers
|
|
--------------------------------------------------------------------------------
|
|
|
|
.. type:: bernoulli_rev_t
|
|
|
|
An iterator object for generating a range of even-indexed Bernoulli numbers
|
|
exactly in reverse order, i.e. computing the exact
|
|
fractions `B_n, B_{n-2}, B_{n-4}, \ldots, B_0`.
|
|
The Bernoulli numbers are generated from scratch, i.e.
|
|
no caching is performed.
|
|
|
|
The Bernoulli numbers are computed by direct summation of the zeta series.
|
|
This is made fast by storing a table of powers (as done by [Blo2009]_).
|
|
As an optimization, we only include the odd powers, and use
|
|
fixed-point arithmetic.
|
|
|
|
The reverse iteration order is preferred for performance reasons,
|
|
as the powers can be updated using multiplications instead of divisions,
|
|
and we avoid having to periodically recompute terms to higher precision.
|
|
To generate Bernoulli numbers in the forward direction without having
|
|
to store all of them, one can split the desired range into smaller
|
|
blocks and compute each block with a single reverse pass.
|
|
|
|
.. function:: void bernoulli_rev_init(bernoulli_rev_t iter, ulong n)
|
|
|
|
Initializes the iterator *iter*. The first Bernoulli number to
|
|
be generated by calling :func:`bernoulli_rev_next` is `B_n`.
|
|
It is assumed that `n` is even.
|
|
|
|
.. function:: void bernoulli_rev_next(fmpz_t numer, fmpz_t denom, bernoulli_rev_t iter)
|
|
|
|
Sets *numer* and *denom* to the exact, reduced numerator and denominator
|
|
of the Bernoulli number `B_k` and advances the state of *iter*
|
|
so that the next invocation generates `B_{k-2}`.
|
|
|
|
.. function:: void bernoulli_rev_clear(bernoulli_rev_t iter)
|
|
|
|
Frees all memory allocated internally by *iter*.
|
|
|
|
|
|
Caching
|
|
-------------------------------------------------------------------------------
|
|
|
|
.. var:: slong bernoulli_cache_num
|
|
|
|
.. var:: fmpq * bernoulli_cache
|
|
|
|
Cache of Bernoulli numbers. Uses thread-local storage if enabled
|
|
in FLINT.
|
|
|
|
.. function:: void bernoulli_cache_compute(slong n)
|
|
|
|
Makes sure that the Bernoulli numbers up to at least `B_{n-1}` are cached.
|
|
Calling :func:`flint_cleanup()` frees the cache.
|
|
|
|
|
|
Bounding
|
|
-------------------------------------------------------------------------------
|
|
|
|
.. function:: slong bernoulli_bound_2exp_si(ulong n)
|
|
|
|
Returns an integer `b` such that `|B_n| \le 2^b`. Uses a lookup table
|
|
for small `n`, and for larger `n` uses the inequality
|
|
`|B_n| < 4 n! / (2 \pi)^n < 4 (n+1)^{n+1} e^{-n} / (2 \pi)^n`.
|
|
Uses integer arithmetic throughout, with the bound for the logarithm
|
|
being looked up from a table. If `|B_n| = 0`, returns *LONG_MIN*.
|
|
Otherwise, the returned exponent `b` is never more than one percent
|
|
larger than the true magnitude.
|
|
|
|
This function is intended for use when `n` small enough that one might
|
|
comfortably compute `B_n` exactly. It aborts if `n` is so large that
|
|
internal overflow occurs.
|
|
|
|
Isolated Bernoulli numbers
|
|
-------------------------------------------------------------------------------
|
|
|
|
.. function:: ulong bernoulli_mod_p_harvey(ulong n, ulong p)
|
|
|
|
Returns the `B_n` modulo the prime number *p*, computed using
|
|
Harvey's algorithm [Har2010]_. The running time is linear in *p*.
|
|
If *p* divides the numerator of `B_n`, *UWORD_MAX* is returned
|
|
as an error code.
|
|
|
|
.. function:: void _bernoulli_fmpq_ui_zeta(fmpz_t num, fmpz_t den, ulong n)
|
|
void _bernoulli_fmpq_ui_multi_mod(fmpz_t num, fmpz_t den, ulong n, double alpha)
|
|
|
|
Sets *num* and *den* to the reduced numerator and denominator
|
|
of the Bernoulli number `B_n`.
|
|
|
|
The *zeta* version computes the denominator `d` using the von Staudt-Clausen
|
|
theorem, numerically approximates `B_n` using :func:`arb_bernoulli_ui_zeta`,
|
|
and then rounds `d B_n` to the correct numerator.
|
|
|
|
The *multi_mod* version reconstructs `B_n` by computing the high bits
|
|
via the Riemann zeta function and the low bits via Harvey's multimodular
|
|
algorithm. The tuning parameter *alpha* should be a fraction between
|
|
0 and 1 controlling the number of bits to compute by the multimodular
|
|
algorithm. If set to a negative number, a default value will be used.
|
|
|
|
.. function:: void _bernoulli_fmpq_ui(fmpz_t num, fmpz_t den, ulong n)
|
|
void bernoulli_fmpq_ui(fmpq_t b, ulong n)
|
|
|
|
Computes the Bernoulli number `B_n` as an exact fraction, for an
|
|
isolated integer `n`. This function reads `B_n` from the global cache
|
|
if the number is already cached, but does not automatically extend
|
|
the cache by itself.
|
|
|