[Edu-sig] Plain vanilla (Stirling numbers)

[Kirby Urner]

Stirling numbers add row-wise sort of like Pascal's Triangle,
except there's a multiplier at work on the first of each
pair of terms to be added.  Generators work well to give
successive rows, while wrappers control iteration, such
that row, column (n,k) value may be returned.

Stirling numbers of the 2nd type give the number of ways you
can divide n elements into k non-empty subsets, while those of
the 1st type give the number of ways you can arrange n elements
into k unique cycles or "necklaces" (order matters).

Examples:

  >>> from concrete import *
  >>> subsets(4,2)  # 4 things maybe be divided into 2 sets 7 ways
  7
  >>> cycles(4,2)   # but more cycles, as subsets of 3 have 2 cycles
  11

The implementation below also makes use of the enumerate() iterator,
new in 2.3, list comprehensions, and the built-in zip function.

====================================

  # Stirling Numbers

  # See: Graham, Knuth, Patashnik, Concrete Mathematics
  #     (Addison-Wesley, 2nd Ed., 1994), pg.258-261

  # Kirby Urner, Oregon Curriculum Network, March 1, 2003

  def stirling2():
     """
     Generate next row of Stirling numbers of type 2
     """
     row = [1]
     yield row
     while True:	
        mrow = [ i*j for (i,j) in enumerate(row) ]
        row  = [ i+j for (i,j) in zip(mrow + [0], [0] + row) ]
        yield row

  def stirling1():
     """
     Generate next row of Stirling numbers of type 1
     """
     row = [1]
     yield row
     while True:
        n = len(row)
        mrow = [ i * (n-1) for i in row ]
        row  = [ i+j for (i,j) in zip(mrow + [0], [0] + row) ]
        yield row

  def subsets(n,k):
     """
     Number of ways to partition a set of n things
     into k nonempty subsets
     """
     if n==0: return 1
     if k==1: return 1
     s = stirling2()
     for i in range(n+1):
	r = s.next()
     return r[k]

  def cycles(n,k):
     """
     Number of ways to arrange n objects into k cycles
     """
     if n==0: return 1
     if k==1: return 1
     s = stirling1()
     for i in range(n+1):
        r = s.next()
     return r[k]