 # Boost :

From: Joel Eidsath (jeidsath_at_[hidden])
Date: 2007-12-22 18:01:57

The poisson_distribution code will fail for any mean larger than 750
or so (on my x86_64 machine). The problem is a comparison with
exp(-_mean) which evaluates as zero right around there.

The easiest fix is to operate with logs instead of exponents. I
changed the operator() code to be something like this:

RealType product = RealType(0);
for(result_type m = 0; ; ++m) {
product += log(eng());
if(product <= -_mean)
return m;
}

This also makes it possible to get rid of the init() function and the
_exp_mean private member variable.

A far better fix would be to use an algorithm other than Knuth's for
this. Knuth's algorithm is simple, but O(n). Here is some code I
wrote, converting some old Fortran code that a friend had got from
Numerical Recipes once upon a time. It works, but would need actual
testing before use in production code:

#include <cmath>

namespace fpoisson
{
double rand_zero_one()
{
double r = 0;
do
{
r = ((double)std::rand() /
((double)(RAND_MAX)+(double)(1)));
} while (r==0.0);
return r;
}

double gammln(const double xx)
{
const double cof = { 76.18009173, -86.50532033, 24.01409822,
-1.231739516,0.12085003e-2,-0.536382e-5};
const double stp = 2.50662827465;
double x = xx - 1;
double tmp = x + 5.5;
tmp = (x+0.5)*log(tmp)-tmp;
double ser = 1;
for(int j=0;j<6;++j)
{
x+=1;
ser+=cof[j]/x;
}
return tmp+std::log(stp*ser);
}

double prand(const int xm)
{
const double pi=3.141592654;
double em;
double t;
if(xm<12)
{
const double g = std::exp(-xm);
em = -1.0;
double t = 1.0;
do
{
em+=1;
t=t*rand_zero_one();
} while (t>g);
return em;
}
do
{
startloop:
const double sq = sqrt(2.0*xm);
const double alxm = log(xm);
const double g = xm*alxm-gammln(xm+1);
const double y = tan(pi*rand_zero_one());
em = sq*y+xm;
if(em < 0) goto startloop;
em=static_cast<int>(em);
t=0.9*(1.0+pow(y,2))*exp(em*alxm-gammln(em+1)-g);
} while(rand_zero_one()>t);
return em;
}
}

#endif

Joel Eidsath