arb/doc/source/bernoulli.rst

**bernoulli.h** -- support for Bernoulli numbers
===============================================================================

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 Bloemen et al.
    http://remcobloemen.nl/2009/11/even-faster-zeta-calculation.html).
    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:: long bernoulli_cache_num

.. var:: fmpq * bernoulli_cache

    Global cache of Bernoulli numbers.

.. function:: void bernoulli_cache_compute(long n)

    Makes sure that the Bernoulli numbers up to at least `B_{n-1}` are cached
    globally. Warning: this function is not currently threadsafe.


Bounding
-------------------------------------------------------------------------------

.. function:: long 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.


Computation of single Bernoulli numbers
-------------------------------------------------------------------------------

.. function:: void bernoulli_fmprb_ui_zeta(fmprb_t b, ulong n, long prec)

    Sets `b` to the numerical value of `B_n` accurate to *prec* bits,
    computed using the formula `B_{2n} = (-1)^{n+1} 2 (2n)! \zeta(2n) / (2 \pi)^n`.

    To avoid potential infinite recursion, we explicitly call the
    Euler product implementation of the zeta function.
    We therefore assume that the precision is small
    enough and `n` large enough for the Euler product to converge
    rapidly (otherwise this function will effectively hang).

.. function:: void bernoulli_fmprb_ui(fmprb_t b, ulong n, long prec)

    Sets `b` to the numerical value of `B_n` accurate to *prec* bits,
    computed by a division of the exact fraction if `B_n` is in
    the global cache or the exact numerator roughly is larger than
    *prec* bits, and using :func:`bernoulli_fmprb_ui_zeta`
    otherwise. This function reads `B_n` from the global cache
    if the number is already cached, but does not automatically extend
    the cache by itself.

.. function:: void _bernoulli_fmpq_ui_zeta(fmpz_t num, fmpz_t den, ulong n)

    Sets *num* and *den* to the reduced numerator and denominator
    of the Bernoulli number `B_n`.

    This function computes the denominator `d` using von Staudt-Clausen
    theorem, numerically approximates `B_n` using :func:`bernoulli_fmprb_ui_zeta`,
    and then rounds `d B_n` to the correct numerator.
    If the working precision is insufficient to determine the numerator,
    the function prints a warning message and retries with increased
    precision (this should not be expected to happen).

.. function:: void _bernoulli_fmpq_ui(fmpz_t num, fmpz_t den, ulong n)

.. function:: 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.