# [R] Stirling numbers

Martin Maechler maechler at stat.math.ethz.ch
Wed Jul 19 16:19:15 CEST 2006

```>>>>> "Robin" == Robin Hankin <r.hankin at noc.soton.ac.uk>
>>>>>     on Wed, 19 Jul 2006 14:25:21 +0100 writes:

Robin> Hi anyone coded up Stirling numbers in R?

"Sure" ;-)

Robin> [I need unsigned Stirling numbers of the first kind]

but with my quick search, I can only see those for which I had
"2nd kind" :

-o<--o<-------------------------------------------------------------------

##-- Stirling numbers of the 2nd kind
##-- (Abramowitz/Stegun: 24,1,4 (p. 824-5 ; Table 24.4, p.835)

##> S^{(m)}_n = number of ways of partitioning a set of \$n\$ elements into \$m\$
##>	non-empty subsets

Stirling2 <- function(n,m)
{
## Purpose:  Stirling Numbers of the 2-nd kind
## 		S^{(m)}_n = number of ways of partitioning a set of
##                      \$n\$ elements into \$m\$ non-empty subsets
## Author: Martin Maechler, Date:  May 28 1992, 23:42
## ----------------------------------------------------------------
## Abramowitz/Stegun: 24,1,4 (p. 824-5 ; Table 24.4, p.835)
## Closed Form : p.824 "C."
## ----------------------------------------------------------------

if (0 > m || m > n) stop("'m' must be in 0..n !")
k <- 0:m
sig <- rep(c(1,-1)*(-1)^m, length= m+1)# 1 for m=0; -1 1 (m=1)
## The following gives rounding errors for (25,5) :
## r <- sum( sig * k^n /(gamma(k+1)*gamma(m+1-k)) )
ga <- gamma(k+1)
round(sum( sig * k^n /(ga * rev(ga))))
}

options(digits=15)
for (n in 11:31) cat("n =", n," S_n(5) =", Stirling2(n,5), "\n")
n <- 25
for(k in c(3,5,10))
cat(" S_",n,"^(",formatC(k,wid=2),") = ", Stirling2(n,k),"\n",sep = "")

for (n in 1:15)
cat(formatC(n,w=2),":", sapply(1:n, Stirling2, n = n),"\n")
##  1 : 1
##  2 : 1 1
##  3 : 1 3 1
##  4 : 1 7 6 1
##  5 : 1 15 25 10 1
##  6 : 1 31 90 65 15 1
##  7 : 1 63 301 350 140 21 1
##  8 : 1 127 966 1701 1050 266 28 1
##  9 : 1 255 3025 7770 6951 2646 462 36 1
## 10 : 1 511 9330 34105 42525 22827 5880 750 45 1
## 11 : 1 1023 28501 145750 246730 179487 63987 11880 1155 55 1
## 12 : 1 2047 86526 611501 1379400 1323652 627396 159027 22275 1705 66 1
## 13 : 1 4095 261625 2532530 7508501 9321312 5715424 1899612 359502 39325 2431 78 1
## 14 : 1 8191 788970 10391745 40075035 63436373 49329280 20912320 5135130 752752 66066 3367 91 1
## 15 : 1 16383 2375101 42355950 210766920 420693273 408741333 216627840 67128490 12662650 1479478 106470 4550 105 1

plot(6:30, sapply(6:30, Stirling2, m=5), log = "y", type = "l")
## Abramowitz says something like  S2(n,m) ~ m^n / m!
##                                 ------------------
nn <- 6:30; sapply(nn, Stirling2, m=5)  / (5^nn / gamma(5+1))
## 0.114 0.21 ...... 0.99033 0.992266 0.993812

-o<--o<-------------------------------------------------------------------

```