[Bug libstdc++/61369] New: std::discrete_distribution::operator() may return event 0 even if its probability is 0.0

[Bug libstdc++/61369] New: std::discrete_distribution::operator() may return event 0 even if its probability is 0.0

[stefan_stuff at web dot de]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=61369

            Bug ID: 61369
           Summary: std::discrete_distribution::operator() may return
                    event 0 even if its probability is 0.0
           Product: gcc
           Version: 4.9.0
            Status: UNCONFIRMED
          Severity: trivial
          Priority: P3
         Component: libstdc++
          Assignee: unassigned at gcc dot gnu.org
          Reporter: stefan_stuff at web dot de

OLD CODE: std::discrete_distribution::operator() in 4.9.0 random.tcc:

2843  __detail::_Adaptor<_UniformRandomNumberGenerator, double>
2844  __aurng(__urng);
2845 
2846  const double __p = __aurng();
2847  auto __pos = std::lower_bound(__param._M_cp.begin(),
2848  __param._M_cp.end(), __p);
2849 
2850  return __pos - __param._M_cp.begin();

PROPOSED FIX: std::lower_bound should be std::upper_bound.

EXAMPLE: Assume the following distribution:
(The specification allows 0-weights, if at least one weight is strictly
positive.)

event:           0      |     1      |  2  |     3 
probability:    0.0     |    0.5     | 0.0 |    0.5
cummulative:    0.0     |    0.5     | 0.5 |    1.0
lower_bound: [0.0, 0.0] | (0.0, 0.5] |  -  | (0.5, 1.0)
upper_bound:     -      | [0.0, 0.5) |  -  | [0.5, 1.0)

__p is in the interval [0.0, 1.0). The lower 2 rows describe the intervals of
__p, in which the according event-index is returned. This shows, that
std::lower_bound produces wrong results, if the first element has probability
0.0 and __p is also 0.0.

(Another positive effect of using std::upper_bound is, that the intervals are
more balanced as all intervals are inclusive-exclusive, in contrast to
lower_bound, which is inclusive-inclusive for the last event.)

This will *probably* be fixed by Bug ID 57925, where this algorithm is
reworked, but I haven't done the math for the new algorithm.