[Bug libstdc++/55047] New: operator() in std::exponential_distribution may call log(0)

[Bug libstdc++/55047] New: operator() in std::exponential_distribution may call log(0)

[hyounes at google dot com]

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: unassigned@gcc.gnu.org
        ReportedBy: hyounes@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.

[Bug libstdc++/55047] operator() in std::exponential_distribution may call log(0)

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=55047

Paolo Carlini <paolo.carlini at oracle dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |ASSIGNED
   Last reconfirmed|                            |2012-10-24
         AssignedTo|unassigned at gcc dot       |paolo.carlini at oracle dot
                   |gnu.org                     |com
   Target Milestone|---                         |4.7.3
     Ever Confirmed|0                           |1

--- Comment #1 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-10-24 09:43:46 UTC ---
You are right, let's fix this asap in mainline and 4.7. I think it went that
the issue already existed many years ago when we integrated the contribution as
TR1 code and then we propagated it in more places :(

[Bug libstdc++/55047] operator() in std::exponential_distribution may call log(0)

[paolo at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=55047

--- Comment #2 from paolo at gcc dot gnu.org <paolo at gcc dot gnu.org> 2012-10-24 12:20:32 UTC ---
Author: paolo
Date: Wed Oct 24 12:20:19 2012
New Revision: 192762

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=192762
Log:
2012-10-24   Haakan Younes  <hyounes@google.com>
         Paolo Carlini  <paolo.carlini@oracle.com>

    PR libstdc++/55047
    * include/bits/random.h (exponential_distribution<>::operator):
    Fix formula to std::log(result_type(1) - __aurng()).
    * include/bits/random.tcc: Likewise, everywhere.

Modified:
    trunk/libstdc++-v3/ChangeLog
    trunk/libstdc++-v3/include/bits/random.h
    trunk/libstdc++-v3/include/bits/random.tcc

[Bug libstdc++/55047] operator() in std::exponential_distribution may call log(0)

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=55047

--- Comment #3 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-10-24 12:21:17 UTC ---
Fixed in mainline. A similar fix will go in 4_7-branch in a few days.

[Bug libstdc++/55047] operator() in std::exponential_distribution may call log(0)

[paolo at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=55047

--- Comment #4 from paolo at gcc dot gnu.org <paolo at gcc dot gnu.org> 2012-11-01 21:09:57 UTC ---
Author: paolo
Date: Thu Nov  1 21:09:51 2012
New Revision: 193070

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=193070
Log:
2012-11-01   Haakan Younes  <hyounes@google.com>
         Paolo Carlini  <paolo.carlini@oracle.com>

    PR libstdc++/55047
    * include/bits/random.h (exponential_distribution<>::operator):
    Fix formula to std::log(result_type(1) - __aurng()).
    * include/bits/random.tcc: Likewise, everywhere.

Modified:
    branches/gcc-4_7-branch/libstdc++-v3/ChangeLog
    branches/gcc-4_7-branch/libstdc++-v3/include/bits/random.h
    branches/gcc-4_7-branch/libstdc++-v3/include/bits/random.tcc

[Bug libstdc++/55047] operator() in std::exponential_distribution may call log(0)

[paolo.carlini at oracle dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=55047

Paolo Carlini <paolo.carlini at oracle dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|ASSIGNED                    |RESOLVED
         Resolution|                            |FIXED

--- Comment #5 from Paolo Carlini <paolo.carlini at oracle dot com> 2012-11-01 21:11:02 UTC ---
Fixed mainline and 4.7.3.