Re: C++0X uniform_int_distribution and uniform_real_distribution
jkherciueh@gmx.net wrote:
Andrei Alexandrescu (See Website For Email wrote:
C++0X defines uniform_int_distribution<some_integral> and
uniform_real_distribution<some_floating_point>. I have two questions:
1. uniform_int_distribution generates numbers in the range [a, b], while
uniform_real_distribution generates numbers in the range [a, b). Why
the difference? My reasoning is:
a. Sometimes I want the maximum number generated to be
numeric_limits<some_integral>::max(), so if the generated range were [a,
b), I'd be unable to express that due to overflow issues.
b. One rarely needs real distributions in a closed interval.
c. It's even more rare to want the right-hand side of that interval to
be numeric_limits<some_floating_point>::max().
Is this in line with the rationale of the design? (I think a and c are
strong, and b not as strong.)
There is also the issue that a uniform real X in [0,1] lies in (0,1) with
probability 1. Modelling X by a double, chances are that a good uniform RNG
will yield either of the boundary values with probability around
0.5^mantissa-length, which is 0 for all practical purposes.
Let's see... with mt19937, the code readily reveals that 1.0 is yielded
at least every 4.2 billion iterations. This may happen more often,
depending on the parameters of the generator and the choice of the word
length. Also, when the difference between a and b is small, the allowed
distinct values in the interval [a, b] can be arbitrarily few. I'm not
so confident that the probability is practically nil for all applications.
The only
benefit that using a closed interval buys you is that you have to add
checks in your code whenever you divide by X or 1-X (to guard against a
division by 0 that might happen but only with probability practically 0).
If you do perform a division.
Furthermore, if one does want a random real number in the range [a, b]
(which I just found myself needing), they face writing some nontrivial
code:
template<class Urng, class some_fp>
some_fp uniform_closed_interval(Urng& g, some_fp a, some_fp b)
{
assert(a <= b);
return a + ((b - a) * (g() - g.min())) / (g.max() - g.min());
}
Parens are strategically placed such that the conversion to floating
point type occurs before the division, as it should. But there might be
overflow issues when b - a is large, so perhaps this is a safer choice:
template<class Urng, class some_fp>
some_fp uniform_closed_interval(Urng& g, some_fp a, some_fp b)
{
assert(a <= b);
some_fp t = g() - g.min();
return a + (b - a) * (t / (g.max() - g.min()));
}
Either code looks a little fishy to me. If b-a has a long and complicated
mantissa
(I understand a "long" mantissa, but when is it "complicated"?)
and g.max()-g.min() is not a power of 2 (but, say 2^32 -1), then
it is not at all clear that if g() == g.max(), the computation yields b.
You multiply b-a by something and then divide by the same number. If
rounding occurs in both steps, you could end up a little above b, a little
below b, or exactly at b.
This only strengthens my point - it's not trivial code to write.
I agree that the first version is fishy, but I can't find a fault with
the second. In the second, when g() == g.max(), then the division (which
is performed first) yields exactly 1, and then (b - a) * 1.000000
remains unchanged; barring for large differences in the values of a and
b, which would imply that a + (b - a) does not yield b. Probably it's
not worth accounting for that case.
The original Mersenne Twister has a very simple way of generating
numbers in (0, 1), [0, 1], and [0, 1) - it just divides the 32-bit
integral value by 4294967295.0 or 4294967296.0 (and to obtain (0, 1) it
first adds 0.5 to the generated value). See
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c
I hope that's correct, and I'm sure there are more efficient ways of
getting the job done. Anyhow, my point is that writing the code is not a
slam dunk and it would be nice if the library provided it, instead of
having me write it. It's a common task, so the size of the random number
library works against the library, as it emphasizes the irony of the
situation.
Could you provide an example where you need [0,1] and [0,1) would not do? By
all likelyhood, the value 1 would not be returned by the RNG in the first
case during a complete run of the application: it _should_ take a lot of
time to generate enough random numbers to get a particular double.
I don't need to provide any example. You can't provide an example that
proves that you need any one particular value in a random range every
once in a few billion iterations.
It's all a simple matter of expected interface. All I do is look up an
algorithm that entails choosing a random number in [a, b], and then code
it in C++. Then I'd naturally want to generate random numbers in [a, b]
and not [a, b), without caring about how often b will occur, or whether
it will occur at all in any given run.
Again, I'm getting back to the same theme: I don't want to become an
expert in random numbers in order to get one. In this particular
instance, I don't want to have to sit down with pen and paper and
convince myself that [a, b) will do for me when I want to encode [a, b].
Whatever task I'm up to mentions [a, b]. Then I should call a function
that yields just that.
Andrei
---
[ comp.std.c++ is moderated. To submit articles, try just posting with ]
[ your news-reader. If that fails, use mailto:std-c++@ncar.ucar.edu ]
[ --- Please see the FAQ before posting. --- ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html ]