dmurdoch@watstat.waterloo.edu (Duncan Murdoch) (01/15/91)
I've just done a little spelunking into the TP 6.0 random number generator, and
have found some things that might be of interest to others.
1. It's based on a linear congruential generator.
The system.randseed longint value is updated by multiplying as an unsigned
32 bit number by
134775813 = $8088405
and then adding 1. Overflow is ignored, so the effective formula is
newvalue := (oldvalue*134775813 + 1) mod 2^32.
2. The integer versions of random are calculated using the high word of the
randseed, after updating it.
The formula is very simple:
random(n) := hiword(newvalue) mod n
One consequence of the interaction of 1 and 2 is that if n is 2^m,
the period of random(n) is 2^(16+m). If n is not an even power of two,
the period is 2^32.
3. The extended (floating point) version of random is calculated by dividing
the absolute value of the seed by 2^31, after updating it.
This has a really nasty consequence. If you choose to calculate your own
integer valued random numbers by the formula
myrandom(n) := trunc(n*random);
you'll have random coming out to 1 once in 2^32 tries (it's the value
you'll get first if you start with randseed = -1498392781), and so
myrandom(n) will give n. If it's indexing into an array[0..n-1], this
will give an out of bounds error.
I hope this is of interest to someone!
Duncan Murdochdmurdoch@watstat.waterloo.edu (Duncan Murdoch) (01/16/91)
In article <1991Jan14.213834.26624@maytag.waterloo.edu> I wrote: > >3. The extended (floating point) version of random is calculated by dividing > the absolute value of the seed by 2^31, after updating it. > >This has a really nasty consequence. If you choose to calculate your own >integer valued random numbers by the formula > > myrandom(n) := trunc(n*random); > >you'll have random coming out to 1 once in 2^32 tries (it's the value >you'll get first if you start with randseed = -1498392781), and so >myrandom(n) will give n. If it's indexing into an array[0..n-1], this >will give an out of bounds error. I've now taken a look at the real version of random (which you get if you compile with the $N- option), and find that it doesn't suffer from this flaw. It appears to calculate the value by dividing an unsigned long version of randseed by 2^32; this means it can never return 1. It does return 0 if the seed is 0 after updating; you can get this by setting randseed to 649090867. Duncan Murdoch dmurdoch@watstat.waterloo.edu