http://gcc.gnu.org/bugzilla/show_bug.cgi?id=55047
Bug #: 55047
Summary: operator() in std::exponential_distribution may call
log(0)
Classification: Unclassified
Product: gcc
Version: 4.6.3
Status: UNCONFIRMED
Severity: normal
Priority: P3
Component: libstdc++
AssignedTo: unassig...@gcc.gnu.org
ReportedBy: hyou...@google.com
The implementation of operator() for std::exponential_distribution is:
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
__aurng(__urng);
return -std::log(__aurng()) / __p.lambda();
}
If I understand this correctly, __aurng() returns a value in [0,1). This
leaves the possibility of computing -log(0)/lambda, which I expect to be Inf.
On the other hand, -log(1)/lambda will never occur, so the resulting value can
never be 0. There are two problems with this implementation:
1. The actual range (0,Inf) U Inf is not consistent with the reported range of
[0,Inf) as computed by the min() and max() member functions.
2. -log(U01)/lambda is not the mathematically correct form for the inverse
transformation for the exponential distribution. -log(1 - U01)/lambda is the
correct form. This form also gives you the correct range of [0,Inf). It is an
incorrect optimization to change 1-U01 to just U01 when U01 is [0,1). It is
only correct if U01 is [0,1] or (0,1), but I do not believe that to be the case
here.
I believe the correct implementation should have the following return
statement:
return -std::log(result_type(1) - __aurng()) / __p.lambda();
The same problem appears in several other distributions. For example,
std::weibull_distribution.
