Isolated Bernoulli numbers

bernoulli_fmpq_ui is the function to use for a single Bernoulli number B_n. For a range of consecutive values, bernoulli_fmpq_vec_no_cache or the cache (bernoulli_cache_compute) is much faster than repeated calls, and for the residue of B_n modulo a prime, bernoulli_mod_p_harvey avoids computing B_n at all.

How B_n is computed#

The denominator of B_n is known in advance: by the theorem of von Staudt and Clausen it is the product of the primes p with (p - 1) \mid n (arith_bernoulli_number_denom). It has few bits, about O(\log n \cdot n^{\epsilon}) in the worst case, so the work is in the numerator, which has about n \log_2 n bits: \log_2 |B_n| \approx n \log_2(n / (2 \pi e)). Two methods find it.

• The zeta function. For even n \ge 2,

  B_n = (-1)^{n/2+1} \frac{2\, n!}{(2\pi)^n} \zeta(n),

  and \zeta(n) = 1 + 2^{-n} + 3^{-n} + \cdots converges so fast that about n^{1 + o(1)} bits of precision need only O(n / \log n) terms of the series, or an Euler product. Multiplying the approximation by the denominator gives a ball containing one integer, the numerator. _bernoulli_fmpq_ui_zeta does this with arb_bernoulli_ui_zeta; the cost is O(n^2 \log^{O(1)} n) bit operations.

• Multimodular. Harvey's algorithm computes B_n \bmod p in time O(p \log^{O(1)} p), independent of n. With primes up to about n \log n this gives the numerator modulo a number larger than itself in total time O(n^2 \log^{O(1)} n) [Har2010], with a better constant and trivial parallelism. _bernoulli_fmpq_ui_multi_mod computes only the low bits this way, a fraction alpha of the total, and the high bits from the zeta function, which is the cheapest combination in practice.

bernoulli_fmpq_ui first looks in the cache, then uses the zeta method for n < 18000 and the combined method above, where most of the time is spent on the residues, which run in parallel with flint_set_num_threads.

Residues modulo primes#

bernoulli_mod_p_harvey returns B_n \bmod p for a prime p < 2^{53}. It reduces n modulo p - 1 with Kummer's congruence B_n/n \equiv B_m/m \pmod p for n \equiv m \not\equiv 0 \pmod{p-1}, then evaluates a sum of O(p) terms whose inner loop is a table lookup on the binary expansion of 1/p. When p - 1 divides n, p divides the denominator of B_n and the function returns UWORD_MAX. A result of 0 means that p divides the numerator, as for p = 691 and B_{12} = -691/2730; primes with this property for some n < p - 1 are the irregular primes.

Example#

#include <flint/flint.h>
#include <flint/fmpz.h>
#include <flint/fmpq.h>
#include <flint/bernoulli.h>

int main(void)
{
    fmpq_t b;
    fmpz_t p, r;

    fmpq_init(b);
    fmpz_init(p);
    fmpz_init(r);

    bernoulli_fmpq_ui(b, 20);
    flint_printf("B_20 = ");
    fmpq_print(b);
    flint_printf("\n");

    bernoulli_fmpq_ui(b, 10000);
    flint_printf("B_10000: numerator of %wu bits, denominator ",
        fmpz_bits(fmpq_numref(b)));
    fmpz_print(fmpq_denref(b));
    flint_printf("\n");

    /* compare with Harvey's algorithm modulo p = 1000003 */
    fmpz_set_ui(p, 1000003);
    fmpq_mod_fmpz(r, b, p);
    flint_printf("B_10000 mod p: %wu (Harvey), ",
        bernoulli_mod_p_harvey(10000, 1000003));
    fmpz_print(r);
    flint_printf(" (exact)\n");

    flint_printf("B_12 mod 691 = %wu\n", bernoulli_mod_p_harvey(12, 691));
    flint_printf("B_12 mod 13 is undefined: %d\n",
        bernoulli_mod_p_harvey(12, 13) == UWORD_MAX);

    fmpq_clear(b);
    fmpz_clear(p);
    fmpz_clear(r);
    return 0;
}

B_20 = -174611/330
B_10000: numerator of 91986 bits, denominator 2338224387510
B_10000 mod p: 932435 (Harvey), 932435 (exact)
B_12 mod 691 = 0
B_12 mod 13 is undefined: 1