Re: [boost] normal_distribution.hpp

From: Jens Maurer (Jens.Maurer_at_[hidden])
Date: 2001-10-01 14:56:26

Hello!

Levente Farkas wrote:
> I'm just look into the normal_distribution.hpp and see it's use sin and
> cos for calculating the result (use the Box-Muller Method), altough it
> can be "faster" with the Polar Method (see eg:
> http://www.dspguru.com/howto/tech/wgn2.htm)

There are quite a few papers dealing this issue.  The Polar Method
is not even the fastest available, as far as I know.

> the strange thing that the current version are faster even it use
> trigonometric functions(?) at least on a PIII with MSVC 6sp4 when
> I compile release (optimezed) code.

Strange.

> (I append the patch which I made for the test at the end of this mail).

Looks ok at first sight.

The reason for not having "advanced" methods in Boost.Random is two-fold:
 - When Boost.Random was discussed on the boost list, some people
felt that giving lots of choice to the user by providing several
implementations of a given interface (normal distribution) would
confuse non-expert users too much.  To keep things simple, only
the most basic implementations were taken.
 - Having a loop in the implementation makes using Boost.Random
difficult to use for real-time systems, since you never know
beforehand how many times you'll spin in the loop (after all, we're
taking about *random* numbers).

For those reasons, I won't replace the current, working code with
a different algorithm.

However, do feel free to propose an additional file
boost/random/polar_normal_distribution.hpp with the appropriate
code.  (If you go that way, please propose the fastest implementation
you can find a paper about, no need to go only half the way.)

> I can't find the reason. but another intersting question if the
> current implementation use the Box-Muller Method while don't have
> a do loop:
> do {
>   _r2 = _rng();
> } while (_r2 == 0);

For the loop, see above.

> and what's the reason of the
> _cached_rho = sqrt(-2 * log(1.0-_r2));
> and not the
> _cached_rho = sqrt(-2 * log(r2));

Those items are related:  If you have the latter log() code, you also
need your suggestion for computing _r2.  With that, _r2 = 1.0 cannot
ever happen (because _rng() returns [0,1[  ).   You've just, at
least in principle, broken your distribution function.

Jens Maurer