On 10/23/2012 04:36 PM, jerro wrote:
I have an implementation of the Ziggurat algorithm at https://github.com/jerro/phobos/blob/master/std/random.d#L2035. It modified the Ziggurat algorithm a bit, so that it doesn't need as many layers to work well, which reduces the memory consumption and makes initialization faster. The cauchy, normal and exponential functions currently use global tables, but the zigguratAlgorithm function can also be used to implement a struct or class based API.
Oh, nice! I will take a look at this in the next days.
Joseph mentions having different engines for generating normally distributed random numbers, so maybe I could write an engine for his API based on what I already have. If we had multiple engines, I don't think Ziggurat algorithm should be the default, though, because it requires relatively expensive initialization.
The caveat here is that the Ziggurat algorithm is both fast and AFAIK optimal from the statistical point of view -- yes, it's costly space-wise, but it's also safe. The Polar method might be a nice alternative that's faster than Box-Muller, that has good statistical properties, and that requires much less space. (Box-Muller actually isn't really that good statistically.)
The choice probably needs to depend on the typical number of normal variates we think people are going to want to be generating. The way my code is set up you should only need to allocate one engine per thread, so it seems like a once-off hit that's trivial in the context of any moderately sized simulation.
