Fastest way to generate binomial coefficients

Question

I need to calculate combinations for a number.

What is the fastest way to calculate nCp where n>>p?

I need a fast way to generate binomial coefficients for an polynomial equation and I need to get the coefficient of all the terms and store it in an array.

(a+b)^n = a^n + nC1 a^(n-1) * b + nC2 a^(n-2) * ............ +nC(n-1) a * b^(n-1) + b^n

What is the most efficient way to calculate nCp ??

Solution 1

If you want complete expansions for large values of n, FFT convolution might be the fastest way. In the case of a binomial expansion with equal coefficients (e.g. a series of fair coin tosses) and an even order (e.g. number of tosses) you can exploit symmetries thus:

Theory

Represent the results of two coin tosses (e.g. half the difference between the total number of heads and tails) with the expression A + A*cos(Pi*n/N). N is the number of samples in your buffer - a binomial expansion of even order O will have O+1 coefficients and require a buffer of N >= O/2 + 1 samples - n is the sample number being generated, and A is a scale factor that will usually be either 2 (for generating binomial coefficients) or 0.5 (for generating a binomial probability distribution).

Notice that, in frequency, this expression resembles the binomial distribution of those two coin tosses - there are three symmetrical spikes at positions corresponding to the number (heads-tails)/2. Since modelling the overall probability distribution of independent events requires convolving their distributions, we want to convolve our expression in the frequency domain, which is equivalent to multiplication in the time domain.

In other words, by raising our cosine expression for the result of two tosses to a power (e.g. to simulate 500 tosses, raise it to the power of 250 since it already represents a pair), we can arrange for the binomial distribution for a large number to appear in the frequency domain. Since this is all real and even, we can substitute the DCT-I for the DFT to improve efficiency.

Algorithm

1. decide on a buffer size, N, that is at least O/2 + 1 and can be conveniently DCTed
2. initialise it with the expression pow(A + A*cos(Pi*n/N),O/2)
3. apply the forward DCT-I
4. read out the coefficients from the buffer - the first number is the central peak where heads=tails, and subsequent entries correspond to symmetrical pairs successively further from the centre

Accuracy

There's a limit to how high O can be before accumulated floating-point rounding errors rob you of accurate integer values for the coefficients, but I'd guess the number is pretty high. Double-precision floating-point can represent 53-bit integers with complete accuracy, and I'm going to ignore the rounding loss involved in the use of pow() because the generating expression will take place in FP registers, giving us an extra 11 bits of mantissa to absorb the rounding error on Intel platforms. So assuming we use a 1024-point DCT-I implemented via the FFT, that means losing 10 bits' accuracy to rounding error during the transform and not much else, leaving us with ~43 bits of clean representation. I don't know what order of binomial expansion generates coefficients of that size, but I dare say it's big enough for your needs.

Asymmetrical expansions

If you want the asymmetrical expansions for unequal coefficients of a and b, you'll need to use a two-sided (complex) DFT and a complex pow() function. Generate the expression A*A*e^(-Pi*i*n/N) + A*B + B*B*e^(+Pi*i*n/N) [using the complex pow() function to raise it to the power of half the expansion order] and DFT it. What you have in the buffer is, again, the central point (but not the maximum if A and B are very different) at offset zero, and it is followed by the upper half of the distribution. The upper half of the buffer will contain the lower half of the distribution, corresponding to heads-minus-tails values that are negative.

Notice that the source data is Hermitian symmetrical (the second half of the input buffer is the complex conjugate of the first), so this algorithm is not optimal and can be performed using a complex-to-complex FFT of half the required size for optimum efficiency.

Needless to say, all the complex exponentiation will chew more CPU time and hurt accuracy compared to the purely real algorithm for symmetrical distributions above.

Solution 2

You cau use dynamic programming in order to generate binomial coefficients

You can create an array and than use O(N^2) loop to fill it

C[n, k] = C[n-1, k-1] + C[n-1, k];

where

C[1, 1] = C[n, n] = 1

After that in your program you can get the C(n, k) value just looking at your 2D array at [n, k] indices

UPDATE smth like that

for (int k = 1; k <= K; k++) C[0][k] = 0;
for (int n = 0; n <= N; n++) C[n][0] = 1;

for (int n = 1; n <= N; n++)
   for (int k = 1; k <= K; k++)
      C[n][k] = C[n-1][k-1] + C[n-1][k];

where the N, K - maximum values of your n, k

Solution 3

If you need to compute them for all n, Ribtoks's answer is probably the best. For a single n, you're better off doing like this:

C[0] = 1
for (int k = 0; k < n; ++ k)
    C[k+1] = (C[k] * (n-k)) / (k+1)

The division is exact, if done after the multiplication.

And beware of overflowing with C[k] * (n-k) : use large enough integers.

Solution 4

I recently wrote a piece of code that needed to call for a binary coefficient about 10 million times. So I did a combination lookup-table/calculation approach that's still not too wasteful of memory. You might find it useful (and my code is in the public domain). The code is at

http://www.etceterology.com/fast-binomial-coefficients

It's been suggested that I inline the code here. A big honking lookup table seems like a waste, so here's the final function, and a Python script that generates the table:

extern long long bctable[]; /* See below */

long long binomial(int n, int k) {
    int i;
    long long b;
    assert(n >= 0 && k >= 0);

    if (0 == k || n == k) return 1LL;
    if (k > n) return 0LL;

    if (k > (n - k)) k = n - k;
    if (1 == k) return (long long)n;

    if (n <= 54 && k <= 54) {
        return bctable[(((n - 3) * (n - 3)) >> 2) + (k - 2)];
    }
    /* Last resort: actually calculate */
    b = 1LL;
    for (i = 1; i <= k; ++i) {
        b *= (n - (k - i));
        if (b < 0) return -1LL; /* Overflow */
        b /= i;
    }
    return b;
}

---

#!/usr/bin/env python3

import sys

class App(object):
    def __init__(self, max):
        self.table = [[0 for k in range(max + 1)] for n in range(max + 1)]
        self.max = max

    def build(self):
        for n in range(self.max + 1):
            for k in range(self.max + 1):
                if k == 0:      b = 1
                elif  k > n:    b = 0
                elif k == n:    b = 1
                elif k == 1:    b = n
                elif k > n-k:   b = self.table[n][n-k]
                else:
                    b = self.table[n-1][k] + self.table[n-1][k-1]
                self.table[n][k] = b

    def output(self, val):
        if val > 2**63: val = -1
        text = " {0}LL,".format(val)

        if self.column + len(text) > 76:
            print("\n   ", end = "")
            self.column = 3
        print(text, end = "")
        self.column += len(text)

    def dump(self):
        count = 0
        print("long long bctable[] = {", end="");

        self.column = 999
        for n in range(self.max + 1):
            for k in range(self.max + 1):
                if n < 4 or k < 2 or k > n-k:
                    continue
                self.output(self.table[n][k])
                count += 1
        print("\n}}; /* {0} Entries */".format(count));

    def run(self):
        self.build()
        self.dump()
        return 0

def main(args):
    return App(54).run()

if __name__ == "__main__":
    sys.exit(main(sys.argv))

Solution 5

This is my version:

def binomial(n, k):
if k == 0:
    return 1
elif 2*k > n:
    return binomial(n,n-k)
else:
    e = n-k+1
    for i in range(2,k+1):
        e *= (n-k+i)
        e /= i
    return e

View more solutions