Well, given Martin's 2nd kind function and the fact that the 1st and 2nd
kind matrices are inverses, you can get at least some of what you want
by:
... # fill a matrix S2 with second kind numbers and
zeros> S2
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 1 0 0 0 0 0 0
[2,] 1 1 0 0 0 0 0
[3,] 1 3 1 0 0 0 0
[4,] 1 7 6 1 0 0 0
[5,] 1 15 25 10 1 0 0
[6,] 1 31 90 65 15 1 0
[7,] 1 63 301 350 140 21 1> abs(round(solve(S2)))
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 1 0 0 0 0 0 0
[2,] 1 1 0 0 0 0 0
[3,] 2 3 1 0 0 0 0
[4,] 6 11 6 1 0 0 0
[5,] 24 50 35 10 1 0 0
[6,] 120 274 225 85 15 1 0
[7,] 720 1764 1624 735 175 21 1
Not sure how big arguments you need.
HTH,
David L. Reiner
Rho Trading Securities, LLC
Chicago IL 60605
312-362-4963
-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Martin Maechler
Sent: Wednesday, July 19, 2006 9:19 AM
To: Robin Hankin
Cc: r-help at stat.math.ethz.ch
Subject: Re: [R] Stirling numbers
>>>>> "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<-----------------------------------------------------------------
--
______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.