implementing Poisson distribution in c++

Question

I'm trying to write a program to calculate probabilty mass function of a poisson distribution, P(x=n) with parameter lambda, using this formula: ( (e^-lambda)*(lambda^n))/n!

This approach works well when I use small lambda and small numbers, but if I want to calculate for example P(x=30) with lambda 20 the result is 4.68903e+006 which is wrong.

I think the problem is for calculating n!. I implemented a function to calculate factorial value and used unsigned long long data type for the result of factorial calculation, but the problem is that the amount of 30! is equal to 265,252,859,812,191,058,636,308,480,000,000 and the maximum number available for unsigned long long is 18,446,744,073,709,551,615, which is less than 30!.

What should I do to handle this issue? Is there any other way or any function to calculate this prabability in c++?

data type

Answer 1

One workaround to deal with large n's is calculating the distribution in the log domain:

X = ((e^-lambda)*(lambda^n))/n!
ln X = -lambda + n*ln(lambda) - Sum (ln(n))
return e^X