[R] Stirling numbers

[Martin Maechler]

>>>>> "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<-------------------------------------------------------------------