how to compute binomial distribution without overflow?

[graciawang]

I need to compute the value of binomial(n, k)=n!/k!(n-k)! * (1-p)^n * p^k.

When n and k is very big (e.g. n > 160), the step to compute n!/k!(n-k)! is always overflow even when I used unsigned long long.

Can anyone help me?

Thank you

unsigned long long com (int n, int m)
{
unsigned long long cnm = 1UL;
int i, f;

if (m*2 >n) m = n-m;
for (i=1 ; i <= m; n--, i++)
{
if ((f=n) % i == 0)
f /= i;
else cnm /= i;
cnm *= f;
}
return cnm;
}

main ()
{
int n=180, k=164;

cout<<com(n, k)*pow(1-p, k)*pow(p, n);

}


Gracia

[ropez]

Can you post the code you are currently using to compute (n!/k!(n-k)!)?

This is a draft suggestion that might help you, it's not tested:

typedef unsigned long long ull;

ull bin(ull n, ull k) {
    ull b = 1;
    for (ull i = 0; i < k; ++i) {
        b *= (n - i) / (i + 1);
    }
    return b;
}

[ropez]

Forget about the version I posted, it doesn't work.

I think your com(n, m) function is correct, but 180C164 is simply too large for unsigned long long. If you're only interested in the final result, you might wan't to use double instead. The final result doesn't need com(n, m) to be precise.

Try this:

#include <cstdlib>
#include <cmath>
#include <iostream>
using namespace std;


double com(int n, int m)
{
    double cnm = 1.0;
    
    if (m * 2 > n) 
        m = n - m;
    
    for (int i = 1; i <= m; n--, i++)
    {
        cnm /= i;
        cnm *= n;
    }
    return cnm;
}

int main ()
{
    double p = 0.5;
    int n=180, k=164;
    
    cout << com(n, k) << endl;
    cout << com(n, k)*pow(1-p, k)*pow(p, n) << endl;
}

If you really need to calculate 180C164 precisely, you need to use a big integer package. You might be interested in the GNU MP Bignum Library, http://gmplib.org/

Hope this helps.