Stirling numbers

The Stirling numbers of the first kind connect the falling and rising factorials with the powers of x, and those of the second kind the powers with the falling factorials:

The unsigned numbers |s(n,k)| = (-1)^{n-k} s(n,k) are the coefficients of the rising factorial x(x+1)\cdots(x+n-1) and count the permutations of n elements with k cycles; S(n,k) counts the partitions of an n-set into k blocks [Comtet1974]. All of them are zero for k > n, and s(n,n) = S(n,n) = 1, s(n,0) = S(n,0) = 0 for n \ge 1.

The module offers three interfaces, which use different algorithms:

• arith_stirling_number_1u, arith_stirling_number_1 and arith_stirling_number_2 compute one value;
• the _vec functions compute a row k = 0, \ldots, \mathit{klen} - 1 for fixed n;
• the matrix functions fill a whole table of rows 0 \le n < r.

To compute many values, use the row or matrix functions: they cost about as much as a few isolated values.

Isolated values#

The values s(n,k) with k near 0 or near n are cheap: |s(n,1)| = (n-1)!, S(n,2) = 2^{n-1} - 1, S(n, n-1) = |s(n, n-1)| = \binom{n}{2}. In general, the functions choose by n and the ratio k/n:

• |s(n,k)| is a coefficient of the rising factorial, which arith_stirling_number_1u expands by binary splitting truncated to the needed length: as (1+x)(2+x)\cdots(n-1+x) to length k for k < n/2, and as (1+x)(1+2x)\cdots(1+(n-1)x) to length n-k+1 otherwise, with only the wanted coefficient extracted from the last product. For n > 140 and k > 0.87 n it uses instead the exponential generating function |s(n,k)| = (n!/k!)\, [x^{n-k}] \bigl(-\log(1-x)/x\bigr)^k.

• S(n,k) is computed by arith_stirling_number_2 from the triangular recurrence S(n,k) = k S(n-1,k) + S(n-1,k-1) in one or two words for n \le 43 (on 64-bit machines), and otherwise from one of

  S(n,k) = \frac{1}{k!} \sum_{j=0}^{k} (-1)^{k-j} \binom{k}{j} j^n \quad\text{or}\quad S(n,k) = \frac{n!}{k!} [x^{n-k}] \left(\frac{e^x - 1}{x}\right)^k :

  the power sum exactly for small k/n (below about 0.95 e^{-0.00022 n}, and at most 0.8), the power series for k/n above about 0.92 + 0.005 \log n, and the power sum modulo word-size primes with Chinese remaindering in between.

The signed version arith_stirling_number_1 is the unsigned one with the sign (-1)^{n-k}.

Rows#

A full row of |s(n,k)| is the coefficient list of the rising factorial, which arith_stirling_number_1u_vec computes by binary splitting with fast polynomial multiplication, truncated to klen coefficients. The row of S(n,k) comes from the triangular recurrence for n \le 80; for larger n, it is the convolution

a polynomial product computed exactly over \mathbb{Z} (with a common denominator) for \mathit{klen} < n/2 and modulo word-size primes, followed by Chinese remaindering, otherwise. klen may exceed n + 1; the extra entries are zero.

Matrices#

arith_stirling_matrix_1u, arith_stirling_matrix_1 and arith_stirling_matrix_2 are wrappers of gr_mat_stirling, which fills row n from row n - 1 with the recurrences

The matrix may have any shape; its entries above the diagonal are zero. By the defining identities, the square matrices of s and S are inverses of each other.

Example#

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

int main(void)
{
    fmpz_t s;
    fmpz * row;
    fmpz_mat_t A, B, C;
    slong k;

    fmpz_init(s);
    arith_stirling_number_1(s, 10, 3);
    flint_printf("s(10,3) = ");
    fmpz_print(s);
    arith_stirling_number_2(s, 10, 3);
    flint_printf(", S(10,3) = ");
    fmpz_print(s);
    flint_printf("\n");

    row = _fmpz_vec_init(7);
    arith_stirling_number_2_vec(row, 6, 7);
    flint_printf("S(6,k):");
    for (k = 0; k < 7; k++)
    {
        flint_printf(" ");
        fmpz_print(row + k);
    }
    flint_printf("\n");

    arith_stirling_number_2(s, 2000, 1000);
    flint_printf("S(2000,1000) has %wu bits\n", fmpz_bits(s));

    fmpz_mat_init(A, 8, 8);
    fmpz_mat_init(B, 8, 8);
    fmpz_mat_init(C, 8, 8);
    arith_stirling_matrix_1(A);
    arith_stirling_matrix_2(B);
    fmpz_mat_mul(C, A, B);
    flint_printf("s S = I: %d\n", fmpz_mat_is_one(C));

    fmpz_mat_clear(A);
    fmpz_mat_clear(B);
    fmpz_mat_clear(C);
    _fmpz_vec_clear(row, 7);
    fmpz_clear(s);
    return 0;
}

s(10,3) = -1172700, S(10,3) = 9330
S(6,k): 0 1 31 90 65 15 1
S(2000,1000) has 11144 bits
s S = I: 1