Binomial Coefficients rounding error

Question

I have to calculate in c binomial coefficients of the expression (x+y)**n, with n very large (order of 500-1000). The first algo to calculate binomial coefficients that came to my mind was multiplicative formula. So I coded it into my program as

long double binomial(int k, int m)
{
    int i,j;
    long double num=1, den=1;
    j=m<(k-m)?m:(k-m);
    for(i=1;i<=j;i++)
    {
        num*=(k+1-i);
        den*=i;
    }
    return num/den; 
}

This code is really fast on a single core thread, compared for example to recursive formula, although the latter one is less subject to rounding errors since involves only sums and not divisions. So I wanted to test these algos for great values and tried to evaluate 500 choose 250 (order 10^160). I have found that the "relative error" is less than 10^(-19), so basically they are the same number, although they differ something like 10^141.

So I'm wondering: Is there a way to evaluate the order of the error of the calculation? And is there some fast way to calculate binomial coefficients which is more precise than the multiplicative formula? Since I don't know the precision of my algo I don't know where to truncate the stirling's series to get better results..

I've googled for some tables of binomial coefficients so I could copy from those, but the best one I've found stops at n=100...

Answer 1

To get exact integer results for small k and m, a better solution might be (a slight variation of your code) :

  unsigned long binomial(int k, int m)
  {
   int i,j; unsigned long num=1;
   j=m<(k-m)?m:(k-m);
   for(i=1;i<=j;i++)
   {
    num*=(k+1-i);
    num/=i;
   }
   return num;
  }

Every time you get a combinatorial number after doing the division num/=i, so you won't get truncated. To get approximate results for bigger k and m, your solution might be good. But beware that long double multiplication is already much slower than the multiplication and division of integers (unsigned long or size_t). If you want to get bigger numbers exact, probably a big integer class must be coded or included from a library. You can also google if there's fast factorial algorithm for n! of extremely big integer n. That may help with combinatorics, too. Stirling's formula is a good approximation for ln(n!) when n is large. It all depends on how accurate you want to be.

Answer 2

If you really want to use the multiplicative formula, I would recommend an exception based approach.

1. Implement the formula with large integers (long long for example)
2. Attempt division operations as soon as possible (as suggested by Zhuoran)
3. Add code to check correctness of every division and multiplication
4. Resolve incorrect divisions or multiplications, e.g.
   • try the division in loop proposed by Zhuoran, but if it fails resort back to the initial algorithm (accumulating the product of divisor in den)
   • store the unresolved multiplier, divisors in additional long integers and try to resolve them in next iteration loops
5. If you really use large numbers then your result might not fit in long integer. then in that case you can switch to long double or use your personal LongInteger storage.

This is a skeleton code, to give you an idea:

long long binomial_l(int k, int m)
{
    int i,j;
    long long num=1, den=1;
    j=m<(k-m)?m:(k-m);
    for(i=1;i<=j;i++)
    {
        int multiplier=(k+1-i);
        int divisor=i;
        long long candidate_num=num*multiplier;
        //check multiplication
        if((candidate_num/multiplier)!=num)
        {
            //resolve exception...
        }
        else
        {
            num=candidate_num;
        }

        candidate_num=num/divisor;
        //check division
        if((candidate_num*divisor)==num)
        {
            num=candidate_num;
        }
        else
        {
            //resolve exception
            den*=divisor;
            //this multiplication should also be checked...
        }
    }
    long long candidate_result= num/den; 
    if((candidate_result*den)==num)
    {
        return candidate_result;
    }
    // you should not get here if all exceptions are resolved
    return 0;
}

Answer 3

If you're just computing individual binomial coefficients C(n,k) with n fairly large but no larger than about 1750, then your best bet with a decent C library is to use the tgammal standard library function:

tgammal(n+1) / (tgammal(n-k+1) * tgammal(k+1))

Tested with the Gnu implementation of libm, that consistently produced results within a few ULP of the precise value, and generally better than solutions based on multiplying and dividing.

If k is small (or large) enough that the binomial coefficient does not overflow 64 bits of precision, then you can get a precise result by alternately multiplying and dividing.

If n is so large that tgammal(n+1) exceeds the range of a long double (more than 1754) but not so large that the numerator overflows, then a multiplicative solution is the best you can get without a bignum library. However, you could also use

expl(lgammal(n+1) - lgammal(n-k+1) - lgammal(k+1))

which is less precise but easier to code. (Also, if the logarithm of the coefficient is useful to you, the above formula will work over quite a large range of n and k. Not having to use expl will improve the accuracy.)

If you need a range of binomial coefficients with the same value of n, then your best bet is iterative addition:

void binoms(unsigned n, long double* res) {
  // res must have (n+3)/2 elements
  res[0] = 1;
  for (unsigned i = 2, half = 0; i <= n; ++i) {
    res[half + 1] = res[half] * 2;
    for (int k = half; k > 0; --k)
      res[k] += res[k-1];
    if (i % 2 == 0)
      ++half;
  }
}

The above produces only the coefficients with k from 0 to n/2. It has a slightly larger round-off error than the multiplicative algorithm (at least when k is getting close to n/2), but it's a lot quicker if you need all the coefficients and it has a larger range of acceptable inputs.

Answer 4

This may not be what OP is looking for, but one can analytically approximate nCr for large n with binary entropy function. It is mentioned in

• Page 10 of http://www.inference.phy.cam.ac.uk/mackay/itprnn/ps/5.16.pdf

• https://math.stackexchange.com/questions/835017/using-binary-entropy-function-to-approximate-logn-choose-k