Stirling and Bell numbers

The Stirling numbers s(n,k) and S(n,k) and the Bell numbers B_n are defined on Stirling numbers and Bell numbers. All of these numbers have \Theta(n \log n) bits for typical k, so the best possible cost of an isolated value is quasi-linear in n \log n; the methods below are simpler and, apart from the binary splitting method for s(n,k), cost about n^2 word operations, which is adequate up to n of the order of 10^5. The references for the combinatorics are [Comtet1974] and [GKP1994, §6.1].

Stirling numbers of the first kind#

The unsigned numbers are the coefficients of the rising factorial:

A row is computed by binary splitting: the product of the n - 1 linear factors x + j, 1 \le j < n, split recursively in halves, with fast polynomial multiplication at the top levels. The coefficients of the product of the factors with j \in [a, b) have at most (b-a) \log_2 \max(b, b - a) bits, so the base cases with small products are done in single words, and the total cost is O(M(n^2 \log n) \log n) bit operations, quasi-linear in the size of the row.

For one coefficient |s(n,k)|, only a truncated product is needed: for k < n/2 the coefficients up to x^{k-1} of \prod_{j=1}^{n-1} (x + j), and otherwise the coefficients up to x^{n-k} of the reversed product \prod_{j=1}^{n-1} (1 + jx). The last multiplication of the binary splitting tree computes only the wanted coefficient, as a dot product. For k very close to n the exponential generating function

is cheaper, since the power series has length only n - k + 1; FLINT uses it for n > 140 and k > 0.87 n. The signed numbers are s(n,k) = (-1)^{n-k} |s(n,k)|.

Stirling numbers of the second kind#

By inclusion–exclusion over the empty blocks,

a sum of k + 1 terms of about n \log k bits. Two observations reduce the work. The powers j^n are needed only for odd j: for j = 2^e j' with j' odd, j^n = j'^n \cdot 2^{en}, a shift. And the binomial coefficients are symmetric, so only half of them are computed. Modulo a word-size prime the powers of composite j are products of the powers of two factors, read from a divisor table, so that only \pi(k) exponentiations are needed; the sum is then a dot product.

Isolated values use:

• the triangular recurrence S(n,k) = k S(n-1, k) + S(n-1, k-1) in one or two words while S(n,k) fits, that is for n \le 43 on 64-bit machines;
• the exact power sum for k \le c_1 n, where the terms are small;

• the exponential generating function

  S(n,k) = \frac{n!}{k!} [x^{n-k}] \left(\frac{e^x - 1}{x}\right)^k

  for k \ge c_2 n, a power of a series of length n - k + 1;

• the power sum modulo \mathtt{FLINT\_BITS} - 1 bit primes followed by Chinese remaindering in between.

The cutoffs are c_1 = c_2 = 0.9 for n < 200; for larger n, c_1 is 0.95 e^{-0.00022 n} (for n < 3000) or 1500/n, kept in [0.0002, 0.8], and c_2 = \min(0.92 + 0.005 \log n, 0.98). The number of primes is determined by the bound

of Rennie and Dobson [RenDob1969].

A row S(n, 0), \ldots, S(n, \mathit{klen} - 1) is a convolution of the sequences j^n/j! and (-1)^i/i!:

FLINT uses the triangular recurrence for n \le 80, the convolution over \mathbb{Z} with the common denominator \ell!^2, \ell = \min(\mathit{klen}, n) - 1, for \mathit{klen} < n/2, and otherwise the convolution modulo word-size primes with Chinese remaindering, where each entry uses only as many primes as its own bound requires.

Bell numbers#

Summing the power sum formula over k and exchanging the sums gives

and letting the inner sum run to infinity gives Dobiński's formula [Dob1877]

Dobiński's formula, exactly. With P = N! \sum_{k=0}^N k^n/k! and Q = N! \sum_{k=0}^N 1/k!, both integers,

For N \ge n each k > N has k^n > n^n > B_n, so P/Q \le B_n; and N \approx n(1 + 1.2/\log n) makes the tail T_P small enough that P/Q > B_n - 1. Hence B_n = \lceil P/Q \rceil, computed with integer arithmetic only (arith_bell_number_dobinski). The terms are accumulated by Horner's rule in k, so the cost is O(N) multiplications of a large number by a word plus the N/2 powers k^n for odd k, about O(n^2 \log^2 n) bit operations.

Multimodular sum. Modulo a prime p > n, the double sum above is evaluated with O(n) multiplications: the factorials are precomputed as n!/i!, so that the inner sums are prefix sums and only one inversion (of n!^2) is needed at the end, and the powers j^n of composite j come from a divisor table. FLINT accumulates the products in three words and reduces once, uses Montgomery arithmetic (Montgomery multiplication (REDC)) with primes of \mathtt{FLINT\_BITS} - 3 bits, and processes the primes in parallel (arith_bell_number_multi_mod). The number of primes comes from de Bruijn's asymptotic expansion [deB1958]

Vectors. The Bell triangle (Aitken's array) a_{0,0} = 1, a_{i,0} = a_{i-1,i-1}, a_{i,j} = a_{i,j-1} + a_{i-1,j-1} has a_{i,0} = B_i and needs n^2/2 additions for B_0, \ldots, B_{n-1}. Modulo m there are two power series methods: the exponential generating function \exp(e^x - 1), computed with one power series exponential in O(M(n)) operations if (n-1)! is invertible modulo m, and the ordinary generating function

which follows from \sum_n S(n,k) x^n = x^k / \prod_{j=1}^k (1 - jx). Truncated at k = n - 1, the sum is a fraction R/Q that binary splitting computes as polynomials of degree n, followed by one power series division: O(M(n) \log n) operations, for any modulus.

Use in FLINT#

The generic ring functions gr_stirling_s1u_uiui, gr_stirling_s2_uiui, gr_stirling_s1_ui_vec (and the other kinds), gr_bellnum_ui and gr_bellnum_vec call these functions over \mathbb{Z} (and arith_bell_number_nmod and arith_bell_number_nmod_vec over \mathbb{Z}/n\mathbb{Z}) and convert the result for other rings.

Bell numbers modulo 2 and 3 are periodic with periods 3 and 13, by Touchard's congruence B_{n+p} \equiv B_n + B_{n+1} \pmod p for prime p [Comtet1974].