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