[boost] [random] poisson_distribution problem and fix

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:

#ifndef F_POISSON_HEADER
#define F_POISSON_HEADER
#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[6] = { 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