New Random Number GeneratorMark4Flies at aol.com Mark4Flies at aol.comSat Jul 27 17:52:55 UTC 2002
Years ago I requested that the original random number generator inherited
from Smalltalk-80 be updated to a modern linear congruential method. In no
time at all, it was done! This time I thought I could help with a newer
generator.
I use Squeak intermittently so when I wanted to bring up a new simulation and
realized that drawing a lot of random numbers would exhaust the current
generator, I saw an opportunity to pick Squeak up again AND finally
contribute something. The method is known as the Mersenne Twister and we use
it in our commercial software at SAS Institute. You can find everything about
it at the source:
http://www.math.keio.ac.jp/~matumoto/emt.html
I got the class to work (no syntax or logic errors) but the numbers are any
thing but uniformly distributed on [0,1) ! I think the problem is related to
how the algorithm coded in C relies on bit operations to work within unsigned
32-bit integers. I tried to constrain the elements of a WordArray (for
example, masking all intermediate results with bitAnd: 16rFFFFFFFF) but so
far it still is not cooking up good numbers.
I was hoping that someone might have some experience with operations on
arrays like this in C or in Squeak and they could give me a clue. In the mean
time I will keep experimenting to learn more about the methods in Squeak. (I
love this stuff!)
Here is my new class, followed by the original C program. (The C program can
be set up for random 32-bit integers or floating point numbers on [0,1),
easily. Just look for two sets of comments related to this behavior. Right
now it is set up for unsigned long integers.)
BIG thanks in advance for any help that you might provide.
-Mark Bailey
'From Squeak3.0 of 4 February 2001 [latest update: #3545] on 27 July 2002 at
11:27:51 am'!
Random subclass: #MTRandom
instanceVariableNames: 'mt mti '
classVariableNames: 'LowerMask M MatrixA N TemperingMaskB TemperingMaskC
UpperMask '
poolDictionaries: ''
category: 'Kernel-Numbers'!
!MTRandom methodsFor: 'supporting' stamp: 'mwb 7/27/2002 11:00'!
mask: aWord1 and: aWord2
"Perform bit logic and maintain unsigned four-word 32-bit integers"
| t1 t2 |
t1 _ (aWord1 bitAnd: UpperMask) bitAnd: 16rFFFFFFFF.
t2 _ (aWord2 bitAnd: LowerMask) bitAnd: 16rFFFFFFFF.
^ (t1 bitOr: t2) bitAnd: 16rFFFFFFFF! !
!MTRandom methodsFor: 'supporting' stamp: 'mwb 7/27/2002 11:25'!
temperingShiftL: anInteger
"Second of four Mersenne Twister shifts"
^ anInteger >> 18 bitAnd: 16rFFFFFFFF! !
!MTRandom methodsFor: 'supporting' stamp: 'mwb 7/27/2002 11:25'!
temperingShiftS: anInteger
"Second of four Mersenne Twister shifts"
^ anInteger << 7 bitAnd: 16rFFFFFFFF! !
!MTRandom methodsFor: 'supporting' stamp: 'mwb 7/27/2002 11:25'!
temperingShiftT: anInteger
"Second of four Mersenne Twister shifts"
^ anInteger << 15 bitAnd: 16rFFFFFFFF! !
!MTRandom methodsFor: 'supporting' stamp: 'mwb 7/27/2002 11:25'!
temperingShiftU: anInteger
"First of four Mersenne Twister shifts"
^ anInteger >> 11 bitAnd: 16rFFFFFFFF! !
!MTRandom methodsFor: 'supporting' stamp: 'mwb 7/27/2002 11:13'!
xor: aWord1 and: aWord2
"comment stating purpose of message"
| t1 mag01 |
mag01 _ WordArray with: 0 with: MatrixA.
t1 _ (aWord1 bitXor: aWord2 >> 1) bitAnd: 16rFFFFFFFF.
^ (t1 bitXor: (mag01 at: (aWord2 bitAnd: 1) + 1)) bitAnd: 16rFFFFFFFF! !
!MTRandom methodsFor: 'initialization' stamp: 'mwb 7/27/2002 10:26'!
initialize
"Set up space for seeds, then fill them in"
| s |
mt _ WordArray new: N.
mti _ N+1.
[s _ Time millisecondClockValue.
self sgenrand: s.
s = 0] whileTrue! !
!MTRandom methodsFor: 'accessing' stamp: 'mwb 7/27/2002 09:47'!
next
"Answer a random Float in the interval [0 to 1)."
^ (self nextValue / 16rFFFFFFFF) asFloat! !
!MTRandom methodsFor: 'accessing' stamp: 'mwb 7/26/2002 17:37'!
nextInt: anInteger
"Answer a random integer in the interval [1, anInteger]."
^ (self next * anInteger) truncated + 1! !
!MTRandom methodsFor: 'private' stamp: 'mwb 7/27/2002 11:16'!
genrand
| y |
self mti >= N
ifTrue: [self mti = (N + 1) ifTrue: [self sgenrand: 4357].
1
to: N - M
do: [:kk |
y _ self mask: (self mt at: kk) and: (self mt at: kk + 1).
self mt
at: kk
put: (self xor: (self mt at: kk + M) and: y)].
N - M + 1
to: N - 1
do: [:kk |
y _ self mask: (self mt at: kk) and: (self mt at: kk + 1).
self mt
at: kk
put: (self xor: (self mt at: kk + (M - N)) and: y)].
y _ self mask: (self mt at: N) and: (self mt at: 1).
self mt
at: N
put: (self xor: (self mt at: M) and: y).
mti _ 1].
y _ self mt at: mti.
mti _ mti + 1.
y _ (y bitXor: (self temperingShiftU: y)) bitAnd: 16rFFFFFFFF.
y _ ((y bitXor: (self temperingShiftS: y)) bitAnd: TemperingMaskB)
bitAnd: 16rFFFFFFFF.
y _ ((y bitXor: (self temperingShiftT: y)) bitAnd: TemperingMaskC)
bitAnd: 16rFFFFFFFF.
y _ (y bitXor: (self temperingShiftL: y)) bitAnd: 16rFFFFFFFF.
^ y! !
!MTRandom methodsFor: 'private' stamp: 'mwb 7/26/2002 11:38'!
mt
^ mt! !
!MTRandom methodsFor: 'private' stamp: 'mwb 7/26/2002 16:23'!
mti
^ mti! !
!MTRandom methodsFor: 'private' stamp: 'mwb 7/26/2002 16:23'!
mti: anInteger
mti _ anInteger! !
!MTRandom methodsFor: 'private' stamp: 'mwb 7/26/2002 20:56'!
nextValue
"This method generates random instances of Integer. The algorithm
is described in detail by M. Matsumoto and T. Nishimura
in 'Mersenne Twister: A 623-dimensionally equidistributed uniform
pseudorandom number generator', ACM Trans. on Modeling and
Computer Simulation Vol. 8, No. 1, January pp.3-30 1998."
^ self genrand! !
!MTRandom methodsFor: 'private' stamp: 'mwb 7/27/2002 11:25'!
seed: anInteger
seed _ anInteger.
self sgenrand: (anInteger bitAnd: 16rFFFFFFFF)! !
!MTRandom methodsFor: 'private' stamp: 'mwb 7/26/2002 17:36'!
sgenrand: anInteger
mt
at: 1
put: (anInteger bitAnd: 16rFFFFFFFF).
2
to: N
do: [:i | self mt
at: i
put: (69069 * (self mt at: i - 1) bitAnd: 16rFFFFFFFF)].
self mti: N
! !
"-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- "!
MTRandom class
instanceVariableNames: ''!
!MTRandom class methodsFor: 'instance creation' stamp: 'mwb 7/26/2002 10:54'!
new
"Create new instance of MTRandom"
^ super new initialize! !
!MTRandom class methodsFor: 'class initialization' stamp: 'mwb 7/27/2002
11:23'!
initialize
" MTRandom initialize "
"Set period parameters"
N _ 624.
M _ 397.
MatrixA _ 16r9908B0DF.
LowerMask _ 16r7FFFFFFF.
UpperMask _ 16r80000000.
"Set tempering parameters"
TemperingMaskB _ 16r9D2C5680.
TemperingMaskC _ 16rEFC60000! !
MTRandom initialize!
C program:
/* A C-program for MT19937: Real number version */
/* genrand() generates one pseudorandom real number (double) */
/* which is uniformly distributed on [0,1]-interval, for each */
/* call. sgenrand(seed) set initial values to the working area */
/* of 624 words. Before genrand(), sgenrand(seed) must be */
/* called once. (seed is any 32-bit integer except for 0). */
/* Integer generator is obtained by modifying two lines. */
/* Coded by Takuji Nishimura, considering the suggestions by */
/* Topher Cooper and Marc Rieffel in July-Aug. 1997. */
/* This library is free software; you can redistribute it and/or */
/* modify it under the terms of the GNU Library General Public */
/* License as published by the Free Software Foundation; either */
/* version 2 of the License, or (at your option) any later */
/* version. */
/* This library is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. */
/* See the GNU Library General Public License for more details. */
/* You should have received a copy of the GNU Library General */
/* Public License along with this library; if not, write to the */
/* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA */
/* 02111-1307 USA */
/* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. */
/* Any feedback is very welcome. For any question, comments, */
/* see http://www.math.keio.ac.jp/matumoto/emt.html or email */
/* matumoto at math.keio.ac.jp */
#include<stdio.h>
/* Period parameters */
#define N 624
#define M 397
#define MATRIX_A 0x9908b0df /* constant vector a */
#define UPPER_MASK 0x80000000 /* most significant w-r bits */
#define LOWER_MASK 0x7fffffff /* least significant r bits */
/* Tempering parameters */
#define TEMPERING_MASK_B 0x9d2c5680
#define TEMPERING_MASK_C 0xefc60000
#define TEMPERING_SHIFT_U(y) (y >> 11)
#define TEMPERING_SHIFT_S(y) (y << 7)
#define TEMPERING_SHIFT_T(y) (y << 15)
#define TEMPERING_SHIFT_L(y) (y >> 18)
static unsigned long mt[N]; /* the array for the state vector */
static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */
/* initializing the array with a NONZERO seed */
void
sgenrand(seed)
unsigned long seed;
{
/* setting initial seeds to mt[N] using */
/* the generator Line 25 of Table 1 in */
/* [KNUTH 1981, The Art of Computer Programming */
/* Vol. 2 (2nd Ed.), pp102] */
mt[0]= seed & 0xffffffff;
for (mti=1; mti<N; mti++)
mt[mti] = (69069 * mt[mti-1]) & 0xffffffff;
}
/* double */ /* generating reals */
unsigned long /* for integer generation */
genrand()
{
unsigned long y;
static unsigned long mag01[2]={0x0, MATRIX_A};
/* mag01[x] = x * MATRIX_A for x=0,1 */
if (mti >= N) { /* generate N words at one time */
int kk;
if (mti == N+1) /* if sgenrand() has not been called, */
sgenrand(4357); /* a default initial seed is used */
for (kk=0;kk<N-M;kk++) {
y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
}
for (;kk<N-1;kk++) {
y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
}
y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
mti = 0;
}
y = mt[mti++];
y ^= TEMPERING_SHIFT_U(y);
y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
y ^= TEMPERING_SHIFT_L(y);
/* return ( (double) y / (unsigned long)0xffffffff ); */ /* reals */
return y; /* for integer generation */
}
/* this main() outputs first 1000 generated numbers */
main(int argc, char* argv[])
{
int j;
sgenrand(4357); /* any nonzero integer can be used as a seed */
for (j=0; j<1000; j++) {
printf("%5lu\n", genrand());
if (j%8==7) printf("\n");
}
printf("\n");
}
More information about the Squeak-dev mailing list |