Binomial Coefficients and Stirling Numbers in C#

As long as I’m on a roll with my posts on number theory in C#, I thought I’d briefly discuss how to generate binomial coefficients, and then move on to Stirling numbers of the second kind, both of which are extremely useful in combinatorics and finite calculus.

Some prerequisites for calculating binomial coefficients include a standard factorial function, nothing special:

1
2
3
4
5
6
7
long factorial(long n)
{
    if (n == 0) return 1;
    long t = n;
    while(n-- > 2) t *= n;
    return t;
}

Note that I use longs whenever possible because, despite the performance hit from 64-bit operations, it’s worth it to be able to work with numbers that are double the magnitude of ints.

We will also need a function for calculating a falling power of a number, which is defined as: $$x^{\underline{n}} = x(x-1)(x-2)\ldots(x-(n-1))$$You’ll see why in a moment.

1
2
3
4
5
6
long fallingPower(long n, long p)
{
    long t = 1;
    for (long i = 0; i < p; i++) t *= n--;
    return t;
}

Notice that the above function handles the case $$n^{\underline{0}}$$ returning 1. However, it does not handle the case of p > n, which would give an incorrect result, but this is not necessary for our purposes.

Binomial Coefficients

Recall that the definition for the binomial coefficient is $${n \choose k} = \frac{n!}{k!(n-k)!}$$ However, using this exact formula to compute binomial coefficients is a bit naive. If we use falling powers (sometimes called falling factorials), the above formula easily reduces to: $$\frac{n^{\underline{k}}}{k!}$$

We can improve the algorithm a bit more by adding the condition:

$${n \choose k} =
\begin{cases}
\frac{n^{\underline{k}}}{k!} \quad \mbox{if } k \leq \lfloor n/2 \rfloor,\\
\frac{n^{\underline{n-k}}}{(n-k)!} \quad \mbox{if } k > \lfloor n/2 \rfloor.
\end{cases}
$$

Not only is this algorithm faster, but it can also handle larger coefficients than the original formula, since neither the falling power nor the factorial ever gets larger than n/2.
The code for this is straightforward:

1
2
3
4
5
6
long binomialCoeff(long n, long k)
{
    if ((k < 0) || (k > n)) return 0;
    k = (k > n / 2) ? n - k : k;
    return fallingPower(n, k) / factorial(k);
}

However, this is still not as optimal as it can be. The most optimal approach would be to accumulate the falling power while dividing by each factor of the factorial in place. This would minimize the chance of overflow errors, and allow for even larger coefficients to be calculated. The disadvantage of this algorithm is the necessary use of floating-point math:

1
2
3
4
5
6
7
long binomialCoeff(long n, long k) {
    if ((k < 0) || (k > n)) return 0;
    k = (k > n / 2) ? n - k : k;
    double a = 1;
    for (long i = 1; i <= k; i++) a = (a * (n-k+i)) / i;
    return a + 0.5;
}

Stirling Numbers

Stirling numbers (of the second kind) are useful for, among other things, enumerating the coefficients of the falling-power expansion of a regular power. For example, how would we express x3 in terms of falling powers of x? That is, how do we arrive at an equation of the form $$x^3 = ax^{\underline 3} + bx^{\underline 2} + cx^{\underline 1}$$ Well, we could just solve the equation directly, but this would get unwieldy for higher powers. A neat way of doing this involves the use of Stirling numbers of the second kind, $$\left\{\begin{matrix} n \\ k \end{matrix}\right\}$$ A useful theorem for computing these numbers is

$$\left\{\begin{matrix} n \\ k \end{matrix}\right\}=\frac{1}{k!} \sum_{i=0}^{k}(-1)^i{k \choose i}(k-i)^n$$

1
2
3
4
5
6
7
8
9
10
11
long stirling(long n, long k)
{
    long sum = 0, neg = 1;
    for (long i = 0; i <= k; i++)
    {
        sum += neg * binomialCoeff(k, i) * pow(k - i, n);
        neg = -neg;
    }
    sum /= factorial(k);
    return sum;
}

With the Stirling numbers in hand, we can now obtain the coefficients for falling power expansions:

$$x^m = \sum_{k=0}^{m}\left\{\begin{matrix} m \\ k \end{matrix}\right\}x^{\underline k}$$

Share this article:
  • Facebook
  • LinkedIn
  • Reddit
  • Twitter

2 thoughts on “Binomial Coefficients and Stirling Numbers in C#

  1. Thialfi

    Hi there,
    there is no need of using float in the second binomial coefficient algorithm. If you make sure that multiplication comes first in the iteration (like you did), the result of the division is always an natural number (integer, long, …).
    In the first step you got one number divided by one. In the next step you got another number. Either the first or the second number is even and therefore the product of these two is even and can be divided by 2. In step 3, one of the three numbers is divisible by 3 and so is the product (even so we already have divided by 2). And so on. Because k is not exceeding n/2, we got no problems in division. The result in each step will always be natural. That’s why you don’t need float. Just stick to long.

    Reply
  2. viktor schausberger

    Using Binomial Coefficient, I would like to see a code able to generate all the possible conditional combination of N elemente taken K at the time, where the conditions are for example,
    how many odds/evens, sum range with exceptions, storage a prexisted list of combinations, comapare the list generate for the code with the storage list and eliminated duplicates. somebody can do that?

    Reply

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>