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()

You must be logged in to post a comment.