Calculating Stirling Numbers using C#

A “Stirling number of the second kind” is the number of ways to group/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 and the seven groupings 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

How did I compute this? I wrote a C# function. First, I got one of many math definitions for S(n,k) from Wikipedia:

StirlingEquationSmall

To code a function that computes S(n,k) because the answers can be so big, I needed to use the C# BigInteger data type which is not visible by default to a C# program. The type lives in namespace System.Numerics. If you examine the equation above you can see you need a 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).

So, here’s Factorial:

static BigInteger Factorial(int n)
{
  if (n == 0) return 1;
  BigInteger ans = 1;
  for (int i = 1; i <= n; ++i)
    ans *= i;
  return ans;
}

Here’s a Power:

static BigInteger Power(int m, int p)
{
  // m raised to the pth power
  BigInteger result = 1;
  for (int i = 0; i < p; ++i)
    result = result * m;
  return result;
}

And here’s a Choose:

static BigInteger Choose(int n, int k)
{
  if (n < 0 || k < 0)
    throw new Exception("Negative argument in Choose");
  if (n < k) return 0; // special
  if (n == k) return 1; // short-circuit

  int delta, iMax;

  if (k < n - k) // ex: Choose(100,3)
  {
    delta = n - k;
    iMax = k;
  }
  else           // ex: Choose(100,97)
  {
    delta = k;
    iMax = n - k;
  }

  BigInteger ans = delta + 1;
  for (int i = 2; i <= iMax; ++i)
    ans = (ans * (delta + i)) / i;

  return ans;
}

The Choose function is not at all obvious, but that’s another topic. Putting it all together, here’s a C# method that computes S(n,k):

static BigInteger Stirling(int n, int k)
{
  BigInteger sum = 0;

  for (int j = 0; j <= k; ++j)
  {
    BigInteger a = Power(-1, k - j);
    BigInteger b = Choose(k, j);
    BigInteger c = Power(j, n);
    sum += a * b * c;
  }

  return sum / Factorial(k);
}

I called the function as

int n = 100;
in k = 3;
Console.WriteLine(Stirling(n, k));

All in all, an interesting function.

StirlingDemo

About these ads
This entry was posted in Machine Learning. Bookmark the permalink.

2 Responses to Calculating Stirling Numbers using C#

  1. mvaneerde says:

    Another common way to calculate Stirling numbers of the second kind is to build a table using the recurrence relation S(n, k) = k S(n – 1, k) + S(n – 1, k – 1).

    The motivation for the recurrence relation can be understood as follows:

    Suppose you are one of n people and you collectively need to form k groups.

    You could be antisocial and form a group by yourself. Then the other (n – 1) people can form (k – 1) groups in S(n – 1, k – 1) ways.

    Or you could decide to be social and be part of a group with other people. Let the other (n – 1) people form k groups; they can do this in S(n – 1, k) ways. You then have k choices for which group to join, so there are k S(n – 1, k) ways to do this.

    Total: S(n, k) = k S(n – 1, k) + S(n – 1, k – 1).

    • (author reply to Matthew) Very nicely explained. It’s be interesting to code up a recursive Stirling number method based on the recurrence relationship and then investigate its performance versus the iterative definition. JM

Comments are closed.