random number gen.
cmo...@uoft02.utoledo.edu
a while back, I thought this might be of interest. It implements
the multiply-with-carry algorithm recommended by Marsaglia ( a guru
of RNGs ). So far as I know, it hasn't been formally published,
but there is a discussion in the numerical analysis FAQs.
The essence of the algorithm is
x[n+1] = x[n] * multiplier + carry mod 2^32
where carry is taken from the previous iteration. I.e., the
top 32 bits of the 64-bit product+sum is used to add to the next
product. (It's clearer in Forth.)
This is written to generate 16-bit (unsigned) integers with a 16-bit
Forth, or 32-bit values with a 32-bit Forth. The periods should be
of the order of 2^29 for 16 bits and 2^58 for 32 bits. If you want
32 bits with the long period on a 16-bit system, just run separate
16-bit generators with different multipliers for the upper and lower
halves.
On an EIGHT-bit system, you can run four 8-bit generators, all
with different multipliers; with multipliers 65FF 56FF 68FF and 317F
the 32-bit values should have a period of about 2^57.5 (~2x10^17).
Good multiplier values are important, although for the 32-bit
version just about any 30-bit multiplier is apt to give a period
of 2^50 or so.
\ <snip>
\ multiply-with-carry pseudo-random number generator cgm Feb1996
\ assumes address units are 8-bit bytes, and a cell is either 16 or 32 bits
HEX
1 CELLS 2 = [IF] 65E8 [ELSE] 7A3FFD4B [THEN]
CONSTANT rmult
\ other good multipliers for 16-bit systems are 61BF 62DC 6594 6363 5E9B
\ for 32-bits 7549D83B 70DCB7DC 6A364C22 62AFDC7A
2VARIABLE rloc 3 1 rloc 2! \ any seed values should be ok
: rndm ( -- u ) \ returns one cell full of random bits
rloc 2@ rmult UM* ROT 0 D+ OVER rloc 2! ;
\ That was pretty easy, wasn't it! ;-)
\ a few related words that might be useful:
\ random integer from 0 to n-1
: random ( n -- 0..n-1 ) rndm UM* NIP ;
\ floating point value uniformly distributed in [0-1)
0 1 D>F FCONSTANT f2^max
: frand ( -- f ) rndm 0 D>F f2^max F/ ;
\ reseed randomly by exercising rndm a few times
: randomize ( -- ) TIME&DATE 2DROP ( sec min hr day )
+ + + 0 DO rndm DROP LOOP
rndm rndm rloc 2! ;
\ <snip>
enjoy... cgm
C. G. Montgomery c...@physics.utoledo.edu cmo...@uoft02.utoledo.edu