The next version of the Venus library will include a reworked random module.

It is centered around an (part of module venus.interfaces):

  interface iRandomGenerator
  {
      uint RetUint();
      void SetSeed(uint seed);
  }

and the original phobos module is refactored to allow multiple instances and the replacement of the generator:

    /*
    ** random1.d
    ** Copyright (C) 2003 Digital Mars
    ** Reworked from Phobos, donated to Digital Mars.
    ** Released under the GNU GPL, see LICENSE.TXT
    ** and http://www.gnu.org/copyleft/gpl.html
    */

    module venus.random1;
    import venus.all;

    class RandomGenerator_Phobos: iRandomGenerator
    {
        static uint static_xormix1[20] =
        [
            0xbaa96887, 0x1e17d32c, 0x03bcdc3c, 0x0f33d1b2,
            0x76a6491d, 0xc570d85d, 0xe382b1e3, 0x78db4362,
            0x7439a9d4, 0x9cea8ac5, 0x89537c5c, 0x2588f55d,
            0x415b5e1d, 0x216e3d95, 0x85c662e7, 0x5e8ab368,
            0x3ea5cc8c, 0xd26a0f74, 0xf3a9222b, 0x48aad7e4
        ];

        static uint static_xormix2[20] =
        [
            0x4b0f3b58, 0xe874f0c3, 0x6955c5a6, 0x55a7ca46,
            0x4d9a9d86, 0xfe28a195, 0xb1ca7865, 0x6b235751,
            0x9a997a61, 0xaa6e95c8, 0xaaa98ee1, 0x5af9154c,
            0xfc8e2263, 0x390f5e8c, 0x58ffd802, 0xac0a5eba,
            0xac4874f6, 0xa9df0913, 0x86be4c74, 0xed2c123b
        ];

        uint seed;      // starting seed
        uint index;     // ith random number
        uint xormix1[20];
        uint xormix2[20];

        this() {
            xormix1[]=static_xormix1[];
            xormix2[]=static_xormix2[];
        }

        uint RetUint() {
            uint hiword, loword, hihold, temp, itmpl, itmph, i;

            loword = seed;
            hiword = index++;
            for (i = 0; i < 4; i++)     // loop limit can be 2..20, we choose 4
            {
                hihold  = hiword;                           // save hiword for later
                temp    = hihold ^  xormix1[i];             // mix up bits of hiword
                itmpl   = temp   &  0xffff;                 // decompose to hi & lo
                itmph   = temp   >> 16;                     // 16-bit words
                temp    = itmpl * itmpl + ~(itmph * itmph); // do a multiplicative mix
                temp    = (temp >> 16) | (temp << 16);      // swap hi and lo halves
                hiword  = loword ^ ((temp ^ xormix2[i]) + itmpl * itmph); //loword mix
                loword  = hihold;                           // old hiword is loword
            }
            return hiword;
        }
        void SetSeed(uint seed) {
            this.seed = seed;
            this.index = seed*33;
        }
    }

It still seems to be not good enough because I don't think that the state can be reset by the SetSeed method.

A number of methods are separated from the generator into an RandomProducer class:

    /*
    ** random2.d
    ** Copyright (C) 2003 Helmut Leitner
    ** Donated to Digital Mars.
    ** Released under the GNU GPL, see LICENSE.TXT
    ** and http://www.gnu.org/copyleft/gpl.html
    */

    module venus.random2;
    import venus.all;

    RandomGenerator_Phobos RandomGeneratorDefault;
    RandomProducer Random;

    static this()
    {
        RandomGeneratorDefault=new RandomGenerator_Phobos;
        Random=new RandomProducer;
    }

    class RandomProducer {

        iRandomGenerator rg;

        this()
        {
            rg=RandomGeneratorDefault;
        }

        this(iRandomGenerator rg)
        {
            this.rg=rg;
        }

        uint RetUint()
        {
            return rg.RetUint();
        }

        int RetInt()
        {
            return (int) rg.RetUint();
        }

        double RetDouble_01()
        {
            return rg.RetUint()/(uint.max+(double)1);
        }

        double RetDouble_11()
        {
            return (RetInt()+0.5)/(-(double)int.min);
        }

        int RetIntCount(int count)
        {
            return RetDouble_01()*count;
        }

        int RetIntRange(int lo,int hi)
        {
            return lo+RetIntCount(hi-lo+1);
        }

        void CvtRandom() {
            uint seed=(uint)TimerRetCount()&0xFFFFFFFF;
            rg.SetSeed(seed);
        }

        void SetSeed(uint seed)
        {
            rg.SetSeed(seed);
        }

        bit state=0;
        double v2;
        double fact;

        double RetDouble_Gauss() /* Box-Mueller */
        {
            double r;
            double v1;

            if(state)
            {
                state=0;
                return v2*fact;
            }
            do
            {
                v1=RetDouble_11();
                v2=RetDouble_11();
                r=v1*v1+v2*v2;
            } while(r>=1.0 || r==0.0);
            fact=sqrt(-2.0*log(r)/r);

            state=1;
            return v1*fact;
        }

    }

You can use it out of the Box like in this sample program:

    import venus.all;

    int main()
    {
      StatSum ss=new StatSum;
      double val;

      for(int i=0; i<100000; i++) {
        val=Random.RetDouble_Gauss();
        ss.AddVal(val);
      }
      PrintLine("mean=",ss.RetMean());
      PrintLine("deviation=",ss.RetDeviation());

      return 0;
    }

With that StatSum module easing the calculation of mean and standard deviation:

    /*
    ** statsum.d
    ** Copyright (C) 2003 Helmut Leitner
    ** Donated to Digital Mars.
    ** Released under the GNU GPL, see LICENSE.TXT
    ** and http://www.gnu.org/copyleft/gpl.html
    */

    module venus.statsum;
    import venus.all;

    class StatSum {
      double sum=0.0;
      double sum_squares=0.0; /* sum of val*val */
      int n=0;

      void Clear() {
        sum=0.0;
        sum_squares=0.0;
        n=0;
      }

      void AddVal(double x) {
        sum+=x;
        sum_squares+=(x*x);
        n+=1;
      }

      void AddValCount(double x,int n) {
        sum+=(n*x);
        sum_squares+=(n*x*x);
        n+=n;
      }

      double RetMean() {
        double mean;
        if(n<1) {
          mean=0.0;
        } else {
          mean=sum/n;
        }
        return mean;
      }

      error GetDeviation(double *pdeviation) {
        if(n<2) {
          *pdeviation=0.0;
          return -1;
        }
        *pdeviation=sqrt((sum_squares - sum*sum/n)/(n-1));
        return 0;
      }

      double RetDeviation()
      {
        double deviation;
        if(n<2) {
          deviation=0.0;
        } else {
          deviation=sqrt((sum_squares - sum*sum/n)/(n-1));
        }
        return deviation;
      }

    }

which will also be in Venus 0.02.

I would be happy about any comments or suggestions for improvement.

--
Helmut Leitner    leitner@hls.via.at Graz, Austria   www.hls-software.com

If you are interested in probability take a look at
http://matringale.berlios.de/Martingale.html.
I have some low discrepancy sequence generators there which could be moved from
C++ to D.

Also I think that the Box-Muller algorithm could be replaced with the Inverse Normal Distribution Function (Random.h,cc: nInverse):

Next Gaussian deviate =
InverseCumulativeNormalDistibutionFunction(
next uniform deviate in (0,1)
);

mjm wrote:
> If you are interested in probability take a look at http://matringale.berlios.de/Martingale.html.
> I have some low discrepancy sequence generators there which could be moved from
> C++ to D.
> 
> Also I think that the Box-Muller algorithm could be replaced with the Inverse
> Normal Distribution Function (Random.h,cc: nInverse):
> 
> Next Gaussian deviate = InverseCumulativeNormalDistibutionFunction(
> next uniform deviate in (0,1)
> );
> 
> 
> 

Isn't the Mersenne Twister algorithm faster?  It also claims to pass the Die Hard tests, which is good enough for me.  It's the default generator in the GNU Scientific Library.

Here's a link to the guy's web site:

http://www.math.keio.ac.jp/matumoto/emt.html

You can use his code for free.  It should be very easy to convert to D.

-- Bill

>Isn't the Mersenne Twister algorithm faster?  It also claims to pass the Die Hard tests, which is good enough for me.  It's the default generator in the GNU Scientific Library.
>
>Here's a link to the guy's web site:
>
>http://www.math.keio.ac.jp/matumoto/emt.html
>
>You can use his code for free.  It should be very easy to convert to D.
>
>-- Bill
>

The Mersenne twister is the best  currently known UNIFORM random number
generator
and it's also very fast. It's definitely something that we want to have in D.

From the Mersenne twister you get uniform (equally likeley) random integers (in
some fixed interval) which you then can turn into random rationals in (0,1).

The Box-Muller and the Inverse Cumulative Normal Distribution Function then make standard normal deviates out of uniform deviates.

mjm wrote:
>>Isn't the Mersenne Twister algorithm faster?  It also claims to pass the Die Hard tests, which is good enough for me.  It's the default generator in the GNU Scientific Library.
>>
>>Here's a link to the guy's web site:
>>
>>http://www.math.keio.ac.jp/matumoto/emt.html
>>
>>You can use his code for free.  It should be very easy to convert to D.
>>
>>-- Bill
>>
> 
> 
> The Mersenne twister is the best  currently known UNIFORM random number
> generator
> and it's also very fast. It's definitely something that we want to have in D.
> 
> From the Mersenne twister you get uniform (equally likeley) random integers (in
> some fixed interval) which you then can turn into random rationals in (0,1).
> 
> The Box-Muller and the Inverse Cumulative Normal Distribution Function then make
> standard normal deviates out of uniform deviates.

Ok.  Thanks for the info.

Bill

mjm wrote:
> 
> If you are interested in probability take a look at
> http://matringale.berlios.de/Martingale.html.
> I have some low discrepancy sequence generators there which could be moved from
> C++ to D.
> 
> Also I think that the Box-Muller algorithm could be replaced with the Inverse Normal Distribution Function (Random.h,cc: nInverse):
> 
> Next Gaussian deviate =
> InverseCumulativeNormalDistibutionFunction(
> next uniform deviate in (0,1)
> );

Thank you for hinting to
  http://martingale.berlios.de/Martingale.html
(repeated for typo)

Typically I would offer different implementations in a module if they have relative merits in a way that they don't bloat the executables.

--
Helmut Leitner    leitner@hls.via.at Graz, Austria   www.hls-software.com

Bill Cox wrote:
> 
> mjm wrote:
> > If you are interested in probability take a look at
> > http://matringale.berlios.de/Martingale.html.
> > I have some low discrepancy sequence generators there which could be moved from
> > C++ to D.
> >
> > Also I think that the Box-Muller algorithm could be replaced with the Inverse Normal Distribution Function (Random.h,cc: nInverse):
> >
> > Next Gaussian deviate =
> > InverseCumulativeNormalDistibutionFunction(
> > next uniform deviate in (0,1)
> > );
> >
> >
> >
> 
> Isn't the Mersenne Twister algorithm faster?  It also claims to pass the Die Hard tests, which is good enough for me.  It's the default generator in the GNU Scientific Library.
> 
> Here's a link to the guy's web site:
> 
> http://www.math.keio.ac.jp/matumoto/emt.html

Thank you for the link.

It seems Martingale also uses the Mersenne Twister. I intend to included perhaps about 5 different generators for different purposes. Sometimes you need maximum quality (physical simulations), sometimes maximum speed (mutation algorithms)...

--
Helmut Leitner    leitner@hls.via.at Graz, Austria   www.hls-software.com