1Function: stirling
2Section: combinatorics
3C-Name: stirling
4Prototype: LLD1,L,
5Help: stirling(n,k,{flag=1}): if flag=1 (default) return the Stirling number
6 of the first kind s(n,k), if flag=2, return the Stirling number of the second
7 kind S(n,k).
8Doc: \idx{Stirling number} of the first kind $s(n,k)$ ($\fl=1$, default) or
9 of the second kind $S(n,k)$ (\fl=2), where $n$, $k$ are nonnegative
10 integers. The former is $(-1)^{n-k}$ times the
11 number of permutations of $n$ symbols with exactly $k$ cycles; the latter is
12 the number of ways of partitioning a set of $n$ elements into $k$ nonempty
13 subsets. Note that if all $s(n,k)$ are needed, it is much faster to compute
14 $$\sum_k s(n,k) x^k = x(x-1)\dots(x-n+1).$$
15 Similarly, if a large number of $S(n,k)$ are needed for the same $k$,
16 one should use
17 $$\sum_n S(n,k) x^n = \dfrac{x^k}{(1-x)\dots(1-kx)}.$$
18 (Should be implemented using a divide and conquer product.) Here are
19 simple variants for $n$ fixed:
20 \bprog
21 /* list of s(n,k), k = 1..n */
22 vecstirling(n) = Vec( factorback(vector(n-1,i,1-i*'x)) )
23
24 /* list of S(n,k), k = 1..n */
25 vecstirling2(n) =
26 { my(Q = x^(n-1), t);
27   vector(n, i, t = divrem(Q, x-i); Q=t[1]; simplify(t[2]));
28 }
29
30 /* Bell numbers, B_n = B[n+1] = sum(k = 0, n, S(n,k)), n = 0..N */
31 vecbell(N)=
32 { my (B = vector(N+1));
33   B[1] = B[2] = 1;
34   for (n = 2, N,
35     my (C = binomial(n-1));
36     B[n+1] = sum(k = 1, n, C[k]*B[k]);
37   ); B;
38 }
39 @eprog
40Variant: Also available are \fun{GEN}{stirling1}{ulong n, ulong k}
41 ($\fl=1$) and \fun{GEN}{stirling2}{ulong n, ulong k} ($\fl=2$).
42