Computing a Stirling Number of the Second Kind from Scratch Using Python

A “Stirling number of the second kind” — S(n, k) — is the number of ways to partition n items into k subsets. For example, the number of ways to group n=4 items into k=2 subsets is 7. Suppose the 4 items are a, b, c, d. Then S(4,2) = 7. The seven partitions are: (a | bcd), (b | acd), (c | abd), (d | abc), (ab |cd), (ac | bd), and (ad | bc).

The value of S(n, k) gets astronomically large, very quickly. For example, S(100, 3) =

85,896,253,455,335,221,205,584,888,180,155,511,368,666,317,646

I’m sure there is certainly a Stirling function in one of the Python libraries (probably SciPy or NumPy) but I decided to implement a Stirling function from scratch using raw Python, without any libraries.

I got one of many math definitions for S(n,k) from Wikipedia:

Next, I had to worry about dealing with huge integer values. But one of the extremely nice things about the Python language is that is has built-in support for arbitrarily large integers. Dealing with huge integers in many languages is a huge pain point.



The math definition above uses the Factorial function (outside the summation), an integer Power function (used in two places), and a Choose function. The Choose(n,k) function returns the number of ways to choose k items from n items, where order does not matter. For example, Choose(4,2) = 6 because there are six ways to select 2 items from 4 items. Suppose the four items are q, r, s, t. Then the six ways are (q r), (q s), (q t), (r s), (r t), and (s t).

I implemented functions my_factorial(n), my_power(n, p), and my_choose(n, k) and then used them to implement a my_stirling_second(n, k) function.

The moral of the story is that you don’t have to be a software deity of some sort to implement functions. Except for a few rare situations, most functions can be implemented by mere mortals.


Three illustrations of ancient Egyptian deities by an artist named “Yliade”. From left to right: Seshat, Nephthys, Serqet. I like the modern style of art applied to an ancient theme.


# stirling_second_demo.py
# note: // is integer division

def my_stirling_second(n, k):
  sum = 0
  for j in range(0, k+1):
    a = my_power(-1, k-j)
    b = my_choose(k, j)
    c = my_power(j, n)
    sum += a * b * c

  return sum // my_factorial(k)

def my_factorial(n):
  ans = 1
  for i in range(1,n+1):
    ans *= i
  return ans

def my_power(n, p):
  # n^p
  ans = 1
  for i in range(p):
    ans *= n
  return ans

def my_choose(n, k):
  if n == k: return 1
  if k "less-than" n-k:
    delta = n - k
    i_max = k
  else:
    delta = k
    i_max = n - k
  
  ans = delta + 1
  for i in range(2, i_max+1):
    ans = (ans * (delta + i)) // i
  return ans


def main():
  print("\nBegin Stirling numbers second kind demo ")

  nf = my_factorial(10)
  print("\nmy_factorial(10) = ")
  print("{:,}".format(nf))

  pwr = my_power(5, 3)
  print("\nmy_power(5,3) = ")
  print("{:,}".format(pwr))

  c = my_choose(10, 3)
  print("\nmy_choose(10,3) = ")
  print("{:,}".format(c))

  sn = my_stirling_second(100, 3)
  print("\nmy_stirling_second(100,3) = ")
  print("{:,}".format(sn))

  print("\nEnd demo \n")

if __name__ == "__main__":
  main()
This entry was posted in Miscellaneous. Bookmark the permalink.

Leave a Reply

Please log in using one of these methods to post your comment:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s