Re: [boost] Bugs in boost random number distributions

From: Charles Karney (ckarney_at_[hidden])
Date: 2006-08-08 15:05:41

Neal Becker wrote:
> Charles Karney wrote:
>
>>* a clumsy (and inexact) interface between the random generators and
>>  the random distributions.
>
> I wonder what you refer to here?

Here's one example.  In bernoulli_distribution we have

    if(_p == RealType(0))
      return false;
    else
      return RealType(eng() - (eng.min)()) <=
             _p * RealType((eng.max)()-(eng.min)());

where eng() returns integers between eng.min() and eng.max() incl.  The
problems are:

* If the integers returned by eng() are signed (e.g.,
  linear_congruential) then eng() - (eng.min)() and (eng.max)() -
  (eng.min)() can overflow and the results are bogus.  (Perhaps this
  can't happen because eng.min() is always zero.  Nevertheless the
  problem seems to be present in the overall design of the library.)

* The results depend on the granularity of the underlying engine.  Why
  not return an exact result?  The answer lies in the overall structure
  of the library where the bernoulli is expected to use the raw numbers
  coming out of the generator.  If there were a
  uniform_float01<RealType>(eng); available which was equivalent to
  rounding u down to the nearest representable RealType, then an EXACT
  implementation of bernoulli would be

    return p < uniform_float01<RealType>(eng);

  But typically we can write a faster implementation because it's
  usually unnecessary to generate all the bits of
  uniform_float01<RealType>(eng) before being able to do the comparison.

> I have been thinking that the use of a member template for
>
>   template<class Engine>
>   result_type operator()(Engine& eng)
>
> is sometimes awkward.  It doesn't allow compile-time inquiry of the
> engine properties by the distribution class.  For example, I have a
> distribution that generates an integer of random bits of a fixed
> width.  It is awkward that the distribution class can't get the engine
> bit width.

With boost::random there's no such thing as "engine bit width", e.g.,
ecuyer1988 returns results in [1,2^31-87].  So rather than make the code
implementing the distributions pay the penalty for extracting a whole
number of random bits, it would make sense to push the responsibility
off onto the generator (or another intermediate layer).  Then, you could
ask for an unsigned "word" of bits and you would need to be able to
query the generator for the width of this word (as a compile-time
constant).  (In the case of ecuyer1998 this width would be 30 and you'd
be paying a factor of 2 penalty in sampling until you got a result in
range.  Presumably for this application you would want to devise a
L'Ecuyer-style generator with a max value of 2^30+n.  Or else you could
stick with mt19937 where you get 32 bits of randomness right out of the
box.)

-- 
Charles Karney <ckarney_at_[hidden]>
Sarnoff Corporation, Princeton, NJ 08543-5300
URL: http://charles.karney.info
Tel: [phone]
Fax: [phone]