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
Bill Cox <bill viasic.com> 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: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