Bug in density calculation std::piecewise_constant_distribution?

Question

It seems that std::piecewise_constant_distribution computes the densities wrongly , at least with GCC and its standard library.

According to http://www.cplusplus.com/reference/random/piecewise_constant_distribution/: The densities should be computed as:

Checking this manually reveals the bug!

This can be seen here: http://coliru.stacked-crooked.com/a/ca171bf600b5148f

The source code related to this is found in /usr/include/c++/4.8/bits/random.tcc (on linux) and the extract of the initialization function _M_initialize called by the constructor shows that there is something incorrect here:

const double __sum = std::accumulate(_M_den.begin(),
                       _M_den.end(), 0.0);

      __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
                __sum);  <----- WRONG

      // THIS is not the cummulative distribution (since the above normalization does not give the probability of the intervalls!)
      _M_cp.reserve(_M_den.size());
      std::partial_sum(_M_den.begin(), _M_den.end(),
               std::back_inserter(_M_cp));


      // Make sure the last cumulative probability is one.
      _M_cp[_M_cp.size() - 1] = 1.0;



      // Dividing here by the interval length is WRONG!!!

      for (size_t __k = 0; __k < _M_den.size(); ++__k)
           _M_den[__k] /= _M_int[__k + 1] - _M_int[__k];

Answer 1

Here's the applicable part of the specification, straight from N4296:

As can be seen clearly, the summation applies only to the weights.

It's easy to see that there's something wrong with your testing code. Reducing the number of intervals to two, the first of length 1 and the second of length 2:

std::array<PREC,3> intervals {0, 1, 3};

and giving each interval weight equal to its length:

std::array<PREC,2> weights   {1, 2};

One would expect the density to be constant. But your code reports:

Probability : 0.200000000000000011102230246252
Probability : 0.400000000000000022204460492503