www.digitalmars.com         C & C++   DMDScript  

D - a try at reworking the Phobos random module

reply Helmut Leitner <helmut.leitner chello.at> writes:
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
May 19 2003
parent reply mjm <mjm_member pathlink.com> writes:
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)
);
Aug 14 2003
next sibling parent reply Bill Cox <bill viasic.com> writes:
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
Aug 14 2003
next sibling parent reply mjm <mjm_member pathlink.com> writes:
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.
Aug 14 2003
parent Bill Cox <bill viasic.com> writes:
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
Aug 14 2003
prev sibling parent Helmut Leitner <helmut.leitner chello.at> writes:
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
Aug 23 2003
prev sibling parent Helmut Leitner <helmut.leitner chello.at> writes:
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
Aug 23 2003