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:
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.
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:
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:
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 $$$
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}$$