## D - a try at reworking the Phobos random module

- Helmut Leitner (229/229) May 19 2003 The next version of the Venus library will include a reworked
- mjm (10/10) Aug 14 2003 If you are interested in probability take a look at
- Bill Cox (8/23) Aug 14 2003 Isn't the Mersenne Twister algorithm faster? It also claims to pass the...
- mjm (7/14) Aug 14 2003 The Mersenne twister is the best currently known UNIFORM random number
- Bill Cox (3/26) Aug 14 2003 Ok. Thanks for the info.
- Helmut Leitner (8/33) Aug 23 2003 Thank you for the link.
- Helmut Leitner (9/22) Aug 23 2003 Thank you for hinting to

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

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

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

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. -- BillThe 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

mjm wrote:Ok. Thanks for the info. BillIsn'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. -- BillThe 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

Bill Cox wrote:mjm wrote: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.comIf 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

Aug 23 2003

mjm wrote:

Aug 23 2003