Boost logo

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;
                        r = ((double)std::rand() /
                } while (r==0.0);
                return r;

        double gammln(const double xx)
                const double cof[6] = { 76.18009173, -86.50532033, 24.01409822,
                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)
                return tmp+std::log(stp*ser);

        double prand(const int xm)
                const double pi=3.141592654;
                double em;
                double t;
                        const double g = std::exp(-xm);
                        em = -1.0;
                        double t = 1.0;
                        } while (t>g);
                        return em;
                        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;
                } while(rand_zero_one()>t);
                return em;


Joel Eidsath

Boost list run by bdawes at, gregod at, cpdaniel at, john at