Francesco Di Lauro - 10 months ago 54

C 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 Source

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.