Bell numbers

The Bell number $B_n$ is the number of partitions of a set of $n$ elements: $B_0, B_1, \ldots = 1, 1, 2, 5, 15, 52, 203, \ldots$ In terms of Stirling numbers of the second kind, $B_n = \sum_{k=0}^n S(n,k)$, and the exponential generating function is $\sum_{n \ge 0} B_n x^n/n! = \exp(e^x - 1)$ [Comtet1974]. The numbers grow slightly slower than $n!$: $\log B_n = n(\log n - \log\log n - 1 + o(1))$, 6403 bits for $n = 1000$.

The module computes an isolated $B_n$ (arith_bell_number), the vector $B_0, \ldots, B_{n-1}$ (arith_bell_number_vec), and the same modulo a word-size integer (arith_bell_number_nmod, arith_bell_number_nmod_vec), which is much cheaper when only the residue is needed.

Isolated values#

arith_bell_number reads $B_n$ from a table for $n \le 25$ (64-bit), runs the Bell triangle in fixed-size arrays of up to three words for $n \le 58$, uses Dobiński's formula

for $n \le 3400$, in an exact form without floating-point arithmetic (arith_bell_number_dobinski), and the multimodular algorithm above (arith_bell_number_multi_mod). The latter computes $B_n$ modulo primes $p > n$ of $\mathtt{FLINT\_BITS} - 3$ bits with the identity

which needs $O(n)$ multiplications modulo $p$ per prime, and combines the residues by Chinese remaindering; the primes are processed in parallel by the threads of flint_set_num_threads. About $n \log_2 n / \mathtt{FLINT\_BITS}$ primes are needed, so the residues cost $O(n^2 \log n)$ bit operations.

Vectors#

The Bell triangle (Aitken's array) computes $B_0, \ldots, B_{n-1}$ with $n^2/2$ additions: each row starts with the last entry of the previous row, and each further entry is the sum of its left neighbour and the entry above it; the rows start with the Bell numbers. Over $\mathbb{Z}$ this costs $O(n^3 \log n)$ bit operations and is used for $n < 5000$; for longer vectors arith_bell_number_vec_multi_mod computes the vector modulo word-size primes and reconstructs each entry with only as many primes as it needs.

Modulo $m$, arith_bell_number_nmod_vec chooses between three methods: the triangle ($O(n^2)$ operations, any $m$, used for $n < 300$), the power series $\exp(e^x - 1)$ ($O(M(n))$ operations, requires $(n-1)!$ to be invertible modulo $m$, so in particular $m \ge n$ for prime $m$), and the ordinary generating function

summed by binary splitting ($O(M(n) \log n)$ operations, any $m$).

Example#

The program prints the first Bell numbers, the size of $B_{1000}$, and checks arith_bell_number_nmod against the residue of the exact value for the prime $m = 10^9 + 7$ (where the series method applies) and for the modulus $m = 2^{10} = 1024$ (where it does not).

#include <flint/flint.h>
#include <flint/fmpz.h>
#include <flint/fmpz_vec.h>
#include <flint/nmod.h>
#include <flint/arith.h>

int main(void)
{
    fmpz * b;
    fmpz_t B;
    nmod_t mod;
    slong i;

    b = _fmpz_vec_init(10);
    arith_bell_number_vec(b, 10);
    flint_printf("B_0..B_9:");
    for (i = 0; i < 10; i++)
    {
        flint_printf(" ");
        fmpz_print(b + i);
    }
    flint_printf("\n");

    fmpz_init(B);
    arith_bell_number(B, 1000);
    flint_printf("B_1000 has %wu bits\n", fmpz_bits(B));

    nmod_init(&mod, 1000000007);
    flint_printf("mod 10^9+7: %wu %wu\n",
        arith_bell_number_nmod(1000, mod), fmpz_fdiv_ui(B, mod.n));
    nmod_init(&mod, 1024);
    flint_printf("mod 1024:   %wu %wu\n",
        arith_bell_number_nmod(1000, mod), fmpz_fdiv_ui(B, mod.n));

    fmpz_clear(B);
    _fmpz_vec_clear(b, 10);
    return 0;
}

B_0..B_9: 1 1 2 5 15 52 203 877 4140 21147
B_1000 has 6403 bits
mod 10^9+7: 465231251 465231251
mod 1024:   443 443