dag@hp-lsd.HP.COM (David Geiser) (07/06/87)
The header pretty much says it. I'm operating on a PDP-11/73, which gives
the time only down to the second, so it can't use time as a seed (actually,
I guess you can get tenths of a second, but that still doesn't help).
.....
May vi rot in /dev/null.
----------------------------------------------------------------------
The random # stuff that follows is from a recent edition of BYTE. I think it
was done in Pascal in the article; I've re-done it in Modula-2. I think it is
pretty slick and can be adapted to many different CPUs -- that's why I'm
posting instead of E-mailing.
I had to do some "Mickey Mouse" stuff with a 'tmp' variable to work around a
compiler defect with REAL numbers. (Sorry Mickey).
If I remembered who the author(s) were, I'd give credit.... sorry.
dag
IMPLEMENTATION MODULE utils;
VAR
x, y, z : INTEGER;
x1,y1,z1: REAL;
PROCEDURE rand ( ) : REAL;
(*
Anyone who considers arithmetical methods of producing random digits is, of
course, in a state of sin.
--John von Neumann (1951)
...
Many generator use a simple multiplicative rule to generate the next number
based on a fixed prime number and a multiplier. For instance, if the prime
number were 30269 and the multiplier 171, you would have
x(i+1) = (171*x(i)) MOD 30269
where 'i' is the iteration number (dag) ...
The solution depends upon two techniques. First, consider the problem of
using a high-level language. The simple generator given above can be written
as x := (171 * x) MOD 30269. In practice this coding will not work adequately
on all machines. The values of x are restricted to the range 1 - 30268 and
therefore will fit into a 16 bit word. This is not true for 171*x, so the
expression will malfunction on many machines. You could perform the above
calculation "double length", but this either requires machine code or is very
expensive in both time and space. If we write x in the form 177*r+s where
0<=s<=176, then
171*x = 171*s +171*177*r
= 171*s + 30267*r
= 171*s + 30269*r - 2*r
Since 30269*r MOD 30269 = 0,
(171*x) MOD 30269 = (171*s - 2*r) MOD 30269
Therefore, the statement above can be written as
x := 171 * (x MOD 177) - 2 * (x DIV 177);
There is one problem, however. It can return a negative number if we do not
add the following:
IF x < 0 THEN x := x + 30269;
We now have a mechanism for making a simple generator for a 16 bit machine.
We choose a prime just slightly less than 32768 and then a multiplier (around
the square root of the prime). By using the above coding trick, the problem
of overflow can be eliminated with the equivalent of five or six additional
instrutions.
But such a simple generator on its own is inadequate. It will repeat itself
too soon, and its statistical properties are not near enough to those of
randomness. To improve it, we need a second technique that combines two or
more generators to make a new generator that is statistically superior to any
of its components. To understand this method, consider the random numbers x1
and x2 in the range 0<x<1. We will call the fractional part of x1+x2 the
combination of x1+x2. It is clear that if x1 and x2 are independent and
uniformly distributed, then the conbination of x1+x2 is also uniformly
distributed over the same range of values. Moreover, if each of x1 and x2 is
nearly random, the combination of x1+x2 will be much more nearly so than
either of them individually.
....
Combining the generators works correctly only if they are independent. In
fact, they cannot be completely independent since, if the three primes are p,
q, and r, the cycle lengths are p-1, q-1, and r-1, and these necessarily have
2 as a common factor. They can, however, be made nearly independent enough
for practical purposes if we make sure that they have no other common factor.
This means that the cycle length of the combined generator is
(p-1)(q-1)(r-1)/4 which in our case is about 6.95E12. This means that if 1000
numbers were calculated every second, it would not repeat itself for over 220
years.
*)
VAR
tmp : REAL;
BEGIN
x := 171 * (x MOD 177) - 2 * (x DIV 177);
IF x < 0 THEN x := x + 30269 END;
y := 172 * (y MOD 176) - 35 * (y DIV 176);
IF y < 0 THEN y := y + 30307 END;
z := 170 * (z MOD 178) - 63 * (z DIV 178);
IF z < 0 THEN z := z + 30323 END;
x1 := FLOAT(x);
y1 := FLOAT(y);
z1 := FLOAT(z);
tmp := x1/30269.0 + y1/30307.0 + z1/30323.0;
RETURN tmp - FLOAT(TRUNC(tmp));
END rand;
BEGIN
(* initialize seeds. For production runs, different values
[between 1 and 30000] should be used each time, preferably by
some automatic method such as from date and time readings if available. *)
(* I frequently just allocate a variable and then read from it to obtain
a ~random #. dag *)
x := 1;
y := 10000;
z := 3000;
END utils.
LONG LIVE vi!!
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
"Usual disclaimer about who's responsible for my mouth."
David Geiser
Hewlett-Packard
Logic System Division
P.O.B. 617
Colorado Springs, CO 80901-0617
UUCP: {hplabs,cbosgd}!hp-lsd!dag
CSNET: hp-lsd!dag@hp-labs.csnet
ARPA: hp-lsd!dag%hp-labs@csnet-relay.arpa
INTERNET: dag@hp-lsd.HP.COM
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=