Re: [boost] Bugs in boost random number distributions

From: Jens Maurer (Jens.Maurer_at_[hidden])
Date: 2006-08-14 17:00:11

(Looks like I forgot to send this answer to the list,
instead I sent to Charles only.  Sorry.)

Charles Karney wrote:
 >> (1) uniform_int<char>(low, high) produces values outside [low, high].
 >> E.g.,
 >>
 >>    variate_generator<random_generator&, uniform_int<char> >
 >>      unichar(r, uniform_int<char>(-100, 100));

As a nit-pick, you should never ever perform arithmetic with
plain "char"s, because it's your compiler's (nonportable) choice
whether they're signed or unsigned.  The example is valid with
"signed char", though.

 >> (2) uniform_int<short>(low, high) produces values outside [low, high].
 >> E.g.,
 >>
 >>    variate_generator<random_generator&, uniform_int<short int> >
 >>      unishort(r, uniform_int<short int>(-20000, 20000));
 >
 >
 >
 > These are both caused by overflow in the statement
 >
 >     _range = _max - _min;

Yes, it's a bug.  Fixing it requires somewhat extended surgery
in this area of the code.  I've checked that into CVS HEAD, but
I haven't tested it extensively, yet.  With the fix, the "signed char"
example delivers reasonable output.

 >> (3) uniform_int<long long>(low, high) produces a distribution with the
 >> wrong mean. E.g.,

"long long" is not (yet) C++, thus I don't consider fixing this
high priority.

 > These are both caused by problems in do_compare<false, true> in
 > random/detail/signed_unsigned_compare.hpp.  In particular the statement

Since signed_unsigned_compare.hpp was inadequate considering
your negative number / overflow issues above, I'm now using something
else.

Probably this fixes that issue as well?  (I won't test that tonight.)

As for using uniform_smallint<>, you're probably right, and the
algorithm should be fixed to use rejection as a last resort only.

 >> (5) uniform_01<random_generator, float> uni01(r) sometimes returns a
 >> result outside [0,1).  Specifically it can return 1.0f because
 >> (232-1)/232 rounds to 1.

 > uniform_01 returns
 >
 >     result_type(_rng() - _rng.min()) /
 >         (result_type(_rng.max() - _rng.min()) + 1)
 >
 > This has two problems.  If result_type is float, the result can overflow
 > to 1 (for mt19937).

This is a bug.  I believe Walter Brown's paper contains a "round down"
solution that might be workable.  For the time being, I've implemented
a rejection loop for the rare case of 1.

 >   If result_type is double or long double, the result
 > only contains 32 bits of randomness (instead of 53 bits or 64 bits).

This is, from a standardese point of view, a "quality of implementation"
issue in my opinion.  Yes, the code could be improved.  No, I haven't
done that yet.

 > I might even suggest going further and requiring that max = 2^b - 1, so
 > that a whole number of bits of randomness are returned.  This greatly
 > simplifies the computation of the basic integer and real distributions.

Yes, but it excludes whole classes of engines.  The (valid) question
is whether those engines (linear congruential, inversive congruential,
etc.) are still "interesting", useful and desirable giving the
existence of mt19937.

All fixes checked in to CVS HEAD, please have a look.

Thanks for your debugging.

Jens Maurer