+ /// @cond show_private
+
+ template<class CharT, class Traits>
+ void read(std::basic_istream<CharT, Traits>& is) {
+ param_type parm;
+ if(is >> parm) {
+ param(parm);
+ }
+ }
+
+ bool use_inversion() const
+ {
+ return _mean < 10;
+ }
+
+ static RealType log_factorial(IntType k)
+ {
+ BOOST_ASSERT(k >= 0);
+ BOOST_ASSERT(k < 10);
+ return detail::poisson_table<RealType>::value[k];
+ }
+
+ void init()
+ {
+ using std::sqrt;
+ using std::exp;
+
+ if(use_inversion()) {
+ _exp_mean = exp(-_mean);
+ } else {
+ _ptrd.smu = sqrt(_mean);
+ _ptrd.b = 0.931 + 2.53 * _ptrd.smu;
+ _ptrd.a = -0.059 + 0.02483 * _ptrd.b;
+ _ptrd.inv_alpha = 1.1239 + 1.1328 / (_ptrd.b - 3.4);
+ _ptrd.v_r = 0.9277 - 3.6224 / (_ptrd.b - 2);
+ }
+ }
+
+ template<class URNG>
+ IntType generate(URNG& urng) const
+ {
+ using std::floor;
+ using std::abs;
+ using std::log;
+
+ while(true) {
+ RealType u;
+ RealType v = uniform_01<RealType>()(urng);
+ if(v <= 0.86 * _ptrd.v_r) {
+ u = v / _ptrd.v_r - 0.43;
+ return static_cast<IntType>(floor(
+ (2*_ptrd.a/(0.5-abs(u)) + _ptrd.b)*u + _mean + 0.445));
+ }
+
+ if(v >= _ptrd.v_r) {
+ u = uniform_01<RealType>()(urng) - 0.5;
+ } else {
+ u = v/_ptrd.v_r - 0.93;
+ u = ((u < 0)? -0.5 : 0.5) - u;
+ v = uniform_01<RealType>()(urng) * _ptrd.v_r;
+ }
+
+ RealType us = 0.5 - abs(u);
+ if(us < 0.013 && v > us) {
+ continue;
+ }
+
+ RealType k = floor((2*_ptrd.a/us + _ptrd.b)*u+_mean+0.445);
+ v = v*_ptrd.inv_alpha/(_ptrd.a/(us*us) + _ptrd.b);
+
+ RealType log_sqrt_2pi = 0.91893853320467267;
+
+ if(k >= 10) {
+ if(log(v*_ptrd.smu) <= (k + 0.5)*log(_mean/k)
+ - _mean
+ - log_sqrt_2pi
+ + k
+ - (1/12. - (1/360. - 1/(1260.*k*k))/(k*k))/k) {
+ return static_cast<IntType>(k);
+ }
+ } else if(k >= 0) {
+ if(log(v) <= k*log(_mean)
+ - _mean
+ - log_factorial(static_cast<IntType>(k))) {
+ return static_cast<IntType>(k);
+ }
+ }
+ }
+ }
+
+ template<class URNG>
+ IntType invert(URNG& urng) const
+ {
+ RealType p = _exp_mean;
+ IntType x = 0;
+ RealType u = uniform_01<RealType>()(urng);
+ while(u > p) {
+ u = u - p;
+ ++x;
+ p = _mean * p / x;
+ }
+ return x;
+ }
+
+ RealType _mean;
+
+ union {
+ // for ptrd
+ struct {
+ RealType v_r;
+ RealType a;
+ RealType b;
+ RealType smu;
+ RealType inv_alpha;
+ } _ptrd;
+ // for inversion
+ RealType _exp_mean;
+ };
+
+ /// @endcond