Normal distributed rand

nobody somebody at somewhere.com
Thu Apr 9 07:30:52 PDT 2009


"bearophile" <bearophileHUGS at lycos.com> wrote in message 
news:grjgtt$2uti$1 at digitalmars.com...
> nobody:
>> Does anyone know of a normal/gaussian distributed random number generator
>> for D1 + phobos?
>
> I have created this for you, not much (unit)tested yet, adapted from 
> Python random module, I'll put it in my dlibs. Maybe there are ways to do 
> this in a faster way (especially if you just need an approximate normal 
> distribution).
> Also take a look at Tango source code, it too surely has a similar 
> function/method.
>
> import std.math: sin, cos, log, sqrt, PI, isnan;
> import std.random: rand;
>
> double random() {
>    return rand() / cast(double)uint.max;
> }
>
> /*************************************
> Normal distribution.
> mu is the mean, and sigma is the standard deviation. Not thread safe.
> */
> double normal(double mu=0.0, double sigma=1.0) {
>    // When x and y are two variables from [0, 1), uniformly
>    // distributed, then
>    //
>    //    cos(2*pi*x)*sqrt(-2*log(1-y))
>    //    sin(2*pi*x)*sqrt(-2*log(1-y))
>    //
>    // are two *independent* variables with normal distribution
>    // (mu = 0, sigma = 1).
>    // (Lambert Meertens)
>    static double gauss_next; // nan
>    auto z = gauss_next;
>    gauss_next = double.init; // nan
>    if (isnan(z)) {
>        double x2pi = random() * PI * 2;
>        double g2rad = sqrt(-2.0 * log(1.0 - random()));
>        z = cos(x2pi) * g2rad;
>        gauss_next = sin(x2pi) * g2rad;
>    }
>    return mu + z * sigma;
> }
>
>
> //--------------------------------
>
> import std.string: repeat;
> import std.stdio: put = writef, putr = writefln;
>
> void main() {
>    auto bins = new int[30];
>
>    for (int i; i < 10000; i++) {
>        auto r = cast(int)normal(bins.length / 2, 5);
>        if (r < 0)
>            r = 0;
>        if (r > (bins.length - 1))
>            r = bins.length - 1;
>        bins[r]++;
>    }
>
>    foreach (count; bins)
>        putr(">", "*".repeat(count / 20));
> }
>
> Bye,
> bearophile

Oh wow, thanks!

I think this will suit my purposes just fine.
And the main output is neat :) 




More information about the Digitalmars-d-learn mailing list