@@ -4592,7 +4592,10 @@ public:
4592
4592
4593
4593
template <class _IntType >
4594
4594
poisson_distribution<_IntType>::param_type::param_type(double __mean)
4595
- : __mean_(__mean)
4595
+ // According to the standard `inf` is a valid input, but it causes the
4596
+ // distribution to hang, so we replace it with the maximum representable
4597
+ // mean.
4598
+ : __mean_(isinf(__mean) ? numeric_limits<double >::max() : __mean)
4596
4599
{
4597
4600
if (__mean_ < 10 )
4598
4601
{
@@ -4610,7 +4613,7 @@ poisson_distribution<_IntType>::param_type::param_type(double __mean)
4610
4613
{
4611
4614
__s_ = _VSTD::sqrt (__mean_);
4612
4615
__d_ = 6 * __mean_ * __mean_;
4613
- __l_ = static_cast <result_type> (__mean_ - 1.1484 );
4616
+ __l_ = std::trunc (__mean_ - 1.1484 );
4614
4617
__omega_ = .3989423 / __s_;
4615
4618
double __b1_ = .4166667E-1 / __mean_;
4616
4619
double __b2_ = .3 * __b1_ * __b1_;
@@ -4627,12 +4630,12 @@ template<class _URNG>
4627
4630
_IntType
4628
4631
poisson_distribution<_IntType>::operator ()(_URNG& __urng, const param_type& __pr)
4629
4632
{
4630
- result_type __x ;
4633
+ double __tx ;
4631
4634
uniform_real_distribution<double > __urd;
4632
4635
if (__pr.__mean_ < 10 )
4633
4636
{
4634
- __x = 0 ;
4635
- for (double __p = __urd (__urng); __p > __pr.__l_ ; ++__x )
4637
+ __tx = 0 ;
4638
+ for (double __p = __urd (__urng); __p > __pr.__l_ ; ++__tx )
4636
4639
__p *= __urd (__urng);
4637
4640
}
4638
4641
else
@@ -4642,19 +4645,19 @@ poisson_distribution<_IntType>::operator()(_URNG& __urng, const param_type& __pr
4642
4645
double __u;
4643
4646
if (__g > 0 )
4644
4647
{
4645
- __x = static_cast <result_type> (__g);
4646
- if (__x >= __pr.__l_ )
4647
- return __x ;
4648
- __difmuk = __pr.__mean_ - __x ;
4648
+ __tx = std::trunc (__g);
4649
+ if (__tx >= __pr.__l_ )
4650
+ return std::__clamp_to_integral<result_type>(__tx) ;
4651
+ __difmuk = __pr.__mean_ - __tx ;
4649
4652
__u = __urd (__urng);
4650
4653
if (__pr.__d_ * __u >= __difmuk * __difmuk * __difmuk)
4651
- return __x ;
4654
+ return std::__clamp_to_integral<result_type>(__tx) ;
4652
4655
}
4653
4656
exponential_distribution<double > __edist;
4654
4657
for (bool __using_exp_dist = false ; true ; __using_exp_dist = true )
4655
4658
{
4656
4659
double __e;
4657
- if (__using_exp_dist || __g < 0 )
4660
+ if (__using_exp_dist || __g <= 0 )
4658
4661
{
4659
4662
double __t ;
4660
4663
do
@@ -4664,31 +4667,31 @@ poisson_distribution<_IntType>::operator()(_URNG& __urng, const param_type& __pr
4664
4667
__u += __u - 1 ;
4665
4668
__t = 1.8 + (__u < 0 ? -__e : __e);
4666
4669
} while (__t <= -.6744 );
4667
- __x = __pr.__mean_ + __pr.__s_ * __t ;
4668
- __difmuk = __pr.__mean_ - __x ;
4670
+ __tx = std::trunc ( __pr.__mean_ + __pr.__s_ * __t ) ;
4671
+ __difmuk = __pr.__mean_ - __tx ;
4669
4672
__using_exp_dist = true ;
4670
4673
}
4671
4674
double __px;
4672
4675
double __py;
4673
- if (__x < 10 )
4676
+ if (__tx < 10 && __tx >= 0 )
4674
4677
{
4675
4678
const double __fac[] = {1 , 1 , 2 , 6 , 24 , 120 , 720 , 5040 ,
4676
4679
40320 , 362880 };
4677
4680
__px = -__pr.__mean_ ;
4678
- __py = _VSTD::pow (__pr.__mean_ , (double )__x ) / __fac[__x ];
4681
+ __py = _VSTD::pow (__pr.__mean_ , (double )__tx ) / __fac[static_cast < int >(__tx) ];
4679
4682
}
4680
4683
else
4681
4684
{
4682
- double __del = .8333333E-1 / __x ;
4685
+ double __del = .8333333E-1 / __tx ;
4683
4686
__del -= 4.8 * __del * __del * __del;
4684
- double __v = __difmuk / __x ;
4687
+ double __v = __difmuk / __tx ;
4685
4688
if (_VSTD::abs (__v) > 0.25 )
4686
- __px = __x * _VSTD::log (1 + __v) - __difmuk - __del;
4689
+ __px = __tx * _VSTD::log (1 + __v) - __difmuk - __del;
4687
4690
else
4688
- __px = __x * __v * __v * (((((((.1250060 * __v + -.1384794 ) *
4691
+ __px = __tx * __v * __v * (((((((.1250060 * __v + -.1384794 ) *
4689
4692
__v + .1421878 ) * __v + -.1661269 ) * __v + .2000118 ) *
4690
4693
__v + -.2500068 ) * __v + .3333333 ) * __v + -.5 ) - __del;
4691
- __py = .3989423 / _VSTD::sqrt (__x );
4694
+ __py = .3989423 / _VSTD::sqrt (__tx );
4692
4695
}
4693
4696
double __r = (0.5 - __difmuk) / __pr.__s_ ;
4694
4697
double __r2 = __r * __r;
@@ -4708,7 +4711,7 @@ poisson_distribution<_IntType>::operator()(_URNG& __urng, const param_type& __pr
4708
4711
}
4709
4712
}
4710
4713
}
4711
- return __x ;
4714
+ return std::__clamp_to_integral<result_type>(__tx) ;
4712
4715
}
4713
4716
4714
4717
template <class _CharT , class _Traits , class _IntType >
0 commit comments